Chapter06:Numerical Differentiation and Integration

Ex6.1:pg-201

In [1]:
#example 6.1
#numerical diffrentiation by newton's difference formula 
#page 210
x=[1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2]
y=[2.7183, 3.3201, 4.0552, 4.9530, 6.0496, 7.3891, 9.0250]
c=0
d1=[0,0,0,0,0,0]
d2=[0,0,0,0,0]
d3=[0,0,0,0]
d4=[0,0,0]
d5=[0,0]
d6=[0]
for i in range(0,6):
    d1[c]=y[i+1]-y[i]
    c=c+1;
c=0
for i in range(0,5):
    d2[c]=d1[i+1]-d1[i]
    c=c+1;
c=0
for i in range(0,4):
    d3[c]=d2[i+1]-d2[i]
    c=c+1;
c=0
for i in range(0,3):
    d4[c]=d3[i+1]-d3[i]
    c=c+1;
c=0
for i in range(0,2):
    d5[c]=d4[i+1]-d4[i]
    c=c+1;
c=0
for i in range(0,1):
    d6[c]=d5[i+1]-d5[i]
    c=c+1;
x0=1.2  #first and second derivative at 1.2
h=0.2
f1=((d1[1]-d2[1]/2+d3[1]/3-d4[1]/4+d5[1]/5)/h)
print "the first derivative of fuction at 1.2 is:%f\n" %(f1)
f2=(d2[1]-d3[1]+(11*d4[1])/12-(5*d5[1])/6)/h**2
print "the second derivative of fuction at 1.2 is:%f\n" %(f2)

the first derivative of fuction at 1.2 is:3.320317

the second derivative of fuction at 1.2 is:3.319167

Ex6.2:pg-211

In [6]:
#example 6.2
#numerical diffrentiation by newton's difference formula 
#page 211
x=[1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2]
y=[2.7183, 3.3201, 4.0552, 4.9530, 6.0496, 7.3891, 9.0250]
c=0
d1=[0,0,0,0,0,0]
d2=[0,0,0,0,0]
d3=[0,0,0,0]
d4=[0,0,0]
d5=[0,0]
d6=[0]
for i in range(0,6):
    d1[c]=y[i+1]-y[i]
    c=c+1;
c=0
for i in range(0,5):
    d2[c]=d1[i+1]-d1[i]
    c=c+1;
c=0
for i in range(0,4):
    d3[c]=d2[i+1]-d2[i]
    c=c+1;
c=0
for i in range(0,3):
    d4[c]=d3[i+1]-d3[i]
    c=c+1;
c=0
for i in range(0,2):
    d5[c]=d4[i+1]-d4[i]
    c=c+1;
c=0
for i in range(0,1):
    d6[c]=d5[i+1]-d5[i]
    c=c+1;
x0=2.2  #first and second derivative at 2.2
h=0.2
f1=((d1[5]+d2[4]/2+d3[3]/3+d4[2]/4+d5[1]/5)/h)
print "the first derivative of fuction at 1.2 is:%f\n" %(f1)
f2=(d2[4]+d3[3]+(11*d4[2])/12+(5*d5[1])/6)/h**2
print "the second derivative of fuction at 1.2 is:%f\n" %(f2)
x1=2.0  # first derivative also at 2.0
f1=((d1[4]+d2[3]/2+d3[2]/3+d4[1]/4+d5[0]/5+d6[0]/6)/h)
print "the first derivative of function at 1.2 is:%f\n" %(f1)

the first derivative of fuction at 1.2 is:9.022817

the second derivative of fuction at 1.2 is:8.992083

the first derivative of function at 1.2 is:7.389633

Ex6.3:pg-211

In [8]:
#example 6.3
#numerical diffrentiation by newton's difference formula 
#page 211
x=[1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2]
y=[2.7183, 3.3201, 4.0552, 4.9530, 6.0496, 7.3891, 9.0250]
c=0
d1=[0,0,0,0,0,0]
d2=[0,0,0,0,0]
d3=[0,0,0,0]
d4=[0,0,0]
d5=[0,0]
d6=[0]
for i in range(0,6):
    d1[c]=y[i+1]-y[i]
    c=c+1;
c=0
for i in range(0,5):
    d2[c]=d1[i+1]-d1[i]
    c=c+1;
c=0
for i in range(0,4):
    d3[c]=d2[i+1]-d2[i]
    c=c+1;
