# Chapter 21 Newtin-cotes integration formula¶

## Ex : 21.1 Pg : 612¶

In [1]:
from __future__ import division
from sympy.mpmath import quad
def f(x):
y=(0.2+25*x-200*x**2+675*x**3-900*x**4+400*x**5)
return y
tval=1.640533#
a=0#
b=0.8#
fa=f(a)#
fb=f(b)#
l=(b-a)*((fa+fb)/2)#
Et=tval-l##error
et=Et*100/tval##percent relative error

#by using approximate error estimate

#the second derivative of f
def g(x):
y=-400+4050*x-10800*x**2+8000*x**3
return y

f2x=quad(g,[0,0.8])/(b-a)##average value of second derivative
Ea=-(1/12)*(f2x)*(b-a)**3#
print "The Error Et=",round(Et,3)
print "The percent relative error et=",round(et,3),"%"
print "The approximate error estimate without using the true value=",Ea

The Error Et= 1.468
The percent relative error et= 89.467 %
The approximate error estimate without using the true value= 2.56


## Ex : 21.2 Pg : 613¶

In [2]:
from __future__ import division
from sympy.mpmath import quad

def f(x):
y=(0.2+25*x-200*x**2+675*x**3-900*x**4+400*x**5)
return y
a=0#
b=0.8#
tval=1.640533#
n=2#
h=(b-a)/n#
fa=f(a)#
fb=f(b)#
fh=f(h)#
l=(b-a)*(fa+2*fh+fb)/(2*n)#
Et=tval-l##error
et=Et*100/tval##percent relative error

#by using approximate error estimate

#the second derivative of f
def g(x):
y=-400+4050*x-10800*x**2+8000*x**3
return y
f2x=quad(g,[0,0.8])/(b-a)##average value of second derivative
Ea=-(1/12)*(f2x)*(b-a)**3/(n**2)#
print "The Error Et=",round(Et,3)
print "The percent relative error et=",round(et,3),"%"
print "The approximate error estimate without using the true value=",Ea

The Error Et= 0.572
The percent relative error et= 34.85 %
The approximate error estimate without using the true value= 0.64


## Ex :21.3 Pg : 614¶

In [3]:
from numpy import arange, exp
g=9.8##m/s**2# acceleration due to gravity
m=68.1##kg
c=12.5##kg/sec# drag coefficient
def f(t):
from numpy import exp
v=g*m*(1-exp(-c*t/m))/c
return v
tval=289.43515##m
a=0#
b=10#
fa=f(a)#
fb=f(b)#

for i in arange(10,21,10):
n=i#
h=(b-a)/n#
print "No. of segments=",i
print "Segment size=",h
j=a+h#
s=0#
while j<b:
s=s+f(j)#
j=j+h#

l=(b-a)*(fa+2*s+fb)/(2*n)#
Et=tval-l##error
et=Et*100/tval##percent relative error
print "Estimated d=",l,"m"
print et,"et(%)"
print "---------------------------------------------------------"

for i in arange(50,101,50):
n=i#
h=(b-a)/n#
print "No. of segments=",i
print "Segment size=",h
j=a+h#
s=0#
while j<b:
s=s+f(j)#
j=j+h#

l=(b-a)*(fa+2*s+fb)/(2*n)#
Et=tval-l##error
et=Et*100/tval##percent relative error
print "Estimated d=",l,"m"
print "et(%)",et
print "---------------------------------------------------------"

for i in arange(100,201,100):
n=i#
h=(b-a)/n#
print "No. of segments=",i
print "Segment size=",h
j=a+h#
s=0#
while j<b:
s=s+f(j)#
j=j+h#

l=(b-a)*(fa+2*s+fb)/(2*n)#
Et=tval-l##error
et=Et*100/tval##percent relative error
print "Estimated d=",l,"m"
print "et(%)",et
print "---------------------------------------------------------"

for i in arange(200,501,300):
n=i#
h=(b-a)/n#
print "No. of segments=",i
print "Segment size=",h
j=a+h#
s=0#
while j<b:
s=s+f(j)#
j=j+h#

l=(b-a)*(fa+2*s+fb)/(2*n)#
Et=tval-l##error
et=Et*100/tval##percent relative error
print "Estimated d=",l,"m"
print "et(%)",et
print "---------------------------------------------------------"

for i in arange(1000,2001,1000):
n=i#
h=(b-a)/n#
print "No. of segments=",i
print "Segment size=",h
j=a+h#
s=0#
while j<b:
s=s+f(j)#
j=j+h#

l=(b-a)*(fa+2*s+fb)/(2*n)#
Et=tval-l##error
et=Et*100/tval##percent relative error
print "Estimated d=",l,"m"
print "et(%)",et
print "---------------------------------------------------------"

for i in arange(2000,5001,3000):
n=i#
h=(b-a)/n#
print "No. of segments=",i
print "Segment size=",h
j=a+h#
s=0#
while j<b:
s=s+f(j)#
j=j+h#

