Chapter 12 - Numerical integration

Example No. 12_01 Pg No. 373

In [1]:
from scipy.misc import derivative
from mpmath import quad
#Trapezoidal Rule

def f(x):
    F = x**3 + 1
    return F

#case(a)
a = 1#
b = 2 #
h = b - a #
It = (b-a)*(f(a)+f(b))/2
d2f = derivative(f,2,n=2)
Ett = h**3*d2f/12.0
Iexact = quad(f,[1,2])
Trueerror = It - Iexact
print 'case(a)'
print 'It = ',It
print 'Ett = ',Ett
print 'Iexact = ',Iexact
print 'True error = ',Trueerror,
print 'Here Error bound is an overestimate of true error'

#case(b)
a = 1#
b = 1.5 #
h = b - a #
It = (b-a)*(f(a)+f(b))/2
Ett = h**3*derivative(f,1.5,n=2)/12
Iexact = quad(f,[1,1.5])
Trueerror = It - Iexact
print 'case(b)'
print 'It = ',It
print 'Ett = ',Ett
print 'Iexact = ',Iexact
print 'True error = ',Trueerror,
print 'Here Error bound is an overestimate of true error'

case(a)
It =  5
Ett =  1.0
Iexact =  4.75
True error =  0.25 Here Error bound is an overestimate of true error
case(b)
It =  1.59375
Ett =  0.09375
Iexact =  1.515625
True error =  0.078125 Here Error bound is an overestimate of true error

Example No. 12_02 Pg No. 376

In [2]:
from math import exp
#Tapezoidal rule


def f(x):
    F = exp(x)
    return F
a = -1 #
b = 1 #

#case(a)
n = 2
h = (b-a)/n 
I = 0
for i in range(1,n+1):
    I = I + f(a+(i-1)*h)+f(a+i*h)

I = h*I/2 #
print 'intergral for case(a),Ia = ',I

#case(b)
n = 4
h = (b-a)/n 
I = 0
for i in range(1,n+1):
    I = I + f(a+(i-1)*h)+f(a+i*h)#

I = h*I/2 #
Iexact = 2.35040
print 'intergral for case(b),Ib = ',I
print 'exact integral,Iexact = ',Iexact
print 'n = 4 case is better than n = 2 case'

intergral for case(a),Ia =  2.54308063482
intergral for case(b),Ib =  0.0
exact integral,Iexact =  2.3504
n = 4 case is better than n = 2 case

Example No. 12_03 Pg No. 381

In [3]:
from math import pi,sin,exp,sqrt
#Simpon's 1/3 rule

#case(a)
def f(x):
    F = exp(x)
    return F
a = -1#
b = 1#
h = (b-a)/2 
x1 = a+h
Is1 = h*( f(a) + f(b) + 4*f(x1) )/3 
print 'h = ',h
print 'Integral for case(a) , Is1 = ',Is1

#case(b)
def f(x):
    F = sqrt(sin(x))
    return F
a = 0
b = pi/2
h = (b-a)/2 
x1 = a+h
Is1 = h*( f(a) + f(b) + 4*f(x1) )/3
print 'h = ',h
print 'Integral for case(b),Is1 = ',Is1

h =  1
Integral for case(a) , Is1 =  2.36205375654
h =  0.785398163397
Integral for case(b),Is1 =  1.14238405466

Example No. 12_04 Pg No.382

In [4]:
from math import pi,sin,exp,sqrt
#Simpon's 1/3 rule

def f(x):
    F = sqrt( sin(x) )
    return F
x0 = 0 #
xa = pi/2 #

#case(a) n = 4
n = 4 #
h = ( xa-x0 )/n
I = 0 
for i in range(1,n/2+1):
    I = I + f(x0 + (2*i-2)*h) + 4*f(x0 + (2*i-1)*h) + f(x0 + 2*i*h) #
I = h*I/3
print 'h = ',h
print 'Integral value for n = 4 is',I 

#case(b) n = 6
n = 6
h = ( xa-x0 )/n
I = 0 
for i in range(1,n/2+1):
    I = I + f(x0 + (2*i-2)*h) + 4*f(x0 + (2*i-1)*h) + f(x0 + 2*i*h) #

I = h*I/3
print 'h = ',h
print 'Integral value for n = 6 is',I

