Numerical approximation of Dtγ u with a controlled precision
Nehemie Nguimbous
May 2023

1

Introduction and Background

As indicated in the title of [2], the primary goal of this study is to numerically solve a one-dimensional
Gierer-Meinhardt model with sub-diffusion. The fractional operator Dtγ u from lemma (2.1) in [2] is an
inherent part of this system. It is essential to compute it with consistent precision to prevent any form of
contamination. Additionally, we chose t instead of y as in [2] as the dependent variable of Dtγ u to symbolize
its time dependency, as it is conceptually the case.
The fractional operator Dtγ u depends on two parameters: p and γ, with p being a parameter related to
the function u. As shown in the next section, numerically computing Dtγ u implies the integration of a spiketype fractional power function whose shape varies significantly, as shown in Figure (2) of Section (6). Since
that function is bell-shaped with a very steep slope, computing its integral requires a significant number
of subdivisions. Moreover, since the spikes are defined in an infinite domain, the number of subdivisions
required increases even further. Even though computing Dtγ u requires a high number of subdivisions in all
cases, their order of magnitude varies significantly depending on the shape of the spikes. For example, a
normal spike with a bell width of order O(ϵ) requires a number of subdivisions of order O(1/ϵ). Moreover,
the spike tail decays exponentially, which makes it easy to be captured in a relatively small interval. For
instance, a spike of width 0.1 centered at the origin on a domain of (-1, 1) will need 1/0.1 = 10 subdivisions
within (-0.05, 0.05) to be properly captured.
By contrast, an anomalous spike’s tail decays algebraically (i.e., considerably more slowly). It will then
need a very large and finely meshed interval to be properly captured. For example, an anomalous spike’s
tail could decay to 10−5 over a domain (-100, 100). The situation can worsen when dealing with interactions
between several spikes, where the domain increases even further, or in the case of chaotic behavior where
extremely fine time resolution is required. Furthermore, the number of subdivisions also increases with the
precision being required.
An important consideration arises: while employing numerical methods to approximate subdivisions,
these methods often provide error-bounded formulas based on subdivisions. It is logical to question why we
do not use such a formula to determine n based on the desired error. Let us explore this using Simpson’s
method as an example, which we will utilize further. According to [1, page 203], considering |f (4) | as the
continuous fourth derivative of a function f , and M as any upper bound for the values of |f (4) | on [a, b], then
the error |Ef | in the Simpson’s rule approximation of the integral of f from a to b satisfies the inequality:
|Ef | ≤

(b − a)5 M
· 4.
180
n

However, this approach faces two primary challenges. Firstly, it involves M , supposedly the maximum for the
fourth derivative of f . Yet, determining a maximum for a complex function like u(4) is unfeasible. Secondly,
1

Nehemie Nguimbous

an issue also arises with the term b − a. Operating within a theoretically infinite domain where a and b can
tend to infinity, the term b − a becomes an indeterminate form.
In summary, to make the best possible use of computer resources, it is essential to find a way to determine
the number of subdivisions required for controlling precision, depending on the shape of those spikes and
the domain in which they are defined. Answering this question will be the subject of our study.

2

Procedure
1. We split the main term of Dtγ u using two integrations by parts to simplify it and replace the original
improper integral by proper computable terms.
2. The results of these integrations is a sum of constants and a yet to be evaluated integral. We evaluate
the integral using the composite Simpson method. We selected the C programming language due to
its speed and efficiency, providing a significant advantage when dealing with memory-intensive and
power-demanding data.
3. We gradually increment the number of subdivisions used by our program to compute these integrals
until a precision of 10−10 is reached. These computations are performed for a predefined set of values
of the parameters p and γ.
4. We attempted to plot the correspondence between the number of subdivisions n and the value of the
integral obtained , but the resulting functions were neither smooth nor continuous. We then shifted
our focus toward the residual function.
5. The shapes of the residual functions are similar to that of the hyperbolic arctangent function, which
is then used as the fitting function.
6. Finally, we use the inverse functions of the fitting curves to approximate the number of subdivisions n
given a residual value r.

In summary, by implementing the above procedure, we expect to build a program capable of computing Dtγ u
for all values of p and γ with a precision of 10−10 .

3

Plan
1. Section 4 covers Step 1 in the procedure for regularizing Dtγ u.
2. Steps 2 and 3, starting from Section 6 up to Subsection 6.1.4, involve determining the number
of subdivisions n needed to compute I. This calculation is for (t, p, γ) values within the ranges:
{0.1, 1, 5} × {1.5, 2, 2.5, . . . , 4.5} × {0.1, 0.2, 0.3, . . . , 0.9}.
3. Steps 4, 5, and 6, spanning from Subsection 5.2 to 5.3.2, describe the fitting process and can be
visualized as follows:
We fit the correspondence
n(t,p,γ) → R(t,p,γ) (n),
using a variant of the hyperbolic-arctangent function, denoted as f . Consequently,
n(t,p,γ) → f(t,p,γ) (n) ∼ R(t,p,γ) .

Nehemie Nguimbous

To find the number of subdivisions based on a specific residual (the inverse path), we determine the
inverse f −1 of f so that
−1
f(t,p,γ)
(R) → n(t,p,γ) ,
where n(t,p,γ) and R(t,p,γ) represent the number of subdivisions and the corresponding residual values
for specific t, p, and γ values.
4. The process is now generalized for (t, p, γ) ∈ [0, 5] × [1.5, 4.5] × [0.1, 0.9] from Section 7 to 9 using a
series of linear and bilinear interpolations.
5. Section 10 involves the verification process where the accuracy of our results is assessed.

4

Regularization of Dtγ u

From section (2.1) in [2], we have
Dtγ u(t) = sign



dxi
dσ



1
Γ(−γ)

Z ∞





u(t) − u t + sign
0

dxi
dσ

  
γ+1
dxi 1
y
−
dy,
dσ y

(1)

with
γ
D(−t)
u(−t)

x′i >0

= Dtγ u(t)

x′i <0

.

From (1), the expression of Dtγ u is not directly computable because of the singularity present in the expression. Hence, we perform a double integration by parts and employ a Taylor expansion of u around a certain
point to get rid of the singularity.
For x′i < 0 and t > 0, we have
Z t∞
1
u(t) − u(t − y)
γ
Dt u(t) = −
dy.
Γ(−γ) 0
y γ+1
The first integration by parts
a(y) = u(t) − u(t − y),

db
1
= γ+1 ,
dy
y

leads to
Dtγ u(t) = −

1
Γ(−γ)

