#15 (03/27/2024)

Numerical differentiation

The only occasion where numerical differentiation is needed is when data are given numerically only, i.e.

x	f (x)
1	1
2	8
3	27
4	64

Forward difference

f (x + h)

f (x) + h f ′(x) +

h²

f " (x) +

h³

f "′(x)+…

(1)

≈

f (x) + h f ′(x),

(2)

hence

f ′(x) ≈

f (x + h) − f (x)

(3)

Backward difference

f ( x − h )

f (x) − h f ′(x) +

h²

f "(x) −

h³

f "′(x) + …

≈

f (x) − h f ′(x),

(4)

hence

f ′(x) ≈

f (x) − f (x − h)

(5)

Central difference

f (x + h)

f (x) + h f ′(x) +

h²

f "(x) +

h³

f "′(x)+…

(6)

f (x − h)

f (x) − h f ′(x) +

h²

f "(x) −

h³

f "′(x)+…

(7)

from which

f ( x + h ) − f ( x − h ) = 2 h f ′(x) + 2

h³

f "′(x)+ …

(8)

hence

f ′(x) ≈

f ( x + h )− f ( x − h )

(9)

It is seen from the above that the central difference should yield better approximation as the truncation error is in the order of h² while the truncation error for the forward and backward difference methods is in the order of h.

Forward difference:

f ′(x) ≈ f ( x + h )−f (x)
h
+ O(h).
Backward difference:

f ′(x) ≈ f (x)−f ( x − h )
h
+ O(h).
Central difference:

f ′(x) ≈ f ( x + h )−f ( x − h )
2h
+ O(h²).

Second order derivative

f (x + h)

f (x)+h f ′(x) +

h²

f "(x) +

h³

f "′(x)+…

(10)

f (x − h)

f (x)−h f ′(x) +

h²

f "(x) −

h³

f "′(x)+…

(11)

Adding the above two yields

f (x + h)+f (x − h) = 2 f (x) + h² f "(x) + …

(12)

f "(x) ≈

f (x + h) + f (x − h) − 2 f (x)

h²

(13)

Example

The program below is to numerically differentiate a function given by the table below (numerical values of sin x for x = 0.0 ∼ 1.0). for x = 0 to x = 1 with the step size of h = 0.1 excluding at x=0 and x=1.0 where the central difference scheme does not work.

x	0.0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1.0
f(x)	0	0.0998	0.1986	0.2955	0.3894	0.4794	0.5646	0.6442	0.7173	0.7833	0.8414


#include <stdio.h>
#define N 11

int main()
{
double y[N]={0, 0.0998, 0.1986, 0.2955, 0.3894, 0.4794, 0.5646, 0.6442, 0.7173, 0.7833, 0.8414};
double central[N], h=0.1;
int i;

/* Excluding x=0 and x=1.0 */
for (i=1;i<N-1;i++) central[i]= (y[i+1]-y[i-1])/(2*h);

printf ("   x    Central \n---------------------------\n");

for (i=1;i<N-1;i++)
 printf ("%lf %lf\n", i*h, central[i]);

return 0;
}

What if you want to obtain f′(x) at x = a with the same accuracy as the central difference ?

Obviously you cannot use the central difference scheme at x=a. Using the forward difference scheme, one can get

f′(x) ∼

f(x+h) − f(x)

which is a poor approximation at best.

Consider

f(x + h)

f(x) + h f′(x) +

h²

f"(x) + …

(14)

f(x + 2 h)

f(x) + 2 h f′(x) +

4 h²

f"(x) + …

(15)

The h² term can be eliminated if one subtracts Eq.(15) from 4 times Eq.(14) as

4 f(x + h) − f(x + 2 h) = 3 f(x) + 2 h f′(x) + (higher order) …

(16)

from which one obtains

f′(x) ∼

4 f(x+h) − f(x+2h) − 3 f(x)

2 h

(17)

Eq.(17) has the same order of accuracy as the central difference scheme.

So the program above can be modified as:


#include <stdio.h>
#define N 11

int main()
{
double y[N]={0, 0.0998, 0.1986, 0.2955, 0.3894, 0.4794, 0.5646, 0.6442, 0.7173, 0.7833, 0.8414};
double central[N], h=0.1;
int i;

for (i=1;i<N-1;i++) central[i]= (y[i+1]-y[i-1])/(2*h);

central[0]=(4*y[1]-y[2]-3*y[0])/(2*h);

printf ("   x    Central \n---------------------------\n");

for (i=0;i<N-1;i++)
 printf ("%lf %lf\n", i*h, central[i]);

return 0;
}

Integration review

f (x)	∫f (x) dx
x^α (α ≠ −1)	x^{α+ 1}/(α+ 1)
1/x	ln x
e^x	e^x
1/(x² + a²)	arctan x
sin k x	− cos k x/k
cos k x	sin k x/k

Some exercise:

⌠
⌡

x² − 1

⌠
⌡

π

0

x sin x dx

⌠
⌡

∞

0

x² + 1

Numerical Integration

Rectangular rule:

∼

h ×f₀ + h ×f₁ + h ×f₂ + …+ h ×f_n−1