c=0
for i in range(0,3):
    d4[c]=d3[i+1]-d3[i]
    c=c+1;
c=0
for i in range(0,2):
    d5[c]=d4[i+1]-d4[i]
    c=c+1;
c=0
for i in range(0,1):
    d6[c]=d5[i+1]-d5[i]
    c=c+1;
x0=1.6  #first and second derivative at 1.6
h=0.2
f1=(((d1[2]+d1[3])/2-(d3[1]+d3[2])/4+(d5[0]+d5[1])/60))/h
print "the first derivative of function at 1.6 is:%f\n" %(f1)
f2=((d2[2]-d4[1]/12)+d6[0]/90)/(h**2)
print "the second derivative of function at 1.6 is:%f\n" %(f2)

the first derivative of function at 1.6 is:4.885975

the second derivative of function at 1.6 is:4.953361

Ex6.4:pg-213

In [9]:
#example 6.4
#estimation of errors 
#page 213
x=[1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2]
y=[2.7183, 3.3201, 4.0552, 4.9530, 6.0496, 7.3891, 9.0250]
c=0
d1=[0,0,0,0,0,0]
d2=[0,0,0,0,0]
d3=[0,0,0,0]
d4=[0,0,0]
d5=[0,0]
d6=[0]
for i in range(0,6):
    d1[c]=y[i+1]-y[i]
    c=c+1;
c=0
for i in range(0,5):
    d2[c]=d1[i+1]-d1[i]
    c=c+1;
c=0
for i in range(0,4):
    d3[c]=d2[i+1]-d2[i]
    c=c+1;
c=0
for i in range(0,3):
    d4[c]=d3[i+1]-d3[i]
    c=c+1;
c=0
for i in range(0,2):
    d5[c]=d4[i+1]-d4[i]
    c=c+1;
c=0
for i in range(0,1):
    d6[c]=d5[i+1]-d5[i]
    c=c+1
x0=1.6      #first and second derivative at 1.6
h=0.2
f1=((d1[1]-d2[1]/2+d3[1]/3-d4[1]/4+d5[1]/5)/h)
print "the first derivative of fuction at 1.2 is:%f\n" %(f1)
f2=(d2[1]-d3[1]+(11*d4[1])/12-(5*d5[1])/6)/h**2
print "the second derivative of fuction at 1.2 is:%f\n" %(f2)
T_error1=((d3[1]+d3[2])/2)/(6*h)   #truncation error
e=0.00005   #corrected to 4D values
R_error1=(3*e)/(2*h)
T_error1=T_error1+R_error1    #total error
f11=(d1[2]+d1[3])/(2*h) #using stirling formula first derivative
f22=d2[2]/(h*h)#second derivative
T_error2=d4[1]/(12*h*h)
R_error2=(4*e)/(h*h)
T_error2=T_error2+R_error2
print "total error in first derivative is %0.4g:\n" %(T_error1)
print "total error in second derivative is %0.4g:" %(T_error2)

the first derivative of fuction at 1.2 is:3.320317

the second derivative of fuction at 1.2 is:3.319167

total error in first derivative is 0.03379:

total error in second derivative is 0.02167:

Ex6.5:pg-214

In [15]:
#cubic spline method
#example 6.5
#page 214
import math
from __future__ import division
x=[0, math.pi/2, math.pi]
y=[0, 1, 0]
M0=0
M2=0
h=math.pi/2
M1=(6*(y[0]-2*y[1]+y[2])/(h**2)-M0-M2)/4
def s1(x):
    return (2/math.pi)*(-2*3*x*x/(math.pi**2)+3/2)
S1=s1(math.pi/4)
print "S1(pi/4)=%f" %(S1)
def s2(x):
    return (-24*x)/(math.pi**3)
S2=s2(math.pi/4)
print "S2(pi/4)=%f" %(S2)

S1(pi/4)=0.716197
S2(pi/4)=-0.607927

Ex6.6:pg-216

In [18]:
#derivative by cubic spline method
#example 6.6
#page 216
x=[-2, -1, 2, 3]
y=[-12, -8, 3, 5] 
def f(x):
    return x**3/15-3*x**2/20+241*x/60-3.9
def s2(x):
    return (((2-x)**3)/6*(14/55)+((x+1)**3)/6*(-74/55))/3+(-8-21/55)*(2-x)/3+(3-(9/6)*(-74/55))*(x+1)/3
