# Chapter 11 - Numerical differentiation¶

## Example No. 11_01 Pg No. 348¶

In :
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 :
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 :
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 :
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 :
#Differentiation of tabulated data

T = range(5,10)
s = [10,  14.5,  19.5,  25.5 , 32 ]
h = T-T
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 :
T = range(5,10)
s = [10,  14.5 , 19.5 , 25.5,  32 ]
h = T-T
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 :
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,C,C)

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 :
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