(

u(t) − u(t − t∞ )
1 u(t) − u(t − t∞ )
− lim
+
t→0 γ
γtγ∞
tγ

Z t∞
0

)
u′ (t − y)
dy ,
γy γ

u(t) − u(t − y)
= 0. The proof for this is straightforward since
tγ


u′′ (t)t2
′
+···+···
u(t) − u(t) + u (t)(−t) +
u(t) − u(t − y)
u′′ (t)t2−γ
2!
′
1−γ
∼
=
u
(t)t
−
+···
tγ
tγ
2

with lim

t→0

with 0 < γ < 1. We obtain

Nehemie Nguimbous

Dtγ u(t) = −

1
Γ(−γ)

(

u(t) − u(t − t∞ )
+
γtγ∞

Z t∞
0

)
u′ (t − y)
dy .
γy γ

By applying a second integration by parts on the last terms in the brackets,
a(y) = u′ (t − y),

db
1
= γ,
dy
y

we obtain
Z t∞
0

u′ (t − y)
t1−γ
1
∞
dy
=
−
u′ (t − t∞ ) −
γy γ
γ−1
γ(γ − 1)

Z t∞

u′′ (t − y)t1−γ dy.

0

This finally leads to
Dtγ u(t) = −

1
Γ(−γ)

(

t−γ
t1−γ
1
∞
∞
(u(t) − u(t − t∞ ) −
u′ (t − t∞ ) −
γ
γ(γ − 1)
γ(γ − 1)

Z t∞

)
′′

1−γ

u (t − y)t

dy .

(2)

0

We have
γΓ(−γ) = Γ(1 − γ) and γ(γ − 1)Γ(−γ) = Γ(2 − γ).
Additionally, from lemma (2.2) in [2],
u′′ (t − y) = u(t − y) − up (t − y),
which leads to
Z t∞

u′′ (t − y)t1−γ dy = −

Z t−t∞

0

u′′ (t)(t − y)1−γ dy.

0

We finally obtain
Dtγ u(t) = −

Z t−t∞


t−γ
t1−γ
1
∞
∞
′
u(t) − u(t − t∞ ) +
u (t − t∞ ) −
u′′ (y)(t − y)1−γ dy. (3)
Γ(1 − γ)
Γ(2 − γ)
Γ(2 − γ) y

Similarly, for x′i > 0, t < 0 and since u is even, we have
Z −t+t∞


t−γ
t1−γ
1
∞
u(−t) − u(−t + t∞ ) + ∞
u′ (t∞ − t) −
u′′ (y)(t + y)1−γ dy
Γ(1 − γ)
Γ(2 − γ)
Γ(2 − γ) −t
Z −t+t∞



1−γ
t−γ
t1−γ
1
∞
∞
′
u(t) − u(t − t∞ ) −
u (t − t∞ ) −
u′′ (−y) t − (−y)
dy
=
Γ(1 − γ)
Γ(2 − γ)
Γ(2 − γ) −t
Z t−t∞


t−γ
t1−γ
1
∞
=
u′′ (y)(t − y)1−γ dy
u(t) − u(t − t∞ ) − ∞
u′ (t − t∞ ) +
Γ(1 − γ)
Γ(2 − γ)
Γ(2 − γ) t

γ
D(−t)
u(−t) =

with 0 < γ < 1, and u defined as

u(t) =

p+1
(p − 1)t
sech2
2
2

1
 p−1

,

Nehemie Nguimbous

and verifying the following differential equation
u′′ − u + up = 0,
u′ (0) = 0,

u(0) > 0,

−∞ < t < ∞.

(4a)

and

(4b)

lim u = 0.

|t|→∞

Let us have
Z t−t∞


−γ
t1−γ
t1−γ
t∞
∞
∞
′
u(t − t∞ ) − u(t) , I 2 =
u′′ (y)(t − y)1−γ dy.
u (t − t∞ ), and I 3 =
I1 =
Γ(1 − γ)
Γ(2 − γ)
Γ(2 − γ) t
In both cases, Dtγ u is expressed as the sum of three terms, the first two being constants, and the third being
an integral. The first two terms are exact values, hence our approximation will focus solely on the integral.
Furthermore, it’s important to note that for both t > 0 and t < 0, the integrand as well as the integration
bounds remain the same. Therefore, the number of subdivisions required to compute the integral with the
desired
Z precision will be the same for both signs of t. Hence, we will concentrate solely on determining
t−t∞

I=

u′′ (y)(t − y)1−γ dy for a specific sign of t, in our case, t > 0.

t

5

Numerical approximation of I with a controlled precision

5.1

Structure of the project

Below is the repository tree of the project. At the root of the project, you have two main repositories, Code
and Data.

Figure 1: Tree diagram of the project repository.

Nehemie Nguimbous

• Repository path: \Code.
• Repository description: Contains all the code for the project, either in C or in Octave.
• Repository content: C, Octave.
– Repository path: \Code\C.
– Repository description: Contains all the C code for the project.
– Repository content: Du, functions, WriteFile.
∗ Repository path: \Code\C\functions.
∗ Repository description: Contains all the functions and their prototypes used to approximate
the number of subdivisions.
∗ Repository content: functions.h, functions.c.
· File name: functions.h.
· Contains all the functions’ prototypes used to approximate the number of subdivisions.
· File name: functions.c.
· File description: Contains the implementation of all the functions used in the approximation process.
· File content: sech, u, u1, u2, integrand, I_1, I_2, I_3, Du, simpson, findPosition_p,
findPosition_g, linearN_T, bilinearN_p_g, n_approx, n_approx_general, N_T, N_logT,
logN_T, logN_logT, even, reverse_n, join.
· File name: functions.h.
· File description: Contains the prototypes of all the functions used in the approximation
process.
∗ Repository path: \Code\C\Du.
∗ Repository description: Contains the C file used to reconstruct the curve of Dtγ u.
∗ Repository content: Du.c.
· File name: Du.c.
· File description: Contains the code used to reconstruct the curve of Dtγ u.
∗ Repository path: \Code\C\WriteFile.
∗ Repository description: Contains the C file used to implement the composite Simpson method
and to print its returned value in text files.
∗ Repository content: WriteFile.c.
· File name: WriteFile.c.
· File description: Contains the C code used to implement the composite Simpson method
and to print its returned value in text files.
– Repository path: \Code\Octave.
– Repository description: contains all the Octave code of the project.
– Repository content: Fit_Curve.
∗ Repository path: \Code\Octave\Fit_Curve.

Nehemie Nguimbous

∗ Repository description: Contains all the Octave files used to fit the residuals.
∗ Repository content: test1.m, fit_curve.m, fit_curvef.m
· file name: test1.m.
· file description: Main file containing the overall logic of the fitting process.
· file name: fit_curve.m.
· file description: Contains the implementation of the fit_curve function which use is
explained later in this document.
· file name: fit_curvef.m.
· file description: Contains the implementation of the fit_curvef function which use is
explained later in this document.
• Repository path: \Data.
• Repository description: Contains all the data used in our project.
• Repository content: P_G, Nmax_Nmin, R, C.
– Repository path: \Data\P_G.
– Repository description: Holds the files containing the returned value of the Simpson method for
t = 0.1, 1, 5.
– Repository content: t0p1, t1, t5.
– Repository path: \Data\C.
– Repository description: Holds the files containing the c(c1 , c2 ) parameters for t = 0.1, 1, 5.
– Repository content: t0p1, t1, t5.
– Repository path: \Data\Nmax_Nmin.
– Repository description: Holds the files containing the minimum and maximum number of subdivisions for t = 0.1, 1, 5.
– Repository content: t0p1, t1, t5.
– Repository path: \Data\R.
– Repository description: Holds the files containing the residuals for t = 0.1, 1, 5.
– Repository content: t0p1, t1, t5.

6

Finding the number of subdivisions required to approximate I
with a controlled precision

The integrand of I depends on three parameters: p, γ, and t∞ . Based on extensive numerical trials and as
demonstrated in the results section, the value of t∞ = 5 has proved to be sufficient to approach I with the
desired precision. In the subsequent sections of this document, when referring to discrete values of p and
γ, we respectively mean p ∈ {1.5, . . . , 4.5} and γ ∈ {0.1, . . . , 0.9}. Conversely, when discussing continuous

Nehemie Nguimbous

values of p and γ, we are referring to p ∈ [1.5, 4.5] and γ ∈ [0.1, 0.9]. As mentioned earlier, the shape
and gradient exhibit significant variation based on these parameters, influencing the number of subdivisions
required for computation. To streamline the process and avoid redundant calculations of I, we seek a method
to predict the necessary number of subdivisions and assess precision as we compute the integral. Below are
the curves of the integrand for some values of the parameters.

(a) u′′ (y)(t − y)1−γ for t = 5

(b) u′′ (y)(t − y)1−γ for t = 0.5

Figure 2: Plots of u′′ (t)(t − y)1−γ for p = 2, γ = 0.1 for t = 5(right) and t = 0.5(left). The x-axis represents
range of y-values and the y-axis the corresponding integrand values.

6.1

Computing I using the composite Simpson method

The following composite Simpson algorithm from [1, page 204] is used to compute the value I.
Z b
n/2−1
i
h X h
b−a
f (x)dx ≈
f (x2i ) + 4f (x2i+1 ) + f (x2i+2 ) with h =
.
3
n
a
i=0
6.1.1

Implementation of the Simpson method in C: located in functions.c

This function computes the integral of the function integrand using the composite Simpson method.
Parameter list:
• lower: lower bound of the integral.
• upper: upper bound of the integral.
• subInterval: number of sub-intervals.
• g: gamma parameter.
• A = pow((p+1.0)/2.0, 1.0/(p-1.0)).
• B = 2.0/(p-1.0).
• C = (p-1.0)/2.0.
Process:
• Calculate the step size.

Nehemie Nguimbous

• Sum the image by integrand of the lower and upper bounds of the integral first.
• Add the images by integrand of the values between the lower and upper bounds depending on whether
the discretization index is odd or even.
• Multiply the result by one-third of the stepSize.
1
2
3
4

double simpson ( double lower , double upper , long int subInterval , double p , double
g) {

5

double stepSize = ( upper - lower / subInterval ) ;

6
7

int i ;
double integration = integrand ( lower , p , g , y ) + integrand ( upper , p , g , y ) ;
for ( i = 1; i <= subInterval - 1; i ++) {
double k = lower + i * stepSize ;
if ( i % 2 == 0) {
integration = integration + 2 * integrand (k , p , g , y ) ;
} else {
integration = integration + 4 * integrand (k , p , g , y ) ;
}
}
integration = integration * stepSize / 3;
return integration ;

8
9
10
11
12
13
14
15
16
17
18
19
20
21

}
// implementation of the integrand

22
23
24

double integrand ( double y , double p , double g , double t ) {

25

return ( u (y , p ) - pow ( u (y , p ) , p ) ) * pow ( t - y , 1 - g ) ;

26
27
28

}

29
30

// implementation of u

31

double u ( double t , double p )
{
34
double A = pow ( ( p +1.0) /2.0 , 1.0/( p -1.0) ) ;
35
double B = 2.0/( p -1.0) ;
36
double C = (p -1.0) /2.0;
32
33

37

return ( A * pow ( sech ( C * t ) , B ) ) ;

38
39

}

Listing 1: Implementation of the composite Simpson method with its dependencies in C
Note: The number of subInterval used by the Simpson method must be even. Lines 48 to 50
from Listing 11 always ensure it is the case.

Nehemie Nguimbous

6.1.2

Printing the result of the Simpson method using the W riteF ile function: located in
W riteF ile.c

For each combination of t, p, and γ, the WriteFile function creates a file and prints, side by side, the return
values of the Simpson method with their corresponding number of subdivisions until the desired precision
is reached. Considering that we are dealing with very large values of n, and only the last values of n are of
interest, we double the increment value each time we run the first loop until we reach a precision of 5×10−10 .
Once that precision is achieved, new values of the increments are used until 1 × 10−10 is reached. The value
of the increments used in the second loop is returned by the getIncrement function and represents 5% of the
number of subdivisions required to reach the 5 × 10−10 precision in the first loop. This WriteFile function
is used to print the returned value of the Simpson method for discrete values of p and γ.
Parameter list:
• p: discrete values of p.
• g: discrete values of γ.
• t: ∈ {0.1, 1, 5}.
• n: max number of subdivisions allowed.
Process:
• Extract the integer and decimal part of p and g.
• Write or create a file for a specific combination of p and g where the returned value of the Simpson
method will be printed.
• Run the first loop and double the increment until 1 × 10−10 < error < 5 × 10−10 .
• Run the second loop with a custom increment until error1 < 10−10 .
• error and error1 have the same expression as Simpson(lower, upper, n) - Simpson(lower, upper, 2n).
Please note that ”Simpson(lower, upper, n)” represents the result of the Simpson method with the given
parameters lower, upper, and n.

1

void writeFile ( double p , double g , double t , long int n ) {

2
3
4
5
6
7
8

double
double
double
double
double
double

A = pow (( p + 1.0) / 2.0 , 1.0 / ( p - 1.0) ) ;
B = 2.0 / ( p - 1.0) ;
C = ( p - 1.0) / 2.0;
upper = t ;
tmax = 5;
lower = t - tmax ;

9
10
11
12
13

int int_g = 10 * g ;
double p1 = ( int ) p ;
double p2 = 10 * ( p - p1 ) ;
long int i = 0;

Nehemie Nguimbous

14
15
16
17
18
19
20

double g1 = ( int ) g ;
double g2 = 10 * ( g - g1 ) ;
long int increment = 1000;
long int next ;
double precision1 = 5e -10;
double precision2 = 1e -10;

21
22
23
24
25

long
char
FILE
FILE

int j ;
str_i [25];
* diag ;
* diag1 ;

26
27

sprintf ( str_i , " p %0.0 fp %0.0 f_g %0.0 fp %0.0 f . txt " , p1 , p2 , g1 , g2 ) ;

28
29
30
31

diag = fopen ( str_i , " a " ) ;
if (! diag ) {
printf ( " Failed to open diagonal dominance file .\ n " ) ;

32
33

}

34
35

for ( i = 1000; i < n ; i += i ) {

36

double value = simpson ( lower , upper , i , g , A , B , C ) ;

37
38

double value5 = simpson ( lower , upper , 2 * i , g , A , B , C ) ;

39
40

double error = fabs ( value - value5 ) ;

41
42

if ( error > precision1 ) {

43
44

fprintf ( diag , " % ld %0.15 f \ n " , i , value5 ) ;
} else if ( error <= precision1 && error > precision2 ) {
next = i ;
break ;
}

45
46
47
48
49
50
51

}

52
53

increment = getIncrement (t , p , g ) ; // Find the increment depending on t , p and
gamma

54
55

for ( j = next ; j < n ; j += increment ) {

56
57

double value1 = simpson ( lower , upper , j , g , A , B , C ) ;

58
59

double value10 = simpson ( lower , upper , 2 * j , g , A , B , C ) ;

60
61
62

double error1 = fabs ( value1 - value10 ) ;

Nehemie Nguimbous

if ( error1 < precision2 ) {

63
64

break ;
} else {

65
66
67

fprintf ( diag , " % ld %0.15 f \ n " , j , value10 ) ;

68
69

}

70
71

}

72
73

fclose ( diag ) ;

74
75

}

Listing 2: Implementation W riteF ile method in C

6.1.3

Calling the W riteF ile inside the main function: located in W riteF ile.c

This section shows how the WriteFile function is called in the main function for every values of p and γ.
We fixed the maximum number of subdivisions allowed to 7 billions.

1

int main () {

2

double p = 1.5;
double g = 0.1;
long int n = 7 e9 ;

3
4
5
6

for ( p = 1.5; p <= 4.5; p += 0.5) {

7
8

for ( g = 0.1; g < 0.9; g += 0.1) {

9
10

writeFile (p , g , t , n ) ;

11
12

}

13
14

}
return 0;

15
16
17

}

Listing 3: Calling the WriteFile method.

6.1.4

Maximum number of subdivisions required for t = 0.1 , 1 and 5

When the desired precision is reached, the last line printed by the WriteFile function holds the returned
value of the Simpson method with the right precision as well as the number of subdivisions required to reach
that precision. Below are tables displaying the maximum number of subdivisions for every t, p, and γ.

Nehemie Nguimbous

p/g
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9

1.5
35290
99900
216940
609778
1693204
4718250
16113048
76333764
383072000

2.0
51740
135800
344940
865778
2362408
7300300
28265572
132864908
730144000

2.5
61610
163900
433880
1101630
3386408
9970550
38166882
182371458
1018144000

3.0
73870
199800
522820
1377778
3721010
12464400
46087930
223956960
1028288000

3.5
80450
217750
567290
1495704
4390214
14066550
54550882
263749554
1028288000

4.0
87030
235700
611760
1731556
5059418
15668700
60491668
301374532
1028288000

4.5
93610
253650
689880
1849482
5394020
17270850
68412716
333058724
1028288000

Table 1: Number of subdivisions for t = 0.1

p/g
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9

1.5
28000
71900
175500
444800
1141700
3547200
11817400
52417200
260332000

2.0
26500
67550
166000
397600
1078730
3172400
10530500
46626150
232504000

2.5
11000
23900
60500
134800
318970
886800
2691450
12460850
44212000

3.0
22000
55700
130500
316800
826850
2335800
7956700
34400600
162934000

3.5
32500
83750
204000
515600
1393580
4109400
14626500
63653600
325868000

4.0
39500
103500
251500
633600
1716670
5233800
18487200
83600550
451094000

4.5
42500
111400
280000
700800
1905580
5983400
21061000
95826100
506750000

Table 2: Number of subdivisions for t = 1

p/g
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9

1.5
17100
42200
98200
234500
596400
1684000
5458000
21119000
101674000

2.0
11600
28000
63800
146500
364800
992000
3070000
11412000
51754000

2.5
9100
21700
48800
110500
268600
716000
2171000
7847000
34392000

3.0
8000
18800
41600
94000
227900
596000
1768000
6357000
27318000

3.5
7300
17300
38200
84500
202000
526000
1566000
5551000
23650000

4.0
6900
16100
35800
78500
187200
486000
1442000
5055000
21292000

4.5
6600
15500
34000
75500
178700
456000
1349000
4714000
19720000

Table 3: Number of subdivisions for t = 5

6.2

Fitting the residual function

Once the WriteFile function is executed, we are left with files containing two columns: the first representing
the number of subdivisions, and the other displaying the returned values of the Simpson method. Table 4
shows a sample file printed by the WriteFile function for t = 1, p = 4, and γ = 0.8. Since our objective is to
find the number of subdivisions required to reach a specific precision or error, it makes sense to establish a
relationship or correspondence between N and the error or residual (difference between the last value printed
by the WriteFile function and the other ones). Table 5 is an update of Table 4, now with residual values
as the second column.

Nehemie Nguimbous

N
75879150
76522600
77166050
77809500
78452950
79096400
79739850
80383300
81026750
81670200
82313650
82957100
83600550
84244000

I
-0.590274039200124
-0.590274039199605
-0.590274039198941
-0.590274039197867
-0.590274039197240
-0.590274039196180
-0.590274039195590
-0.590274039194982
-0.590274039194301
-0.590274039193123
-0.590274039192075
-0.590274039191671
-0.590274039191116
-0.590274039189668

Table 4: table4
6.2.1

N
75879150
76522600
77166050
77809500
78452950
79096400
79739850
80383300
81026750
81670200
82313650
82957100
83600550
84244000

R
1.2251e-11
1.0456e-11
9.9369e-12
9.2729e-12
8.1990e-12
7.5719e-12
6.5120e-12
5.9219e-12
5.3140e-12
4.6330e-12
3.4550e-12
2.4070e-12
2.0030e-12
1.4480e-12

Table 5: table5

Inside the test1.m file

The next step is to fit this correspondence with a usual function such that given any value of R, we could
find its corresponding N . During the fitting process of the residual function, we opted for the logarithm of
the residual instead of the actual residual function. In fact, the residual function is a power function, which
implies that its values are either extremely small or large over most of its domain. Then, attempting to get
insights from its plot has shown to be quite challenging. Also, the fit would be highly non-uniform. For
example, how do we ensure that the small values are fit as well as the large values? How do we measure what
is a good fit when the disparity in magnitude is so great? The log is the answer: it turns a very strongly
varying power function into quasi-straight lines. Straight lines are simpler to handle in terms of quantifying
their properties. If the original function is not a single power (which is our case), the log would not be a
straight line, but the magnitude disparity issue is resolved nonetheless. The other aspect of it all is that we
are indeed interested in the error magnitude, not its absolute value. We will always describe the residual as
a magnitude. Thus, the log treatment also makes sense conceptually, beyond being a technical convenience.
In the subsequent sections of the document, when we mention the residual function, we are actually referring
to its logarithm. Process:
• Extract the integer and decimal parts of p and g and convert them into strings.
• Open the file printed by the writeFile function for that specific combination of p and g.
• Load this file into the vector U .
• Compute the residual function.
• Fit the residual function.
• Plot the results.
In the context of the test1.m file, the goal is to establish a functional relationship between the number of
subdivisions N and the error or residual R for a given combination of p and γ. By fitting the plot of the

Nehemie Nguimbous

residual against the number of subdivisions with a continuous and easily invertible function, the program
can determine the number of subdivisions required to achieve a specific precision.
1
2
3

for p =1.5:0.5:4.5
for g = 0.1:0.1:0.9

4
5

% extract the integer and decimal part of p

6
7
8

p1 = fix ( p ) ;
p2 =10*( p - p1 ) ;

9

% extract the integer and decimal part of g
g1 = fix ( g ) ;
12
g2 =10*( g - g1 ) ;
10
11

13

% convert them into strings
sp1 = num2str ( p1 ) ;
16
sp2 = num2str ( p2 ) ;
17
sg1 = num2str ( g1 ) ;
18
sg2 = num2str ( g2 ) ;
14
15

19
20
21
22

%
open the file containing the return value of the WriteFile
method depending on a sepcific combination of p and g
%

23
24
25

file =[ ’C :/ Users / nguim / OneDrive / Documents / research_papers
/ Code / newdata / P_G / t5 / p ’ sp1 ’p ’ sp2 ’ _g ’ sg1 ’p ’ sg2 ’. txt ’ ];

26
27
28

% load the file as a 2 dimoentional vector
U = load ( file ) ;

29
30
31

N = U (: ,1) ; % the vector N contains the subdivisions
G = U (: ,2) ; % the vector G contains the actual return value the

WriteFile method

32
33

34
35
36

% R is the residual vector containing the difference between the last value of G
and the other ones
%
R = log10 ( abs ( G (1: end -1) -G ( end ) ) ) ;
N = N (1: end -1) ;

37

% c is a parameter vector that will minimizes the euclidian distance between
the Residual functions and their fitted curves
40 %
41
c0 =[ -1 -6]; % initial parameter of the fminunc function
42
c = fminunc ( @ ( c ) fit_curve (c ,R , N ) , c0 ) ;
38
39

43
44
45
46

[f , imin , imax ]= fit_curvef (c , N ) ;
N1 = N ;
N1 ([ imin imax ]) =[];

Nehemie Nguimbous

47
48

plot (N ,R , N1 ,f , ’ -- ’) ;
hold on

49
50
51

endfor
endfor

Listing 4: Fitting the residual function.

6.2.2

Inside the fit curvef.m file

The fit_curvef function is responsible for performing the fitting task of the residual using atanh.
Parameters list:
• c: a two-dimensional vector.
• n: the number of subdivisions.
Process:
• Get the minimum value of the vector n.
• Get the maximum value of the vector n.
• Remove these extreme values from the vector n.
• Determine the fitting curve using the parameter c.
1
2

function [f , imin , imax ]= fit_curvef (c , n )

3

[ minn , imin ]= min ( n ) ;
[ maxn , imax ]= max ( n ) ;
6 n ([ imin imax ]) =[];
7 n1 =2*( n - minn ) /( maxn - minn ) -1;
8 f = c (1) * atanh ( n1 ) + c (2) ;
4
5

Listing 5: Implementation of the f it curvef method .

6.2.3

Inside the fit curve.m file

The fit_curve function returns the Euclidean distance between the residuals and their fitting curves.
Parameters list:
• c: a two-dimensional vector returned by the fminunc function.
• u: the function being fitted; in our case, it represents the residual function.
• n: represents the number of subdivisions.
Process:
• Call the fit_curvef function and pass into it the vector c.

Nehemie Nguimbous

• fit_curvef returns the fitting curve f, as well as Imax and Imin, which respectively represent the
maximum and minimum values of n.
• Remove Imax and Imin from n.
• Determine the Euclidean distance between the fitting curve f and u.
1
2

function R = fit_curve (c ,u , n )

3

[f , imin , imax ]= fit_curvef (c , n ) ;
u ([ imin imax ]) =[];
6 R = sqrt ( sum (( f - u ) .^2) ) ;
4
5

Listing 6: Implementation of the f it curve method.

Figure 3: specific case for t = 5, p = 2, and γ = 0.8. The x-axis represents the number of subdivisions and
the y-axis the logarithm of the corresponding residuals. The solid curve represents the numerical results
while the dashed one represents the fitting curve.

6.3

Approximating N using the fitting curve

Once completing the fitting with the atanh function , its inverse can now be utilised to determine the
2 × (n − minn)
corresponding R for every value of n. From the fit_curvef.m file, we have n1 =
− 1 and
(maxn − minn)
tanh(f − c2 )
(n1 + 1) × (maxn − minn)
f = c1 ×atanh(n1 )+c2 . The inverse process leads to n1 =
and n =
+
c1
2
1. This expression is used to approximate n in the n_approx function present in the functions.c file.
6.3.1

Inside the n approxf unction: located in functions.c

This function computes and returns the number of subdivisions n given t, p, g, and the residual R.
Parameters list:

Nehemie Nguimbous

• t: upper bound of the integral.
• p: between 1.5 and 4.5 with an increment of 0.5.
• g: between 0.1 and 0.9 with an increment of 0.1.
• R: the residual value.
Process:
• Load the c.txt containing the different values of the c vector for a given t.
• Load the Nmax_Nmin.txt file containing the maximum and minimum number of subdivisions for a
given t.
• Compute n using the new formula.
1

double n_approx ( double t , double p , double g ) {

2
3
4
5
6
7
8
9
10
11
12

double n1 , n , pos_p , pos_g ;
FILE * fp1 ;
FILE * fp2 ;
FILE * fp3 ;
char file_name1 [25];
char file_name2 [33];
char file_name3 [25];
float c1 = 0 , c2 = 0 , nMax = 0 , nMin = 0 , R ;
char * param ;
int i = 0 , location = 0;

13
14
15

pos_p = 1 + ( p - 1.5) / 0.5;
pos_g = g / 0.1;

16
17

// find the location of c1 and c2 in the C . txt file as well as the location of
Nmax and Nmin in the N . txt file knowing p and g

18
19

location = ( int ) (( pos_p - 1) * 9 + pos_g ) ;

20
21

// assign a value to param based on the value of t

22
23
24
25
26
27
28
29

if ( t == 0.1) {
param = " 0 p1 " ;
} else if ( t == 1) {
param = " 1 " ;
} else {
param = " 5 " ;
}

30
31

// create a path to the C . txt , N . txt , and R . txt files using the param variable

32
33
34

sprintf ( file_name1 , " ../../ Data / C / t % s / C . txt " , param ) ;
sprintf ( file_name2 , " ../../ Data / Nmax_Nmin / t % s / N . txt " , param ) ;

Nehemie Nguimbous

sprintf ( file_name3 , " ../../ Data / R / t % s / R . txt " , param ) ;

35
36

// open N . txt , C . txt , and R . txt

37
38

fp1 = fopen ( file_name1 , " r " ) ;
fp2 = fopen ( file_name2 , " r " ) ;
fp3 = fopen ( file_name3 , " r " ) ;

39
40
41
42

if (! fp2 ) {
printf ( " fail_2 " ) ;
}

43
44
45
46

if (! fp1 ) {
printf ( " fail_1 " ) ;
}
if (! fp3 ) {
printf ( " fail_3 " ) ;
}

47
48
49
50
51
52
53

for ( i = 0; i < location ; i ++) {
fscanf ( fp1 , " % f % f " , & c1 , & c2 ) ; // find c1 and c2 using the location
variable
fscanf ( fp2 , " % f % f " , & nMin , & nMax ) ; // find nMin and nMax using the
location variable
fscanf ( fp3 , " % f " , & R ) ; // find R using the location variable

54
55

56

57
58

}

59

;

60

R = R - 1.5;

61
62

// close the files

63
64

fclose ( fp1 ) ;
fclose ( fp2 ) ;

65
66
67

n1 = tanh (( R - c2 ) / c1 ) ;
n = (( n1 + 1) * ( nMax - nMin ) ) / 2 + 1; // find n

68
69
70

return n ;

71
72
73

}

Listing 7: Implementation of the n approx method.

6.3.2

Approximated number of subdivisions required for t = 0.1, 1, and 5

Below are the approximate number of subdivisions returned by the n_approx_general function for t = 0.1,
t = 1 and t = 5.

Nehemie Nguimbous

g/p
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9

1.5
35065
99002
214565
602747
1674913
1674913
4662043
15930379
76155108

2.0
51414
134710
342453
858551
2341552
7241978
28071526
132684487
727802772

2.5
61242
162937
430919
1093366
3365120
9904136
37953058
182123946
1015373018

3
73570
198630
519457
1370345
3698389
12403810
45866233
223677191
1025834015

3.5
80108
216504
563799
1487766
4365537
13999256
54342976
263586831
1025841827

4
86651
234382
608143
1722571
5033725
15597980
60260601
301086078
1025918580

4.5
93201
252262
686814
1840192
5366720
17196034
68207443
332831576
1025859852

Table 6: Approximated number of subdivisions for t = 0.1

g/p
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9

27873
27873
71597
174554
442094
1134337
1134337
3527550
11753526
52352342

26378
26378
67150
165106
395096
1071579
3153872
10468585
46563914
231613175

10927
10927
23662
59775
132836
313754
871089
2640180
12411139
43366292

21895
21895
55344
129525
314216
820575
2317025
7900759
34340174
162079821

32362
32362
83377
202917
512656
1385528
4088223
14564452
63593768
324972385

39371
39371
103044
250255
630341
1708829
5210846
18417018
83514848
449967670

42358
42358
110919
278688
697983
1897254
5958598
20991875
95754535
505825475

Table 7: Approximated number of subdivisions for t = 0.1

g/p
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9

17083
17083
42163
98128
234320
595947
595947
1682968
5453895
21115912

11585
11585
27968
63737
146338
364389
990878
3066607
11407642
51737121

9085
9085
21670
48741
110344
268200
714928
2168010
7844868
34370523

7986
7986
18772
41544
93857
227526
595009
1765265
6353968
27302158

7286
7286
17273
38147
84362
201649
525021
1563165
5548342
23631082

6886
6886
16074
35748
78364
186865
485046
1439166
5052182
21275601

6586
6586
15474
33949
75366
178355
455084
1346239
4711303
19703359

Table 8: Approximated number of subdivisions for t = 5

7

Number of subdivisions required to compute I for continuous
values of t and discrete values of p and γ

Now that we have a program capable of predicting the number of subdivisions n required to compute I with
controlled precision, our goal is to expand it for continuous values of t, discrete values of p and γ. For that,
we implemented four different variants of linear interpolation formulas and selected the most conservative
one, which returns the highest value of n. For t ∈ [t1 , t2 ] their general formulas are

Nehemie Nguimbous

1. nt = nt1 +

(t−t1 )(nt2 −nt1 )
,
(t2 −t1 )

2. nt = nt1 +

(log(t)−log(t1 ))(nt2 −nt1 )
,
(log(t2 )−log(t1 ))

3. log(nt ) = log(nt1 ) +

(t−t1 )(log(nt2 )−log(nt1 ))
,
(t2 −t1 )

4. log(nt ) = log(nt1 ) +

(log(t)−log(t1 ))(log(nt2 )−log(nt1 ))
,
(log(t2 )−log(t1 ))

with log referring to the log10 function.

7.1

Implementation in C: located in the functions.c file.

Below are their respective implementation in C.
1
2

int N_T ( double t , double t1 , double t2 , int n1 , int n2 )
{
5
int n = ( int ) ( n1 + ( t - t1 ) *( n2 - n1 ) /( t2 - t1 ) ) ;
6
return n ;
7 }
3
4

8

int N_logT ( double t , double t1 , double t2 , int n1 , int n2 )
{
11
int n = ( int ) ( n1 + ( log10 ( t ) - log10 ( t1 ) ) *( n2 - n1 ) /( log10 ( t2 ) - log10 ( t1 ) ) ) ;
12
return n ;
13 }
9

10

14

int logN_T ( double t , double t1 , double t2 , int n1 , int n2 )
{
17
double logn = log10 ( n1 ) + ( t - t1 ) *( log10 ( n2 ) - log10 ( n1 ) ) /( t2 - t1 ) ;
18
return ( int ) ( pow (10 , logn ) ) ;
19 }
15
16

20

int logN_logT ( double t , double t1 , double t2 , int n1 , int n2 )
{
23
double logn = log10 ( n1 ) + ( log10 ( t ) - log10 ( t1 ) ) *( log10 ( n2 ) - log10 ( n1 ) ) /
24
( log10 ( t2 ) - log10 ( t1 ) ) ;
25
return ( int ) ( pow (10 , logn ) ) ;
26 }
21
22

Listing 8: Implementation of the previous interpolation functions in C.

7.2

Determining the most conservative linear interpolation method: inside the
linearNT function, located in the functions.c file

Since we have four variants of linear interpolation, the task of linearN_T is to find the most conservative
one. n1 and n2 are respectively the number of subdivisions required for t1 and t2 with t1 <= t <= t2.
The linearN_T function evaluates the number of subdivisions required for a given range of t values and
then selects the most conservative interpolation method among the four available variants.

Nehemie Nguimbous

int linearN_T ( double t , double t1 , double t2 , int n1 , int n2 )
{
3
int
temp1 = N_T (t , t1 , t2 , n1 , n2 ) ;
4
int temp2 = N_logT (t , t1 , t2 , n1 , n2 ) ;
5
int
temp3 = logN_T (t , t1 , t2 , n1 , n2 ) ;
6
int
temp4 = logN_logT (t , t1 , t2 , n1 , n2 ) ;
1
2

7

int n = temp1 ;

8
9

if ( n < temp2 )
{
n = temp2 ;
}

10
11
12
13
14

if ( n < temp3 )
{
n = temp3 ;
}

15
16
17
18
19

if ( n < temp4 )
{
n = temp4 ;
}

20
21
22
23
24

return n ;

25
26

}

Listing 9: Code used to compare the different interpolation methods.

8

Number of subdivisions required to compute I with a control
precision for continuous values of t, p and and γ

In the last section, we were able to reconstruct the curve of Dtγ u for p = 2 and discrete values of γ. Moreover,
we can do the same for all discrete values of p as well. However, we face a problem for continuous values of
p and γ. For example, for p = 2.38 and γ = 0.46, we need to employ a different approach. Again, we use
interpolation, but this time it is bi-linear instead of linear. In fact, for a given (p, γ) such that p1 ≤ p ≤ p2
and γ1 ≤ γ ≤ γ2 , we want to find its corresponding number of subdivisions n. Below are the bi-linear
interpolation formula and its implementation in our code.

8.1

Bi-linear interpolation formula

2
n = γγ−γ
1 −γ2

8.2



np1 ,γ1 (p−p2 )
n
(p−p1 )
1
+ p2 ,γ
p1 −p2
p2 −p1



γ1
+ γγ2−−γ
1



np1 ,γ2 (p−p2 )
n
(p−p1 )
2
+ p2 ,γ
p1 −p2
p2 −p1



.

Implementation in C: inside the billinearN p g function: inside the functions.c file.

Nehemie Nguimbous

1
2

int bilinearN_p_g ( double p , double g , double t ) {

3

int n ;
double p1 , p2 , g1 , g2 ;
double n_p1g1 , n_p1g2 , n_p2g1 , n_p2g2 ;

4
5
6
7

findPosition_p (p , & p1 , & p2 ) ; // find p1
findPosition_g (g , & g1 , & g2 ) ; // find p2

8
9
10

n_p1g1
n_p1g2
n_p2g1
n_p2g2

11
12
13
14

=
=
=
=

n_approx (t ,
n_approx (t ,
n_approx (t ,
n_approx (t ,

p1 ,
p1 ,
p2 ,
p2 ,

g1 ) ;
g2 ) ;
g1 ) ;
g2 ) ;

//
//
//
//

find
find
find
find

n
n
n
n

for
for
for
for

p = p1
p = p1
p = p2
p = p2

and
and
and
and

g = g1
g = g2
g = g1
g = g2

15

// find n

16
17

n = ( int ) (( g - g2 ) / ( g1 - g2 ) * ( n_p1g1 * ( p - p2 ) / ( p1 - p2 ) + n_p2g1 * ( p p1 ) / ( p2 - p1 ) )
+ ( g - g1 ) / ( g2 - g1 ) * ( n_p1g2 * ( p - p2 ) / ( p1 - p2 ) + n_p2g2 * ( p - p1 ) / (
p2 - p1 ) ) ) ;

18

19

20

return n ;

21
22

}

Listing 10: Bi-linear interpolation of n in C.

9

Summary of the interpolation process

On one hand, n_approx allows us to predict the number of subdivisions required to approximate I for
discrete values of t, p, and γ. On the other hand, n_approx_general performs the same task but over
their continuous range. It achieves this by employing both linear and bi-linear interpolation methods. The
process to find the number of subdivisions required for an unknown triplet (t, p, γ) involves finding two
triplets (t1 , p1 , γ1 ) and (t2 , p2 , γ2 ) such that t1 ≤ t ≤ t2 , p1 ≤ p ≤ p2 , and γ1 ≤ γ ≤ γ2 . Moreover, it finds
n1 required for (t1 , p1 , γ1 ) and n2 required for (t2 , p2 , γ2 ), and then computes n required for (t, p, γ) using
interpolation. This comprehensive approach allows us to efficiently estimate the number of subdivisions
needed to achieve the desired precision for the integral over a wider range of input values.

Nehemie Nguimbous

9.1

Algorithm of the n approx general function

Algorithm 1: n approx general Algorithm
Result: Find n for continuous values of y, p and γ
;
if y ∈ {0.1, 1, 5} then
n = bilinear(t, p, γ) ;
end
else
find (t1 , t2 ) ∈ {0.1, 1, 5} × {0.1, 1, 5} such that t1 ≤ t ≤ t2 ;
n1 = bilinear(t1, p, γ) ;
n2 = bilinear(t2, p, γ) ;
n = linear(t, t1, n1, t2, n2) ;
end

9.2

1

Implementation in C of the n approx general function: : inside the functions.c file

int n_approx_general ( double t , double p , double g ) {

2
3

double t1 , t2 ;

4
5
6
7
8
9

double p1 = ( int ) p ; // extract the integer part of p
double p2 = 10 * ( p - p1 ) ; // extract the decimal part of p
long int i = 0;
double g1 = ( int ) g ; // extract the integer part of g
double g2 = 10 * ( g - g1 ) ; // extract the decimal part of p

10
11

int n , n1 , n2 ;

12
13

// if t is not in the set {0 , 0.1 , 1 , 5} , find t1 and t2 such that t1 <= t <= t2

14
15
16
17

if ( t != 0.1 && t != 5 && t != 1 && t != 0) {
t1 = 0.1;
t2 = 10;

18
19
20
21
22
23
24
25
26
27
28
29

if ( t < 1) {
t1 = 0.1;
t2 = 1;
} else if (1 < t < 5) {
t1 = 1;
t2 = 5;
} else {
t1 = 0.1;
t2 = 5;
}

Nehemie Nguimbous

n1 = bilinearN_p_g (p , g , t1 ) ; // find n for p , g and t1 using bilinear
interpolation
n2 = bilinearN_p_g (p , g , t2 ) ; // find n for p , g and t2 using bilinear
interpolation

30

31

32

n = linearN_T (t , t1 , t2 , n1 , n2 ) ; // find n knowing n1 and n2 using linear
interpolation

33

34

} else if ( t == 0) // find n for t =0
{
t = 0.1;
n = bilinearN_p_g (p , g , t ) + 1;

35
36
37
38
39

} else {

40
41

// if t is in the set {0.1 , 1 , 5} use bilinear interpolation to compute it
directly
n = bilinearN_p_g (p , g , t ) ;

42

43
44

}
// verify that n is even , because the composite Simpson method only works with
even values of n .

45
46

47

if ( n % 2 != 0) {
n += 1;
}

48
49
50
51

return n ;

52
53

}

Listing 11: Implementation of the n approx general function

10

Assessment of the accuracy of the results

As a primary test for the quality of the results, we aim to reproduce the curve of Dtγ u on [−t∞ , t∞ ], and
compare its shape with the one from [1, page 12]. The newly computed curve is at the left and the former
one at the right.

Nehemie Nguimbous

error
-4.744730e-008
-5.048329e-009
6.482569e-008
-4.435881e-009
2.060108e-007
-3.872901e-009
3.676208e-007
-3.364300e-009
-3.129406e-009
6.138696e-007
-2.703101e-009
7.505957e-007
-2.328432e-009
8.405479e-007
-2.001432e-009

y
3.6
3.7
3.8
3.9
4.0
4.1
4.2
4.3
4.4
4.5
4.6
4.7
4.8
4.9
5.0

error
-1.784137e-010
-1.721758e-010
-1.664409e-010
-1.605218e-010
-1.547047e-010
-1.494032e-010
-1.443059e-010
-1.395319e-010
-1.351115e-010
-1.307046e-010
-1.269774e-010
-1.235263e-010
-1.206676e-010
-1.180839e-010
-1.161406e-010

y
3.6
3.7
3.8
3.9
4.0
4.1
4.2
4.3
4.4
4.5
4.6
4.7
4.8
4.9
5.0

Table 9: Error file for γ = 0.7

Table 10: New error file for γ = 0.7

(a) Plot of the reconstructed Dtγ u.

(b) Plot of the reconstructed Dtγ u .

Figure 4: Comparison between the former and reconstructed Plot of Dtγ u for sign
The x-axis represents range of t-values and the y-axis the associated Dtγ u.



dxi
dσ


<0

and p = 2.

The reproduction of Dtγ u was just the visual part of the assessment process. In fact, we still need to verify
the precision of Dtγ u values. To achieve this, we compute the absolute difference between Dtγ u for N = n,
N = 2n and p = 2. Initially, the error variable used in the stopping condition of the writeFile function
was defined as error = abs(Simpson(lower, upper, n) - Simpson(lower, upper, n + 10)). Due to
the non-monotonic and the uncontrolled character of the errors’ magnitude depicted in table 9, we adjusted
the error variable to error = abs(Simpson(lower, upper, n) - Simpson(lower, upper, 2*n)). This
modification yielded satisfactory results, as illustrated in table 10.

10.1

Implementation of Dtγ u in C: Inside the Du.c file

Below is the implementation of the Dtγ u function.

Nehemie Nguimbous

10.1.1

Implementation of Dtγ u

The function Du calculates the value of Dtγ u for −5 ≤ t ≤ 5.
Parameters list:
• t: the function’s variable, between -5 and 5.
• nt: the number of subdivisions required to compute I3 .
Process:
• Calculate the first operand of Dtγ u and store its value in I1.
• Calculate the second operand of Dtγ u and store its value in I2.
• Calculate the third operand of Dtγ u, the one with the integral, and store its value in I3.
1

double Du ( double t , long int nt , double tmax , double p , double g ) {

2

double I1 = I_1 (t , tmax , p , g ) ;
double I2 = I_2 (t , tmax , p , g ) ;
double I3 = I_3 (t , tmax , nt , p , g ) ;

3
4
5
6

return I1 + I2 - I3 ;

7
8
9

}

Listing 12: Implementation of the Dtγ u function in C

10.1.2

Implementation of I1 , I2 , and I3

I1 , I2 and I3 are used to implement respectively the first, second and third operands of Dtγ u
1

// Implementation of I_1

2
3

double I_1 ( double t , double tmax , double p , double g ) {

4

return ( pow ( tmax , -g ) / ( tgamma (1 - g ) ) ) * ( u (t - tmax , p ) - u (t , p ) ) ;

5
6
7

}

8
9

// Implementation of I_2

10
11

double I_2 ( double t , double tmax , double p , double g ) {

12

return ( pow ( tmax , 1 - g ) / tgamma (2 - g ) ) * ( u1 (t - tmax , p ) ) ;

13
14
15

}

16
17

// Implementation of I_3

Nehemie Nguimbous

18
19

double I_3 ( double t , double tmax , long int nt , double p , double g ) {

20

return (1 / tgamma (2 - g ) ) * simpson (t , t - tmax , nt , p , g ) ;

21
22
23

}

Listing 13: Implementation of I1 , I2 and I3 in C.

References
[1] Richard L Burden, Douglas J Faires, and Annette M Burden. Numerical Analysis. Brooks Cole, 10
edition, 2015.
[2] Y. Nec and M. J. Ward. Dynamics and stability of spike-type solutions to a one dimensional gierermeinhardt model with sub-diffusion. Physica D, 241:947—963, 2012.