h=0.0001
x0=1.0
y1=(s2(x0+h)-s2(x0))/h
print "the value y1(%0.2f) is : %f" %(x0,y1)

the value y1(1.00) is : 3.527232

Ex6.7:pg-218

In [26]:
#maximun and minimun of functions
#example 6.7
#page 218
x=[1.2, 1.3, 1.4, 1.5, 1.6]
y=[0.9320, 0.9636, 0.9855, 0.9975, 0.9996]
d1=[0,0,0,0]
d2=[0,0,0]
for i in range(0,4):
    d1[i]=y[i+1]-y[i]
for i in range(0,3):
    d2[i]=d1[i+1]-d1[i]
p=(-d1[0]*2/d2[0]+1)/2;
print "p=%f" %(p)
h=0.1
x0=1.2
X=x0+p*h
print " the value of X correct to 2 decimal places is : %0.2f" %(X)
Y=y[4]-0.2*d1[3]+(-0.2)*(-0.2+1)*d2[2]/2
print "the value Y=%f" %(Y)

p=3.757732
 the value of X correct to 2 decimal places is : 1.58
the value Y=0.999972

Ex6.8:pg-226

In [25]:
#example 6.8
#trapezoidal method for integration
#page 226
from __future__ import division
x=[7.47, 7.48, 7.49, 7.0, 7.51, 7.52]
f_x=[1.93, 1.95, 1.98, 2.01, 2.03, 2.06]
h=x[1]-x[0]
l=6
area=0
for i in range(0,l):
    if i==0:
       area=area+f_x[i]
    elif i==l-1:
        area=area+f_x[i]
    else:
       area=area+2*f_x[i]
area=area*(h/2)
print "area bounded by the curve is %f" %(area)

area bounded by the curve is 0.099650

Ex6.9:pg-226

In [8]:
#example 6.9
#simpson 1/3rd  method for integration
#page 226
from __future__ import division
import math
x=[0,0.00, 0.25, 0.50, 0.75, 1.00]
y=[0,1.000, 0.9896, 0.9589, 0.9089, 0.8415]
h=x[2]-x[1]
area=0
for i in range(0,6):
    y[i]=y[i]**2
for i in range(1,6):
    if i==1:
        area=area+y[i]
    elif i==5:
        area=area+y[i]
    elif i%2==0:
        area=area+4*y[i]
    elif i%2!=0: 
        area=area+2*y[i]
area=(area/3)*(h*math.pi)
print "area bounded by the curve is %f" %(area)

area bounded by the curve is 2.819247

Ex6.10:pg-228

In [36]:
#example 6.10
#integration by trapezoidal and simpson's method
#page 228
from __future__ import division
def f(x):
    return 1/(1+x)
h=0.5
x=[0,0.0,0.5,1.0]
y=[0,0,0,0]
l=4
for i in range(0,l):
    y[i]=f(x[i])
area=0   #trapezoidal method
for i in range(1,l):
    if i==1:
       area=area+y[i]
    elif i==l-1:
        area=area+y[i]
    else:
       area=area+2*y[i]
area=area*(h/2)
print "area bounded by the curve by trapezoidal method with h=%f is %f\n \n" %(h,area)
area=0  #simpson 1/3rd rule
for i in range(1,l):
    if i==1: 
       area=area+y[i]
    elif i==l-1:
        area=area+y[i]
    elif i%2==0:
        area=area+4*y[i]
    elif i%2!=0:
        area=area+2*y[i]
area=(area*h)/3
print "area bounded by the curve by simpson 1/3rd method with h=%f is %f\n \n" %(h,area)
h=0.25
x=[0,0.0,0.25,0.5,0.75,1.0]
y=[0,0,0,0,0,0]
l=6
for i in range(0,l):
    y[i]=f(x[i])
area=0   #trapezoidal method
for i in range(1,l):
    if i==1: 
       area=area+y[i]
    elif i==l-1:
        area=area+y[i]
    else:
        area=area+2*y[i]
area=area*(h/2)
print "area bounded by the curve by trapezoidal method with h=%f is %f\n \n" %(h,area)
area=0  #simpson 1/3rd rule
for i in range(1,l):
    if i==1:
       area=area+y[i]
    elif i==l-1:
        area=area+y[i]
    elif i%2==0:
        area=area+4*y[i]
    elif i%2!=0:
        area=area+2*y[i]