l=(b-a)*(fa+2*s+fb)/(2*n)#
Et=tval-l##error
et=Et*100/tval##percent relative error
print "Estimated d=",l,"m"
print "et(%)",et
print "---------------------------------------------------------"

for i in arange(5000,10001,5000):
n=i#
h=(b-a)/n#
print "No. of segments=",i
print "Segment size=",h
j=a+h#
s=0#
while j<b:
s=s+f(j)#
j=j+h#

l=(b-a)*(fa+2*s+fb)/(2*n)#
Et=tval-l##error
et=Et*100/tval##percent relative error
print "Estimated d=",l,"m"
print "et(%)",et
print "---------------------------------------------------------"

No. of segments= 10
Segment size= 1.0
Estimated d= 288.749146143 m
0.237014701487 et(%)
---------------------------------------------------------
No. of segments= 20
Segment size= 0.5
Estimated d= 289.263574224 m
0.0592795228803 et(%)
---------------------------------------------------------
No. of segments= 50
Segment size= 0.2
Estimated d= 298.382319223 m
et(%) -3.09125177877
---------------------------------------------------------
No. of segments= 100
Segment size= 0.1
Estimated d= 293.915596452 m
et(%) -1.54799665905
---------------------------------------------------------
No. of segments= 100
Segment size= 0.1
Estimated d= 293.915596452 m
et(%) -1.54799665905
---------------------------------------------------------
No. of segments= 200
Segment size= 0.05
Estimated d= 289.43343055 m
et(%) 0.000594070904571
---------------------------------------------------------
No. of segments= 200
Segment size= 0.05
Estimated d= 289.43343055 m
et(%) 0.000594070904571
---------------------------------------------------------
No. of segments= 500
Segment size= 0.02
Estimated d= 290.332334709 m
et(%) -0.309977799375
---------------------------------------------------------
No. of segments= 1000
Segment size= 0.01
Estimated d= 289.883809248 m
et(%) -0.155012011658
---------------------------------------------------------
No. of segments= 2000
Segment size= 0.005
Estimated d= 289.435129352 m
et(%) 7.13401428866e-06
---------------------------------------------------------
No. of segments= 2000
Segment size= 0.005
Estimated d= 289.435129352 m
et(%) 7.13401428866e-06
---------------------------------------------------------
No. of segments= 5000
Segment size= 0.002
Estimated d= 289.435143766 m
et(%) 2.15393877364e-06
---------------------------------------------------------
No. of segments= 5000
Segment size= 0.002
Estimated d= 289.435143766 m
et(%) 2.15393877364e-06
---------------------------------------------------------
No. of segments= 10000
Segment size= 0.001
Estimated d= 289.480018962 m
et(%) -0.0155022506708
---------------------------------------------------------


## Ex : 21.4 Pg : 618¶

In [4]:
from __future__ import division
from sympy.mpmath import quad
def f(x):
y=(0.2+25*x-200*x**2+675*x**3-900*x**4+400*x**5)
return y
a=0#
b=0.8#
tval=1.640533#
n=2#
h=(b-a)/n#
fa=f(a)#
fb=f(b)#
fh=f(h)#
l=(b-a)*(fa+4*fh+fb)/(3*n)#
print"l=", round(l,2)
Et=tval-l##error
et=Et*100/tval##percent relative error

#by using approximate error estimate

#the fourth derivative of f
def g(x):
y=-21600+48000*x
return y

f4x=quad(g,[0,0.8])/(b-a)##average value of fourth derivative
Ea=-(1/2880)*(f4x)*(b-a)**5#
print "The Error Et=",round(Et,2)
print "The percent relative error et=",round(et,3),"%"
print "The approximate error estimate without using the true value=",round(Ea,3)

l= 1.37
The Error Et= 0.27
The percent relative error et= 16.645 %
The approximate error estimate without using the true value= 0.273


## Ex : 23.5 Pg : 620¶

In [16]:
from __future__ import division
from sympy.mpmath import quad
def f(x):
y=(0.2+25*x-200*x**2+675*x**3-900*x**4+400*x**5)
return y
a=0#
b=0.8#
tval=1.640533#
n=4#
h=(b-a)/n#
fa=f(a)#
fb=f(b)#
j=a+h#
s=0#
count=1#
while j<b:
if (-1)**count==-1:
s=s+4*f(j)#
else:
s=s+2*f(j)#

count=count+1#
j=j+h#

l=(b-a)*(fa+s+fb)/(3*n)#
print"l=", round(l,2)
Et=tval-l##error
et=Et*100/tval##percent relative error

#by using approximate error estimate

#the fou:rth derivative of f
def g(x):
y=-21600+48000*x
return y
f4x=quad(g,[0,0.8])/(b-a)##average value of fourth derivative
Ea=-(1/(180*4**4))*(f4x)*(b-a)**5#
print "The Error Et=",round(Et,2)
print "The percent relative error et=",round(et,3),"%"
print "The approximate error estimate without using the true value=",round(Ea,3)

l= 1.62
The Error Et= 0.02
The percent relative error et= 1.04 %
The approximate error estimate without using the true value= 0.017


