Chapter 11 - Numerical differentiation

Example No. 11_01 Pg No. 348

In [1]:
from __future__ import division
from scipy.misc import derivative

#First order forward difference

def f(x):
    F=x**2
    return F
                 

def df(x,h):
    DF=(f(x+h)-f(x))/h
    return DF

dfactual = derivative(f,1)
h = [0.2, 0.1, 0.05, 0.01 ]
print 'h\ty\terr'
for hh in h:
    y = df(1,hh)
    err = y - dfactual
    print hh,'\t',y,'\t',err

h	y	err
0.2 	2.2 	0.2
0.1 	2.1 	0.1
0.05 	2.05 	0.05
0.01 	2.01 	0.01

Example No. 11_02 Pg No. 350

In [2]:
from __future__ import division
from scipy.misc import derivative

#Three-Point Formula

def f(x):
    F = x**2
    return F
def df(x,h):
    DF = (f(x+h)-f(x-h))/(2*h)
    return DF
dfactual = derivative(f,1)
h = [0.2, 0.1, 0.05, 0.01 ]
print 'h\ty\terr'
for hh in h:
    y = df(1,hh)
    err = y - dfactual
    print hh,'\t',y,'\t',err

h	y	err
0.2 	2.0 	-4.4408920985e-16
0.1 	2.0 	4.4408920985e-16
0.05 	2.0 	4.4408920985e-16
0.01 	2.0 	1.7763568394e-15

Example No. 11_03 Pg No. 353

In [1]:
from __future__ import division
from scipy.misc import derivative
from numpy import arange,sin,cos,sqrt
x = 0.45#
def f(x):
    F = sin(x)
    return F
def df(x,h):
    DF = (f(x+h) - f(x))/h
    return DF
dfactual =  cos(x)#
h = arange(0.01,0.005+0.04,0.005)
n = len(h)#
print 'h\ty\t\terr'
for hh in h:
    y = df(x,hh)
    err = y - dfactual 
    print hh,'\t',y,'\t',err

#using 16 significant digits so the bound for roundoff error is 0.5*10**(-16)
e = 0.5*10**(-16)
M2 = max(sin(x+h))#
hopt = 2*sqrt(e/M2)#
print '\nM2 = %.3e & hopt = %.3e'%(hopt,M2)

h	y		err
0.01 	0.898257285429 	-0.00218981692372
0.015 	0.897151155627 	-0.00329594672573
0.02 	0.896037563392 	-0.00440953896077
0.025 	0.894916522709 	-0.00553057964367
0.03 	0.893788047675 	-0.00665905467759
0.035 	0.892652152498 	-0.00779494985432
0.04 	0.891508851498 	-0.00893825085448

M2 = 2.061e-08 & hopt = 4.706e-01

Example No. 11_04 Pg No. 354

In [2]:
from __future__ import division
from numpy import arange,sin,cos,sqrt

#Approximate second derivative

x = 0.75#
h = 0.01#
def f(x):
    F = cos(x)
    return F
def d2f(x,h):
    D2F = ( f(x+h) - 2*f(x) + f(x-h) )/h**2
    return D2F
y = d2f(0.75,0.01)#
d2fexact = -cos(0.75)
err = d2fexact - y #
print ' d2fexact = %.3f \n err = %.3e \n y = %.3f '%(d2fexact,err,y)

 d2fexact = -0.732 
 err = -6.097e-06 
 y = -0.732

Example No. 11_05 Pg No. 358

In [1]:
#Differentiation of tabulated data

T = range(5,10)
s = [10,  14.5,  19.5,  25.5 , 32 ]
h = T[1]-T[0]
n = len(T)
def v(t):
    if t in T and T.index(t) == 0:
        V = (-3*s[T.index(t) ] + 4*s[T.index(t+h)]  - s[T.index(t+2*h) ] ) /(2*h) #Three point forward difference formula
    elif t in T and T.index(t) == n-1:
        V = ( 3*s[T.index(t)] - 4*s[ T.index(t-h)]  + s[T.index(t-2*h) ]  )/(2*h) #Backward difference formula
    else:
        V = ( s[T.index(t+h)]  - s[T.index(t-h)]  )/(2*h) #Central difference formula
    return V

