Chapter-23 : Numerical Differentiation

Ex:23.1 Pg: 656

In [1]:
from numpy import arange
#f(x) = -0.1*x**4 - 0.15*x**3 - 0.5 * x**2 - 0.25 *x + 1.2
h = 0.25#
t = -0.9125#
x = arange(0,1.1,h)
print "x = ",x
fx=[0]
for xx in x:
    fx.append(-0.1*xx**4 - 0.15*xx**3 - 0.5 * xx**2 - 0.25 *xx + 1.2)
print "f(x) = ",fx
fd = (- fx[(5)] + 4*fx[(4)] - 3 * fx[(3)])/(2 * h)
efd = (t - fd) * 100 / t#
print "by forward difference : ",fd
print "error in forward difference method = ",efd,"%"
bd = (3 * fx[(3)] - 4 * fx[(2)] + fx[(1)])/ (2*h)
ebd = (t - bd) * 100 / t#
print "by backward difference",bd
print "error in backward difference method = ",ebd,"%"
cdm = (-fx[(5)] + 8*(fx[(4)]) -8*fx[(2)] + fx[(1)] ) / (12*h)
ecdm = (t - cdm) * 100 / t
print "by central difference : ",cdm
print "error in central difference method = ",ecdm,"%"
x =  [ 0.    0.25  0.5   0.75  1.  ]
f(x) =  [0, 1.2, 1.103515625, 0.92499999999999993, 0.63632812499999991, 0.19999999999999996]
by forward difference :  -0.859375
error in forward difference method =  5.82191780822 %
by backward difference -0.878125
error in backward difference method =  3.76712328767 %
by central difference :  -0.9125
error in central difference method =  -2.43336553343e-14 %

Ex:23.2 Pg: 657

In [3]:
from numpy import arange
#f(x) = -0.1*x**4 - 0.15*x**3 - 0.5 * x**2 - 0.25 *x + 1.2
h = 0.5#
t = -0.9125#
x = arange(0,1.1,h)
print "x with h = 0.5 is",x
fx=[0]
for xx in x:
    fx.append(-0.1*xx**4 - 0.15*xx**3 - 0.5 * xx**2 - 0.25 *xx + 1.2)
print "f(x) with h = 0.5 is",fx
cdm = (fx[(3)] - fx[(1)])/ 1#
ecdm = (t - cdm) * 100 / t#
print "by central difference ( h = 0.5 ) ",cdm
print "error in central difference method ( h = 0.5 ) = ",ecdm,"%"
h1 = 0.25#
x1 = arange(0,1.1,h1)
print "x with h = 0.25 is",x1
fx1=[0]
for xx in x1:
    fx1.append(-0.1*xx**4 - 0.15*xx**3 - 0.5 * xx**2 - 0.25 *xx + 1.2)
print "fx with h = 0.25 is",fx1
cdm1 = (fx1[(4)] - fx1[(2)])/ (2*h1)
ecdm1 = (t - cdm1) * 100 / t#
print "by central difference ( h = 0.25 )  = ",cdm1,
print "error in central difference method ( h = 0.25 ) = ",ecdm1,"%"
D = 4 * cdm1 /3 - cdm / 3#
print "improved estimate =",D
x with h = 0.5 is [ 0.   0.5  1. ]
f(x) with h = 0.5 is [0, 1.2, 0.92499999999999993, 0.19999999999999996]
by central difference ( h = 0.5 )  -1.0
error in central difference method ( h = 0.5 ) =  -9.58904109589 %
x with h = 0.25 is [ 0.    0.25  0.5   0.75  1.  ]
fx with h = 0.25 is [0, 1.2, 1.103515625, 0.92499999999999993, 0.63632812499999991, 0.19999999999999996]
by central difference ( h = 0.25 )  =  -0.934375 error in central difference method ( h = 0.25 ) =  -2.39726027397 %
improved estimate = -0.9125

Ex:23.3 Pg: 658

In [4]:
#q(z = 0) = -k*p*C*(dT/dz)/(z = 0)
k = 3.5 * 10** - 7##m**2/s
p = 1800##kg/m**3
C = 840# #(J/(kg.C))
x = 0#
fx0 = 13.5#
fx1 = 12#
fx2 = 10#
x0 = 0#
x1 = 1.25#
x2 = 3.75#
dfx = fx0 *(2*x - x1 - x2)/((x0 - x1)*(x0 - x2)) + fx1 *(2*x - x0 - x2)/((x1 - x0)*(x1 - x2)) + fx2 *(2*x - x1 - x0)/((x2 - x1)*(x2 - x0))#
print "(dT/dz) = ",dfx,"C/cm"
q = - k * p *C * dfx*100#
print "q(z = 0) =",q,"W/m**2"
(dT/dz) =  -1.33333333333 C/cm
q(z = 0) = 70.56 W/m**2

Ex:23.4 Pg: 662

In [12]:
from sympy.mpmath import quad
from numpy import trapz,diff
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#
Q=quad(f,[0,0.8])
print "Q=",Q
x=[0, 0.12 ,0.22 ,0.32 ,0.36 ,0.4 ,0.44 ,0.54 ,0.64, 0.7 ,0.8]
y=[]
for xx in x:
    y.append(f(xx))
integral=trapz(x,y)
print "intergral=",integral
print "diff(x)=",diff(x)
d=diff(y)/diff(x)#
print "d=",d
Q= 1.64053333333333
intergral= -1.40920096
diff(x)= [ 0.12  0.1   0.1   0.04  0.04  0.04  0.1   0.1   0.06  0.1 ]
d= [  9.247744  -0.04488    4.38152    8.287744   9.527424   9.674624
   6.64312   -3.25368  -13.648816 -21.31    ]

Ex:23.5 Pg: 664

In [13]:
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#
Qt=1.640533#
Q=quad(f,[0,0.8])
print "Computed=",Q
Er=abs(Q-Qt)*100/Qt
print "Error estimate=",Er
Computed= 1.64053333333333
Error estimate= 2.03185997099678e-5