h =  0.392699081699
Integral value for n = 4 is 1.17822754435
h =  0.261799387799
Integral value for n = 6 is 1.18728120346

Example No. 12_05 Pg No. 386

In [5]:
from math import pi,sin,exp,sqrt

#Simpson's 3/8 rule

#case(a)
def f(x):
    F = x**3 + 1
    return F
a = 1 #
b = 2 #
h = (b-a)/3 
x1 = a + h 
x2 = a + 2*h
Is2 = 3*h*( f(a) + 3*f(x1) + 3*f(x2) + f(b) )/8 #
print 'Integral of x**3 +1 from 1 to 2 : ',Is2
#case(b)
def f(x):
    F = sqrt( sin(x) )
    return F
a = 0 #
b = pi/2 #
h = (b-a)/3 
x1 = a + h 
x2 = a + 2*h
Is2 = 3*h*( f(a) + 3*f(x1) + 3*f(x2) + f(b) )/8 #
print 'Integral of sqrt(sin(x)) from 0 to pi/2:',Is2

Integral of x**3 +1 from 1 to 2 :  0
Integral of sqrt(sin(x)) from 0 to pi/2: 1.16104132669

Example No. 12_06 Pg No. 387

In [6]:
from math import sin,sqrt,pi
#Booles's Five-Point formula

def f(x):
    F = sqrt( sin(x) )
    return F
x0 = 0#
xb = pi/2 #
n = 4 #
h = (xb - x0)/n
Ib = 2*h*(7*f(x0) + 32*f(x0+h) + 12*f(x0 + 2*h) + 32*f(x0+3*h) + 7*f(x0+4*h))/45#
print 'Ib = ',Ib

Ib =  1.18061711033

Example No. 12_07 Pg No. 391

In [7]:
from __future__ import division
from numpy import zeros
def f(x):
    F = 1.0/x
    return F
#since we can't have (0,0) element in matrix we start with (1,1)
a = 1 ;
b = 2 ;
h = b-a ;
R=zeros([3,3])
R[0,0] = h*(f(a)+f(b))/2 
print 'R(0,0) = ',R[0,0]

h=[h]
for i in range(2,4):
    h.append((b-a)/2**(i-1))
    s = 0
    for k in range(1,2**(i-2)+1):
        s = s + f(a + (2*k - 1)*h[i-1]);
    
    R[i-1,0] = R[i-1,0]/2 + h[i-1]*s;
    print '\nR(%i,0) = %f \n'%(i-1,R[i-1,0])

for j in range(2,4):
    for i in range(j,4):
        R[i-1,j-1] = (4**(j-1)*R[i-1,j-2] - R[i-2,j-2])/(4**(j-1)-1);
        print '\nR(%i,%i) = %f \n'%(i-1,j-1,R[i-1,j-1])

R(0,0) =  0.75

R(1,0) = 0.333333 


R(2,0) = 0.342857 


R(1,1) = 0.194444 


R(2,1) = 0.346032 


R(2,2) = 0.356138

Example No. 12_08 Pg No. 397

In [8]:
from math import sqrt,exp 
#Two Point Gauss -Legedre formula

def f(x):
    F = exp(x)
    return F
x1 = -1/sqrt(3)
x2 = 1/sqrt(3)
I = f(x1) + f(x2)
print 'x1 = ',x1
print 'x2 = ',x2
print 'I = ',I, 

x1 =  -0.57735026919
x2 =  0.57735026919
I =  2.34269608791

Example No. 12_09 Pg No. 398

In [9]:
from math import sqrt,exp
#Gaussian two point formula

a = -2 #
b = 2 #
def f(x):
    F = exp(-x/2)
    return F
A = (b-a)/2 
B = (a+b)/2
C = (b-a)/2
def g(z):
    G = exp(-1*(A*z+B)/2)
    return G
w1 = 1 #
w2 = 1 #
z1 = -1/sqrt(3)
z2 = 1/sqrt(3)
Ig = C*( w1*g(z1) + w2*g(z2) )
print ' g(z) = exp(-(%f*z + %f)/2) \n C = %f \n Ig = %f \n'%(A,B,C,Ig)

 g(z) = exp(-(2.000000*z + 0.000000)/2) 
 C = 2.000000 
 Ig = 4.685392