(18)

h ( f₀ + f₁ + …+ f_n−1 ).

(19)

Sample program Implement

⌠
⌡

1

0

1+x²

dx = 4 arctan x

⎢
⎢

1
0

= π.

to approximate the value of π.


/* Rectangular rule */
#include <stdio.h>

double f(double x)
 {return 4.0/(1.0+x*x) ; }

int main()
{
 int i, n;
 double a=0.0, b=1.0, h, s=0.0 , x ;

 printf("Number of partitions = ");
 scanf("%d", &n) ;

 h = (b-a)/n ;

 for (i= 0;i<n;i++) s = s + f(a + i*h) ;

  s=s*h ;

  printf("Result =%lf\n", s) ;
  return 0;
}

Trapezoidal rule:

∼

( f₀ + f₁ ) +

( f₁ + f₂ )+

( f₂ + f₃ )+ …+

( f_n−1 + f_n )

(f₀ + 2 f₁ + 2f₂ +…+2f_n−1+f_n).

(20)

Implementation

I ∼

(f₀+f_n) + h ×(f₁+f₂+f₃+…+ f_n−1).

/* Trapezoidal rule */
#include <stdio.h>

double f(double x)
 {return 4.0/(1.0+x*x);}

int main()
{
 int i, n ;
 double a=0.0, b=1.0 , h, s=0.0, x;

 printf("Enter number of partitions = ");
 scanf("%d", &n) ;

 h = (b-a)/n ;

 for (i=1;i<=n-1;i++) s = s + f(a + i*h);

  s=h/2*(f(a)+f(b))+ h* s;

  printf("%20.12f\n", s) ;
  return 0;
}

Simpson's rule:

Approximate the curve that passes (−h, f(−h)), (0, f(0)) and (h, f(h)) by

y = a x²+ b x + c.

(21)

f(−h)

a h² − b h + c,

(22)

f(0)

(23)

f(h)

a h²+ b h + c,

(24)

from which

f(h)+ f(−h) − 2 f(0)

2h²

(25)

f(h)−f(−h)

(26)

f(0).

(27)

Hence,

⌠
⌡

h

−h

(a x² + b x+ c) dx

⌠
⌡

h

−h

⎛
⎝

f(h)+ f(−h) − 2 f(0)

2h²

x² +

f(h)−f(−h)

x+ f(0)

⎞
⎠

…

⎛
⎝

f(−h)+ 4 f(0)+ f(h)

⎞
⎠

(28)

This is like a weighted average of f(−h), f(0) and f(h) as

S ≈

f(−h)+ 4 f(0)+ f(h)

×(2h).

(29)

By adding this for the interval [a,b], one obtains

∼

⎛
⎝

(f₀ + 4f₁ + f₂ )+(f₂ + 4f₃ + f₄ )+(f₄ + 4f₅ + f₆ )+…+(f_2n−2+ 4 f_2n−1+f_2n )

⎞
⎠

∼

⎛
⎝

f₀ +4 f₁+2 f₂ + 4f₃ + 2f₄ + …+ 2 f_2n−2 + 4 f_2n−1 + f_2n

⎞
⎠

(30)

where

h =

b−a

Note that the number of partitioning for the Simpson rule must be an even number (=2n).

For the coding purpose, it is convenient to write Eq.(30) as

⎛
⎝

f₀+f_2n

⎞
⎠

×4

⎛
⎝

f₁ + f₃ + f₅ + …+f_2n−1

⎞
⎠

×2

⎛
⎝

f₂+f₄+ f₆ + …+ f_2n−2

⎞
⎠

(31)

Implementation


/* Simpson's rule */
#include <stdio.h>
#include <math.h>

double f(double x)
 {return 4.0/(1.0+x*x);}

int main()
{
 int i, n ;
 double a=0.0, b=1.0 , h, s1=0.0, s2=0.0, s3=0.0, x;

 printf("Enter number of partitions  = ");
 scanf("%d", &n) ;

  h = (b-a)/(2.0*n) ;

  s1 = (f(a)+ f(b));

 for (i=1; i<2*n; i=i+2) s2 = s2 + f(a + i*h);

 for (i=2; i<2*n; i=i+2) s3 = s3 + f(a + i*h);

    printf("%20.12lf\n", (h/3.0)*(s1+ 4.0*s2 + 2.0*s3)) ;
 return 0;
}

According to error analysis, the truncation error for each rule is as follows:

Rectangular rule

I ∼ A + f′(ξ) h
Trapezoidal rule

I ∼ A + f"(ξ) h²
Simpson's rule

I ∼ A + f"′(ξ) h³

where h is the step size and ξ is a value within the interval.

Computation of π

⌠
⌡

1

0

1+x²

dx =

⌠
⌡

1

0

√

1−x²

dx =

It is important to know that the accuracy of numerical integration generally depends on the derivatives of the function, f(x), to be integratated. Although the two examples above yield approximation to π, the convergence speed (and accuracy) varies greatly. For further information, read the following link:
Error analysis

File translated from T_EX by T_TH, version 4.03.
On 26 Mar 2024, 12:38.