#16 (03/24/2025)

Numerical Integration

Some of important integrals:


⌠
⌡ 		xα dx
xα+ 1
α+ 1
⌠
⌡ 		1
x
 dx
ln x
⌠
⌡ 		dx
1 + x2
arctan x
⌠
⌡ 		dx

√
a2 − x2
arcsin (x/a)
⌠
⌡ 		f(x) g′(x) dx
f(x)g(x) −⌠
⌡ 		f′(x) g(x) dx
⌠
⌡ 		f(x) dx
x f(x) − ⌠
⌡ 		x f′(x) dx

Exercises

Evaluate the following integrals:

(a) ⌠
⌡ 		dx
ex + e−x
(b)
⌠
⌡
√
 
a2 − x2
 
 dx
(c) ⌠
⌡ 		eax cos b x dx
(d)
⌠
⌡
√
 
x2 + 1
 
 dx
It is more practical these days to carry out integrations via computer.
Enter the following line 
Sqrt[x^2 + 1]
exactly (case sensitive, square bracket) in the first box and x (lower case) in the second box, then press the query button.
Integrate with respect to .

Numerical integrals:

numint.gif
  • Rectangular rule:

    I
     ∼ 
    h ×f0 + h ×f1 + h ×f2 + …+ h ×fn−1
    		(1)
    =
    h ( f0 + f1 + …+ fn−1 ).
    		(2)
    Sample program Implement
    4 ⌠
    ⌡     		1

    0 
    		 1
    1+x2
    		 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:
    trap.gif

    I
     ∼ 
    h
    2
    		 ( f0 + f1 ) + h
    2
    		 ( f1 + f2 )+ h
    2
    		 ( f2 + f3 )+ …+ h
    2
    		 ( fn−1 + fn )
    =
    h
    2
    		 (f0 + 2 f1 + 2f2 +…+2fn−1+fn).
    		(3)
    Implementation

    I  ∼  h
    2
    		 (f0+fn) + h ×(f1+f2+f3+…+ fn−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:
    simpson.gif
    Approximate the curve that passes (−h, f(−h)), (0, f(0)) and (h, f(h)) by
    y = a x2+ b x + c.
    		(4)

    f(−h)
    =
    a h2 − b h + c,
    		(5)
    f(0)
    =
    c,
    		(6)
    f(h)
    =
    a h2+ b h + c,
    		(7)
    from which

    a
    =
    f(h)+ f(−h) − 2 f(0)
    2h2
    		 ,
    		(8)
    b
    =
    f(h)−f(−h)
    2h
    		 ,
    		(9)
    c
    =
    f(0).
    		(10)
    Hence,

    ⌠
    ⌡     		h

    −h 
    		 (a x2 + b x+ c) dx
    =
    ⌠
    ⌡     		h

    −h 
    		 ⎛
    ⎝     		f(h)+ f(−h) − 2 f(0)
    2h2
    		 x2 + f(h)−f(−h)
    2h
    		 x+ f(0)⎞
    ⎠     		dx
    =
    …
    =
    h
    3
    		 ⎛
    ⎝     		f(−h)+ 4 f(0)+ f(h)⎞
    ⎠     		.
    		(11)
    This is like a weighted average of f(−h), f(0) and f(h) as

    S ≈ f(−h)+ 4 f(0)+ f(h)
    6
    		 ×(2h).
    		(12)
    By adding this for the interval [a,b], one obtains

    I
     ∼ 
    h
    3
    		 ⎛
    ⎝     		(f0 + 4f1 + f2 )+(f2 + 4f3 + f4 )+(f4 + 4f5 + f6 )+…+(f2n−2+ 4 f2n−1+f2n )⎞
    ⎠
     ∼ 
    h
    3
    		 ⎛
    ⎝     		f0 +4 f1+2 f2 + 4f3 + 2f4 + …+ 2 f2n−2 + 4 f2n−1 + f2n⎞
    ⎠     		,
    		(13)
where
h = b−a
2n
 .
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.(13) as

I
=
h
3
 ⎛
⎝ 		f0+f2n⎞
⎠
+
h
3
 ×4 ⎛
⎝ 		f1 + f3 + f5 + …+f2n−1⎞
⎠
+
h
3
 ×2 ⎛
⎝ 		f2+f4+ f6 + …+ f2n−2⎞
⎠ 		.
(14)
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"(ξ) h2
  • Simpson's rule
    I  ∼ A + f"′(ξ) h3
where h is the step size and ξ is a value within the interval.

Computation of π


⌠
⌡ 		1

0 
 1
1+x2
 dx = π
4
 .
1overx2plus1.gif

⌠
⌡ 		1

0 
√
 
1−x2
 
 dx = π
4
 .
sqrt1minusx2.gif
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 TEX by TTH, version 4.03.
On 23 Mar 2025, 10:02.