area=(area*h)/3
print "area bounded by the curve by simpson 1/3rd method with h=%f is %f\n \n" %(h,area)
h=0.125
x=[0,0.0,0.125,0.25,0.375,0.5,0.625,0.75,0.875,1.0]
y=[0,0,0,0,0,0,0,0,0,0]
l=10
for i in range(0,l):
    y[i]=f(x[i])
area=0   #trapezoidal method
for i in range(1,l):
    if i==1:
        area=area+y[i]
    elif i==l-1:
        area=area+y[i]
    elif i%2==0:
        area=area+2*y[i]
    elif i%2!=0:
        area=area+2*y[i]
area=area*(h/2)
print "area bounded by the curve by trapezoidal method with h=%f is %f\n \n" %(h,area)
area=0  #simpson 1/3rd rule
for i in range(1,l):
    if i==1:
       area=area+y[i]
    elif i==l-1:
        area=area+y[i]
    elif i%2==0:
        area=area+4*y[i]
    elif i%2!=0:
        area=area+2*y[i]
area=(area*h)/3
print "area bounded by the curve by simpson 1/3rd method with h=%f is %f\n \n" %(h,area)





      

area bounded by the curve by trapezoidal method with h=0.500000 is 0.708333
 

area bounded by the curve by simpson 1/3rd method with h=0.500000 is 0.694444
 

area bounded by the curve by trapezoidal method with h=0.250000 is 0.697024
 

area bounded by the curve by simpson 1/3rd method with h=0.250000 is 0.693254
 

area bounded by the curve by trapezoidal method with h=0.125000 is 0.694122
 

area bounded by the curve by simpson 1/3rd method with h=0.125000 is 0.693155

Ex6.11:pg-229

In [43]:
#example 6.11
#rommberg's method
#page 229
from __future__ import division
def f(x):
    return 1/(1+x)
k=0
h=0.5
x=[0,0.0,0.5,1.0]
y=[0,0,0,0]
I=[0,0,0]
I1=[0,0]
T2=[0]
l=4
for i in range(0,l):
    y[i]=f(x[i])
area=0   #trapezoidal method
for i in range(1,l):
    if i==1:
       area=area+y[i]
    elif i==l-1:
        area=area+y[i]
    else:
        area=area+2*y[i]
area=area*(h/2)
I[k]=area
k=k+1
h=0.25
x=[0,0.0,0.25,0.5,0.75,1.0]
y=[0,0,0,0,0,0]
l=6
for i in range(0,l):
    y[i]=f(x[i])
area=0   #trapezoidal method
for i in range(1,l):
    if i==1:
       area=area+y[i]
    elif i==l-1:
        area=area+y[i]
    else:
        area=area+2*y[i]
area=area*(h/2)
I[k]=area
k=k+1
h=0.125
x=[0,0.0,0.125,0.25,0.375,0.5,0.625,0.75,0.875,1.0]
y=[0,0,0,0,0,0,0,0,0,0]
l=10
for i in range(0,l):
    y[i]=f(x[i])
area=0   #trapezoidal method
for i in range(1,l):
    if i==1:
       area=area+y[i]
    elif i==l-1:
        area=area+y[i]
    else:
        area=area+2*y[i]
area=area*(h/2)
I[k]=area
k=k+1
print "results obtained with h=0.5 0.25 0.125 is %f %f %f\n \n" %(I[0],I[1],I[2])
for i in range(0,2):
    I1[i]=I[i+1]+(I[i+1]-I[i])/3
for i in range(0,1):
    T2[i]=I1[i+1]+(I1[i+1]-I1[i])/3
print "the area is %f" %(T2[0])

results obtained with h=0.5 0.25 0.125 is 0.708333 0.697024 0.694122
 

the area is 0.693121

Ex6.13:pg-230

In [45]:
#area using cubic spline method
#example 6.13
#page 230
x=[0, 0.5, 1.0]
y=[0, 1.0, 0.0]
h=0.5
M0=0
M2=0
M=[0,0,0]
M1=(6*(y[2]-2*y[1]+y[0])/h**2-M0-M2)/4
M=[M0, M1, M2]
I=0
for i in range(0,2):
    I=I+(h*(y[i]+y[i+1]))/2-((h**3)*(M[i]+M[i+1])/24)
print "the value of the integrand is : %f" %(I)

