Numerical Integration in Python
Rhythm Adhikari
080BCT065
1. Trapezoidal Rule
Code
1 import numpy as np
2 import matplotlib . pyplot as plt
3
4 def f ( x ) :
5 return np . sin ( x )
6
7 def trapezoidal_rule (a , b , n ) :
8 h = (b - a) / n
9 x = np . linspace (a , b , n + 1)
10 y = f(x)
11 result = h * (0.5 * y [0] + np . sum ( y [1: -1]) + 0.5 * y
[ -1])
12 return result , x , y
13
14 a = 0
15 b = np . pi
16 n = 100
17
18 approximation , x , y = trapezoidal_rule (a , b , n )
19 print ( " Trapezoidal rule approximation is : " , approximation )
Output
Trapezoidal rule approximation is: 1.9998355038874436
2. Simpson’s 1/3 and 3/8 Rule
Code
1 def simpsons_one_third (a , b , n ) :
2 if n % 2 != 0:
1
3 raise ValueError ( " n must be even for Simpson ’s 1/3
Rule . " )
4 h = (b - a) / n
5 result = f ( a ) + f ( b )
6 for i in range (1 , n ) :
7 x = a + i * h
8 weight = 4 if i % 2 != 0 else 2
9 result += weight * f ( x )
10 return result * ( h / 3)
11
12 def simps ons_ three _eig hth (a , b , n ) :
13 if n % 3 != 0:
14 raise ValueError ( " n must be a multiple of 3 for
Simpson ’s 3/8 Rule . " )
15 h = (b - a) / n
16 result = f ( a ) + f ( b )
17 for i in range (1 , n ) :
18 x = a + i * h
19 weight = 2 if i % 3 == 0 else 3
20 result += weight * f ( x )
21 return result * (3 * h / 8)
22
23 n1 = 100
24 n2 = 99
25 approx_1_3 = simpsons_one_third (a , b , n1 )
26 approx_3_8 = simp sons _thre e_ei ghth (a , b , n2 )
27
28 print ( f " Simpson ’s 1/3 Rule Approximation ( n ={ n1 }) : {
approx_1_3 } " )
29 print ( f " Simpson ’s 3/8 Rule Approximation ( n ={ n2 }) : {
approx_3_8 } " )
Output
Simpson’s 1/3 Rule Approximation (n=100): 2.0000000108245044
Simpson’s 3/8 Rule Approximation (n=99): 2.0000000253572914
3. Boole’s Rule
Code
1 def boole_rule (f , a , b , n ) :
2 if n < 4:
3 raise ValueError ( " n must be at least 4 " )
4 h = (b - a) / n
5 n_boole = n - ( n % 4)
6 x = np . linspace (a , a + n_boole * h , n_boole + 1)
2
7 integral = 0.0
8 for i in range (0 , n_boole , 4) :
9 x0 , x1 , x2 , x3 , x4 = x [ i : i +5]
10 fx0 , fx1 , fx2 , fx3 , fx4 = f ( x0 ) , f ( x1 ) , f ( x2 ) , f ( x3 )
, f ( x4 )
11 segment_integral = (2* h /45) *(7* fx0 + 32* fx1 + 12* fx2
+ 32* fx3 + 7* fx4 )
12 integral += segment_integral
13 return integral
14
15 result = boole_rule (f , a , b , n )
16 print ( f " Approximate integral of sin ( x ) from { a } to { b } with
n ={ n } intervals : { result } " )
Output
Approximate integral of sin(x) from 0 to 3.141592653589793 with n=100 intervals: 1.999999999
4. Weddle’s Rule
Code
1 def weddle_rule (f , a , b , n ) :
2 if n < 6:
3 raise ValueError ( " n must be at least 6 " )
4 h = (b - a) / n
5 n_weddle = n - ( n % 6)
6 x = np . linspace (a , a + n_weddle * h , n_weddle + 1)
7 integral = 0.0
8 for i in range (0 , n_weddle , 6) :
9 pts = x [ i : i +7]
10 fpts = f ( pts )
11 segment_integral = (3* h /10) *( fpts [0] + 5* fpts [1] +
fpts [2] +
12 6* fpts [3] + fpts [4] + 5* fpts [5]
+ fpts [6])
13 integral += segment_integral
14 return integral
15
16 n = 60
17 result = weddle_rule (f , a , b , n )
18 print ( f " Approximate integral of sin ( x ) from { a } to { b } with
n ={ n } intervals using Weddle ’s Rule : { result } " )
Output
Approximate integral of sin(x) from 0 to 3.141592653589793 with n=60 intervals using Weddle’
3
5. Gauss-Legendre 3-Point Rule
Code
1 def gauss_legendre_3pt (f , a , b ) :
2 t = np . array ([ - np . sqrt (3/5) , 0 , np . sqrt (3/5) ])
3 w = np . array ([5/9 , 8/9 , 5/9])
4 x = (( b - a ) / 2) * t + ( a + b ) / 2
5 integral = (( b - a ) / 2) * np . sum ( w * f ( x ) )
6 return integral
7
8 result = gauss_legendre_3pt (f , a , b )
9 print ( f " Approximate integral of sin ( x ) from { a } to { b } using
3 - point Gauss - Legendre : { result } " )
Output
Approximate integral of sin(x) from 0 to 3.141592653589793 using 3-point Gauss-Legendre: 2.0