v_5 = v(5)
v_7 = v(7)
v_9 = v(9)

print 'v(5) = ',v_5
print 'v(7) = ',v_7
print 'v(9) = ',v_9

v(5) =  4.25
v(7) =  5.5
v(9) =  6.75

Example No. 11_06 Pg No. 359

In [2]:
T = range(5,10)
s = [10,  14.5 , 19.5 , 25.5,  32 ]
h = T[1]-T[0]
def a(t):
    A = ( s[T.index(t+h) ] - 2*s[ T.index(t) ] + s[T.index(t-h)  ] )/h**2
    return A
a_7 = a(6)

print 'a(7) = ',a_7

a(7) =  0.5

Example No. 11_7 Pg No. 359

In [3]:
from numpy import mat,divide,linalg
h = 0.25 #
#y''(x) = e**(x**2)
#y(0) = 0 , y(1) = 0
# y''(x) = y(x+h) - 2*y(x) + y(x-h)/h**2  = e**(x**2)
#( y(x + 0.25) - 2*y(x) + y(x-0.25) )/0.0625 = e**(x**2)
#y(x+0.25) - 2*y(x) + y(x - 0.25) = 0.0624*e**(x**2)
#y(0.5) - 2*y(0.25) + y(0) = 0.0665
#y(0.75) - 2*y(0.5) + y(0.25) = 0.0803
#y(1) - 2*y(0.75) + y(0.5) = 0.1097
#given y(0) = y(1) = 0
#
#0 + y2 - 2y1 = 0.06665
#y3 - 2*y2 + y1 = 0.0803
#-2*y3 + y2 + 0 = 0.1097
#Therefore
A = mat([[0, 1, -2],[1, -2, 1],[-2, 1, 0 ]])
B = mat([[ 0.06665],[0.0803] ,[0.1097 ]])
C = linalg.solve(A,B)
print 'solving the above equations we get \n\n y1 = y(0.25) = %f \n y2 = y(0.5) = %f \n y3 = y(0.75) = %f \n '%(C[2],C[1],C[0])

solving the above equations we get 

 y1 = y(0.25) = -0.117562 
 y2 = y(0.5) = -0.168475 
 y3 = y(0.75) = -0.139088

Example No. 11_08 Pg No. 362

In [1]:
from numpy import arange,exp
from scipy.misc import derivative

x = arange(-0.5,0.25+1.5,0.25)
h = 0.5 ;
r = 1.0/2 ;

def f(x):
    F = exp(x)
    return F
def D(x,h):
    D = (f(x + h) - f(x-h) )/(2*h) 
    return D
def df(x,h,r):
    Df = (D(x,r*h) - r**2*D(x,h))/(1-r**2)
    return Df

df_05 = df(0.5,0.5,1.0/2);
print 'richardsons technique -\n\ndf(0.5) = ',df_05
print 'D(rh) = D(0.25) = ',D(0.5,0.5)
print 'D(0.5) = ',D(0.5,0.25)

dfexact = derivative(f,0.5)
print 'Exact df(0.5) = ',dfexact
print 'The result by richardsons technique is much better than other results'

#r = 2
print 'for r = 2'
print 'df(x) = ',df(0.5,0.5,2)
print 'D(rh) = ',D(0.5,2*0.5)

richardsons technique -

df(0.5) =  1.64850499031
D(rh) = D(0.25) =  1.71828182846
D(0.5) =  1.66594919985
Exact df(0.5) =  1.93757920531
The result by richardsons technique is much better than other results
for r = 2
df(x) =  1.64518270284
D(rh) =  1.93757920531