# Chapter06:Numerical Differentiation and Integration¶

## Ex6.1:pg-201¶

In :
#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=
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-d2/2+d3/3-d4/4+d5/5)/h)
print "the first derivative of fuction at 1.2 is:%f\n" %(f1)
f2=(d2-d3+(11*d4)/12-(5*d5)/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 :
#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=
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+d2/2+d3/3+d4/4+d5/5)/h)
print "the first derivative of fuction at 1.2 is:%f\n" %(f1)
f2=(d2+d3+(11*d4)/12+(5*d5)/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+d2/2+d3/3+d4/4+d5/5+d6/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 :
#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=
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+d1)/2-(d3+d3)/4+(d5+d5)/60))/h
print "the first derivative of function at 1.6 is:%f\n" %(f1)
f2=((d2-d4/12)+d6/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 :
#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=
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-d2/2+d3/3-d4/4+d5/5)/h)
print "the first derivative of fuction at 1.2 is:%f\n" %(f1)
f2=(d2-d3+(11*d4)/12-(5*d5)/6)/h**2
print "the second derivative of fuction at 1.2 is:%f\n" %(f2)
T_error1=((d3+d3)/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+d1)/(2*h) #using stirling formula first derivative
f22=d2/(h*h)#second derivative
T_error2=d4/(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 :
#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-2*y+y)/(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 :
#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 :
#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*2/d2+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-0.2*d1+(-0.2)*(-0.2+1)*d2/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 :
#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-x
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 :
#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-x
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 :
#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 :
#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=
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,I,I)
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)

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 :
#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*y+y)/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 :
#euler's maclaurin formula
#example 6.15
#page 233
import math
y=[0, 1, 0]
h=math.pi/4
I=h*(y+2*y+y)/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+2*y+2*y+2*y+y)/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 :
# 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=(f(a,b)*((b-a)/2)/3)
I=(f(a,c)*((c-a)/2)/3)
I=(f(c,b)*((b-c)/2)/3)
Area=I+I
Error_estimate=((I-I-I)/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 :
# 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=(f(a,b)*((b-a)/2)/3)
I=(f(a,c)*((c-a)/2)/3)
I=(f(c,b)*((b-c)/2)/3)
Area=I+I
Error_estimate=((I-I-I)/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 :
#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 :
#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+4*x+4*x+6*x+x)/4     #trapezoidal method
print "the integration value by trapezoidal method  is %f\n " %(T_area)
S_area=h0*k0*((x+x+x+x+4*(x+x+x+x)+16*x))/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