## Ex :23.6 Pg : 625¶

In [18]:
from __future__ import division
from sympy.mpmath import quad
def f(x):
y=(0.2+25*x-200*x**2+675*x**3-900*x**4+400*x**5)
return y
a=0#
b=0.8#
tval=1.640533#
#part a
n=3#
h=(b-a)/n#
fa=f(a)#
fb=f(b)#
j=a+h#
s=0#
count=1#
while j<b:
s=s+3*f(j)#
count=count+1#
j=j+h#

l=(b-a)*(fa+s+fb)/(8)#
print "Part A:"
print "l=",round(l,3)
Et=tval-l##error
et=Et*100/tval##percent relative error

#by using approximate error estimate

#the fourth derivative of f
def g(x):
y=-21600+48000*x
return y
f4x=quad(g,[0,0.8])/(b-a)##average value of fourth derivative
Ea=-(1/6480)*(f4x)*(b-a)**5#
print "The Error Et=",round(Et,2)
print "The percent relative error et=",round(et,3),"%"
print "The approximate error estimate without using the true value=",round(Ea,3)
#part b
n=5#
h=(b-a)/n#
l1=(a+2*h-a)*(fa+4*f(a+h)+f(a+2*h))/6#
l2=(a+5*h-a-2*h)*(f(a+2*h)+3*(f(a+3*h)+f(a+4*h))+fb)/8#
l=l1+l2#
print "---------------------------------------------------"
print "Part B:"
print "l=", round(l,3)
Et=tval-l##error
et=Et*100/tval##percent relative error
print "The Error Et=", round(Et,3)
print "The percent relative error et=", round(et,3), "%"

Part A:
l= 1.519
The Error Et= 0.12
The percent relative error et= 7.398 %
The approximate error estimate without using the true value= 0.121
---------------------------------------------------
Part B:
l= 1.645
The Error Et= -0.005
The percent relative error et= -0.277 %


## Ex : 23.7 Pg : 626¶

In [1]:
from __future__ import division
def f(x):
y=(0.2+25*x-200*x**2+675*x**3-900*x**4+400*x**5)
return y
tval=1.640533#
x=[0, 0.12, 0.22, 0.32, 0.36, 0.4 ,0.44 ,0.54 ,0.64 ,0.7 ,0.8]
func=[]
for i in range(0,11):
func.append(f(x[i]))#

l=0#
for i in range(0,10):
l=l+(x[i+1]-x[i])*(func[i]+func[i+1])/2#

print "l=",l
Et=tval-l##error
et=Et*100/tval##percent relative error
print "The Error Et=",Et
print "The percent relative error et=",et,"%"

l= 1.59480096
The Error Et= 0.04573204
The percent relative error et= 2.78763304365 %


## Ex : 23.8 Pg : 230¶

In [2]:
def f(x):
y=(0.2+25*x-200*x**2+675*x**3-900*x**4+400*x**5)
return y
tval=1.640533#
x=[0, 0.12, 0.22, 0.32, 0.36, 0.4 ,0.44 ,0.54, 0.64, 0.7, 0.8]
func =[]
for i in range(0,11):
func.append(f(x[i]))

l1=(x[1]-x[0])*((f(x[0])+f(x[1]))/2)#
l2=(x[3]-x[1])*(f(x[3])+4*f(x[2])+f(x[1]))/6#
l3=(x[6]-x[3])*(f(x[3])+3*(f(x[4])+f(x[5]))+f(x[6]))/8#
l4=(x[8]-x[6])*(f(x[6])+4*f(x[7])+f(x[8]))/6
l5=(x[9]-x[8])*((f(x[9])+f(x[8]))/2)#
l6=(x[10]-x[9])*((f(x[10])+f(x[9]))/2)#
l=l1+l2+l3+l4+l5+l6#
print "l=",l
Et=tval-l##error
et=Et*100/tval##percent relative error
print "The Error Et=",Et
print "The percent relative error et=",et,"%"

l= 1.60364091733
The Error Et= 0.0368920826667
The percent relative error et= 2.2487863802 %


## Ex : 23.9 Pg : 629¶

In [3]:
def f(x,y):
t=2*x*y+2*x-x**2-2*y**2+72
return t
Len=8##m,length
wid=6##m,width
a=0#
b=Len#
n=2#
h=(b-a)/n#
a1=0#
b1=wid#
h1=(b1-a1)/n#

fa=f(a,0)#
fb=f(b,0)#
fh=f(h,0)#
lx1=(b-a)*(fa+2*fh+fb)/(2*n)#

fa=f(a,h1)#
fb=f(b,h1)#
fh=f(h,h1)#
lx2=(b-a)*(fa+2*fh+fb)/(2*n)#

fa=f(a,b1)#
fb=f(b,b1)#
fh=f(h,b1)#
lx3=(b-a)*(fa+2*fh+fb)/(2*n)#

l=(b1-a1)*(lx1+2*lx2+lx3)/(2*n)#

avg_temp=l/(Len*wid)#
print"The average termperature is=", avg_temp

The average termperature is= 53.0