the value of the integrand is : 0.625000

Ex6.15:pg-233

In [47]:
#euler's maclaurin formula
#example 6.15
#page 233
import math
y=[0, 1, 0]
h=math.pi/4
I=h*(y[0]+2*y[1]+y[2])/2+(h**2)/12+(h**4)/720
print "the value of integrand with h=%f is : %f\n\n" %(h,I)
h=math.pi/8
y=[0, math.sin(math.pi/8), math.sin(math.pi*2/8), math.sin(math.pi*3/8), math.sin(math.pi*4/8)]
I=h*(y[0]+2*y[1]+2*y[2]+2*y[3]+y[4])/2+(h**2)/2+(h**2)/12+(h**4)/720
print " the value of integrand with h=%f is : %f" %(h,I)

the value of integrand with h=0.785398 is : 0.837331


 the value of integrand with h=0.392699 is : 1.077106

Ex6.17:pg-236

In [49]:
# example 6.17
# error estimate in evaluation of the integral
# page 236
import math
def f(a,b):
    return math.cos(a)+4*math.cos((a+b)/2)+math.cos(b)
a=0
b=math.pi/2
c=math.pi/4
I=[0,0,0]
I[0]=(f(a,b)*((b-a)/2)/3)
I[1]=(f(a,c)*((c-a)/2)/3)
I[2]=(f(c,b)*((b-c)/2)/3)
Area=I[1]+I[2]
Error_estimate=((I[0]-I[1]-I[2])/15)
Actual_area=math.sin(math.pi/2)-math.sin(0)
Actual_error=abs(Actual_area-Area)
print "the calculated area obtained is:%f\n" %(Area)
print "the actual area obtained is:%f\n" %(Actual_area)
print "the actual error obtained is:%f\n" %(Actual_error)

the calculated area obtained is:1.000135

the actual area obtained is:1.000000

the actual error obtained is:0.000135

Ex6.18:pg-237

In [50]:
# example 6.18
# error estimate in evaluation of the integral
# page 237
import math
def f(a,b):
    return 8+4*math.sin(a)+4*(8+4*math.sin((a+b)/2))+8+4*math.sin(b)
a=0
b=math.pi/2
c=math.pi/4
I=[0,0,0]
I[0]=(f(a,b)*((b-a)/2)/3)
I[1]=(f(a,c)*((c-a)/2)/3)
I[2]=(f(c,b)*((b-c)/2)/3)
Area=I[1]+I[2]
Error_estimate=((I[0]-I[1]-I[2])/15)
Actual_area=8*math.pi/2+4*math.sin(math.pi/2)
Actual_error=abs(Actual_area-Area)
print "the calculated area obtained is:%f\n" %(Area)
print "the actual area obtained is:%f\n" %(Actual_area)
print "the actual error obtained is:%f\n" %(Actual_error)

the calculated area obtained is:16.566909

the actual area obtained is:16.566371

the actual error obtained is:0.000538

Ex6.19:pg-242

In [51]:
#gauss' formula
#example 6.19
#page 242
u=[-0.86113, -0.33998, 0.33998, 0.86113]
W=[0.34785, 0.65214, 0.65214, 0.34785]
I=0
for i in range(0,4):
    I=I+(u[i]+1)*W[i]
I=I/4
print " the value of integrand is : %0.5f" %(I)

 the value of integrand is : 0.49999

Ex6.20:pg-247

In [55]:
#example 6.20
#double integration
#page 247
import math
def f(x,y):
    return exp(x+y)
h0=0.5
k0=0.5
x=[[0,0,0],[0,0,0],[0,0,0]]
h=[0, 0.5, 1]
k=[0, 0.5, 1]
for i in range(0,3):
    for j in range(0,3):
      x[i][j]=f(h[i],k[j])
T_area=h0*k0*(x[0][0]+4*x[0][1]+4*x[2][1]+6*x[0][2]+x[2][2])/4     #trapezoidal method
print "the integration value by trapezoidal method  is %f\n " %(T_area)
S_area=h0*k0*((x[0][0]+x[0][2]+x[2][0]+x[2][2]+4*(x[0][1]+x[2][1]+x[1][2]+x[1][0])+16*x[1][1]))/9
print "the integration value by Simpson method  is %f" %(S_area)

the integration value by trapezoidal method  is 3.076274
 
the integration value by Simpson method  is 2.954484