# Chapter 13 - Numerical solution of ordinary differential equations¶

## Example No. 13_01 Pg No. 414¶

In [1]:
from __future__ import division
from math import sqrt
#Taylor method

def f(x,y):
F = x**2 + y**2
return F
def d2y(x,y):
D2Y = 2*x + 2*y*f(x,y)
return D2Y
def d3y(x,y):
D3Y = 2 + 2*y*d2y(x,y) + 2*f(x,y)**2
return D3Y
def y(x):
Y = 1 + f(0,1)*x + d2y(0,1)*x**2/2 + d3y(0,1)*sqrt(x)
return Y
print ' y(0.25) = ',y(0.25)
print ' y(0.5) = ',y(0.5)

 y(0.25) =  5.3125
y(0.5) =  7.40685424949


## Example No. 13_02 Pg No. 415¶

In [1]:
from __future__ import division
from scipy.misc import factorial
#Recursive Taylor Method

def f(x,y):
F = x**2 + y**2
return F
def d2y(x,y):
D2Y = 2*x + 2*y*f(x,y)
return D2Y
def d3y(x,y):
D3Y = 2 + 2*y*d2y(x,y) + 2*f(x,y)**2
return D3Y
def d4y(x,y):
D4Y = 6*f(x,y)*d2y(x,y) + 2*y*d3y(x,y)
return D4Y
h = 0.2 #
def y(x,y):
Y = y + f(x,y)*h + d2y(x,y)*h**2/2 + d3y(x,y)*h**3/6 + d4y(x,y)*h**4/factorial(4)
return Y
x0 = 0#
y0 = 0 #
y_=[]
for i in range(1,3):
y_.append(y(x0,y0))
print ' Iteration-%d\n\n dy(0) = %f\n d2y(0) = %f\n d3y(0) = %f\n d4y(0) = %f\n '%(i,f(x0,y0),d2y(x0,y0),d3y(x0,y0),d4y(x0,y0))
x0 = x0 + i*h
y0 = y_[i-1]
print 'y(0) = %f\n\n'%(y_[i-1])

 Iteration-1

dy(0) = 0.000000
d2y(0) = 0.000000
d3y(0) = 2.000000
d4y(0) = 0.000000

y(0) = 0.002667

Iteration-2

dy(0) = 0.040007
d2y(0) = 0.400213
d3y(0) = 2.005336
d4y(0) = 0.106763

y(0) = 0.021353



## Example No. 13_03 Pg No. 417¶

In [2]:
from __future__ import division
from math import exp
#Picard's Method

#y'(x) = x**2 + y**2,y(0) = 0
#y(1) = y0 + integral(x**2 + y0**2,x0,x)
#y(1) = x**3/3
#y(2) = 0 + integral(xY2 + y1**2,x0,x)
#     = integral(x**2 + x**6/9,0,x) = x**3/3 + x**7/63
#therefore y(x) = x**3 /3 + x**7/63
def y(x):
Y = x**3/3 + x**7/63
return Y
print 'for dy(x) = x**2 + y**2 the results are :'
print 'y(0.1) = ',y(0.1)
print 'y(0.2) = ',y(0.2)
print 'y(1) = ',y(1)

#y'(x) = x*e**y, y(0) = 0
#y0 = 0 , x0 = 0
#Y(1) = 0 + integral(x*e**0,0,x) = x**2/2
#y(2) = 0 + integral( x*e**( x**2/2 ) ,0,x) = e**(x**2/2)-1
#therefore y(x) = e**(x**2/2) - 1
def y(x):
Y = exp(x**2/2) - 1
return Y
print 'for dy(x) = x*e**y the results are '
print 'y(0.1) = ',y(0.1)
print 'y(0.2) = ',y(0.2)
print 'y(1) = ',y(1)

for dy(x) = x**2 + y**2 the results are :
y(0.1) =  0.000333334920635
y(0.2) =  0.00266686984127
y(1) =  0.349206349206
for dy(x) = x*e**y the results are
y(0.1) =  0.0050125208594
y(0.2) =  0.0202013400268
y(1) =  0.6487212707


## Example No. 13_04 Pg No. 417¶

In [3]:
from __future__ import division

#Euler's Method

def dy(x):
DY = 3*x**2 + 1
return DY
x0 = 1
y = [0,2] #
#h = 0.5
h = 0.5
print 'for h = %f\n'%(h)
for i in range(2,4):
y.append(y[i-1] + h*dy(x0+(i-2)*h))
print 'y(%f) = %f\n'%(x0+(i-1)*h,y[i])

#h = 0.25
h = 0.25
print '\nfor h = %f\n'%(h)
y = [0,2] #
for i in range(2,6):
y.append(y[i-1] + h*dy(x0+(i-2)*h))
print 'y(%f) = %f\n'%(x0+(i-1)*h,y[i])

for h = 0.500000

y(1.500000) = 4.000000

y(2.000000) = 7.875000

for h = 0.250000

y(1.250000) = 3.000000

y(1.500000) = 4.421875

y(1.750000) = 6.359375

y(2.000000) = 8.906250



## Example No. 13_05 Pg No. 422¶

In [1]:
from numpy import nditer
from __future__ import division
#Error estimation in Euler's Method

def dy(x):
DY = 3*x**2 + 1
return DY
def d2y(x):
D2Y = 6*x
return D2Y
def d3y(x):
D3Y = 6
return D3Y
def exacty(x):
Exacty = x**3 + x
return Exacty
x0 = 1
y = 2
h = 0.5
X=[];ESTY=[];TRUEY=[];ET=[];GERR=[];
for i in range(2,4):
x = x0 + (i-1)*h
y = y + h*dy(x0+(i-2)*h)
print '\n Step %d \n x(%d) = %f\n y(%f) = %f\n'%(i-1,i-1,x,x,y)
Et = d2y(x0+(i-2)*h)*h**2/2 +  d3y(x0+(i-2)*h)*h**3/6
print '\n Et(%d) = %f\n'%(i-1,Et)
truey = exacty(x0+(i-1)*h)
gerr = truey - y
X.append(x)
ESTY.append(y)
TRUEY.append(truey)
ET.append(Et)
GERR.append(gerr)

#table = [x y(2:3) truey Et gerr]
print 'x      Est y   true y    Et      Globalerr'
for a,b,c,d,e in nditer([X,ESTY,TRUEY,ET,GERR]):
print a,'\t',b,'\t',c,'\t',d,'\t',e

 Step 1
x(1) = 1.500000
y(1.500000) = 4.000000

Et(1) = 0.875000

Step 2
x(2) = 2.000000
y(2.000000) = 7.875000

Et(2) = 1.250000

x      Est y   true y    Et      Globalerr
1.5 	4.0 	4.875 	0.875 	0.875
2.0 	7.875 	10.0 	1.25 	2.125


## Example No. 13_06 Pg No. 427¶

In [1]:
from numpy import nditer,sqrt
from __future__ import division

#Heun's Method

def f(x,y):
F = 2*y/x
return F
def exacty(x):
Exacty = 2*sqrt(x)
return Exacty
x = [0,1] #
y = [0,2] #
h = 0.25 #
#Euler's Method
print 'EULERS METHOD'
for i in range(2,6):
x.append(x[i-1] + h )
y.append(y[i-1] + h*f(x[i-1],y[i-1]))
print 'y(%f) = %f \n '%(x[i],y[i])

eulery = y
#Heun's Method
print 'HEUNS METHOD'
ye=[0,0]
y = [0,2] #
for i in range(2,6):
m1 = f(x[i-1],y[i-1]) #
ye.append(y[i-1] + h*f(x[(i-1)],y[(i-1)]))
m2 = f(x[(i)],ye[(i)])
y.append(y[(i-1)] + h*(m1 + m2)/2)
print '\nIteration %d \n m1 = %f\n ye(%f) = %f \n m2 = %f \n y(%f) = %f \n'%(i-1,m1,x[(i)],ye[(i)],m2,x[(i)],y[(i)])

truey = exacty(x)
print 'x     Eulers   Heuns   \t\tAnalytical'
for a,b,c,d in nditer([x,eulery,y,truey]):
print a,'\t',b,'\t%.6f'%c,'\t',d


EULERS METHOD
y(1.250000) = 3.000000

y(1.500000) = 4.200000

y(1.750000) = 5.600000

y(2.000000) = 7.200000

HEUNS METHOD

Iteration 1
m1 = 4.000000
ye(1.250000) = 3.000000
m2 = 4.800000
y(1.250000) = 3.100000

Iteration 2
m1 = 4.960000
ye(1.500000) = 4.340000
m2 = 5.786667
y(1.500000) = 4.443333

Iteration 3
m1 = 5.924444
ye(1.750000) = 5.924444
m2 = 6.770794
y(1.750000) = 6.030238

Iteration 4
m1 = 6.891701
ye(2.000000) = 7.753163
m2 = 7.753163
y(2.000000) = 7.860846

x     Eulers   Heuns   		Analytical
0.0 	0.0 	0.000000 	0.0
1.0 	2.0 	2.000000 	2.0
1.25 	3.0 	3.100000 	2.2360679775
1.5 	4.2 	4.443333 	2.44948974278
1.75 	5.6 	6.030238 	2.64575131106
2.0 	7.2 	7.860846 	2.82842712475


## Example No. 13_07 Pg NO. 433¶

In [1]:
from __future__ import division

#Polygon Method

def f(x,y):
F = 2*y/x
return F
x=[1]
y=[2]
h = 0.25 #
for i in range(1,3):
x.append(x[(i-1)] + h )
y.append(y[(i-1)] + h*f(  x[(i-1)]+ h/2  ,  y[(i-1)] + h*f( x[(i-1)]  ,  y[(i-1)] )/2  ))
print 'y(%f) = %f \n '%(x[(i)],y[(i)])

y(1.250000) = 3.111111

y(1.500000) = 4.468687



## Example No. 13_08 Pg No. 439¶

In [1]:
from __future__ import division

#Classical Runge Kutta Method

def f(x,y):
F = x**2 + y**2
return F
h = 0.2
x=[0]
y=[0]

for i in range(0,2):
m1 = f(  x[(i)]  ,  y[(i)]  ) #
m2 = f(  x[i] + h/2  ,  y[(i)] + m1*h/2  ) #
m3 = f(  x[(i)] + h/2  ,  y[(i)] + m2*h/2  ) #
m4 = f(  x[(i)] + h  ,  y[(i)] + m3*h  ) #
x.append(x[(i)] + h )
y.append(y[(i)] +  (m1 + 2*m2 + 2*m3 + m4)*h/6 )

print '\nIteration - %d\n m1 = %f\n m2 = %f \n m3 = %f \n m4 = %f \n y(%f) = %f \n'%(i+1,m1,m2,m3,m4,x[(i+1)],y[(i+1)])

Iteration - 1
m1 = 0.000000
m2 = 0.010000
m3 = 0.010001
m4 = 0.040004
y(0.200000) = 0.002667

Iteration - 2
m1 = 0.040007
m2 = 0.090044
m3 = 0.090136
m4 = 0.160428
y(0.400000) = 0.021360



## Example No. 13_09 Pg No. 444¶

In [1]:
from __future__ import division

#Optimum step size

x = 0.8 #
h1 = 0.05 #
y1 = 5.8410870 #
h2 = 0.025 #
y2 = 5.8479637 #

#d = 4
h = ((h1**4 - h2**4)*10**(-4)/(2*(y2 - y1)))**(1/4)
print 'for four decimal places'
print 'h = ',h

#d = 6
h = ((h1**4 - h2**4)*10**(-6)/(2*(y2 - y1)))**(1/4)
print 'for six decimal places'
print 'h = ' ,h
print 'Note-We can use h = 0.01 for four decimal places and h = 0.004 for six decimal places'

for four decimal places
h =  0.0143668085653
for six decimal places
h =  0.00454318377739
Note-We can use h = 0.01 for four decimal places and h = 0.004 for six decimal places


## Example No. 13_10 Pg NO. 446¶

In [1]:
from __future__ import division

#Milne-Simpson Predictor-Corrector method

def f(x,y):
F = 2*y/x
return F
x0 = 1 #
y0 = 2 #
h = 0.25 #
#Assuming y1 ,y2 and y3(required for milne-simpson formula) are estimated using Fourth- order Runge kutta method
x1 = x0 + h
y1 = 3.13 #
x2 = x1 + h
y2 = 4.5 #
x3 = x2 + h
y3 = 6.13 #
#Milne Predictor formula
yp4 = y0 + 4*h*(2*f(x1,y1) - f(x2,y2) + 2*f(x3,y3))/3
x4 = x3 + h
fp4 = f(x4,yp4) #
print 'yp4 = ',yp4
print 'fp4 = ',fp4,
#Simpson Corrector formula
yc4 = y2 + h*( f(x2,y2) + 4*f(x3,y3) + fp4)/3
f4 = f(x4,yc4)
print 'yc4 = ',yc4
print 'f4 = ',f4

yc4 = y2 + h*( f(x2,y2) + 4*f(x3,y3) + f4)/3
print 'yc4 = ',yc4
print 'Note- the exact solution is y(2) = 8'

yp4 =  8.00914285714
fp4 =  8.00914285714 yc4 =  8.00266666667
f4 =  8.00266666667
yc4 =  8.00212698413
Note- the exact solution is y(2) = 8


## Example No. 13_11 Pg NO. 446¶

In [1]:
from __future__ import division

def f(x,y):
F = 2*y/x
return F
x0 = 1 #
y0 = 2 #
h = 0.25 #
x1 = x0 + h
y1 = 3.13 #
x2 = x1 + h
y2 = 4.5 #
x3 = x2 + h
y3 = 6.13 #
yp4 = y3 + h*(55*f(x3,y3) - 59*f(x2,y2) + 37*f(x1,y1) - 9*f(x0,y0))/24
x4 = x3 + h
fp4 = f(x4,yp4)
print 'yp4 = ',yp4
print 'fp4 = ',fp4
yc4 = y3 + h*( f(x1,y1) - 5*f(x2,y2) + 19*f(x3,y3) + 9*fp4 )/24
f4 = f(x4,yc4)
print 'yc4 = ',yc4
print 'f4 = ',f4

yc4 = y3 + h*( f(x1,y1) - 5*f(x2,y2) + 19*f(x3,y3) + 9*f4 )/24
print '\nrefined-yc4 = ',yc4

Adams Predictor formula
yp4 =  8.01135714286
fp4 =  8.01135714286

yc4 =  8.00727901786
f4 =  8.00727901786

refined-yc4 =  8.00689669364


## Example No. 13_12 Pg No. 453¶

In [1]:
from __future__ import division

#Milne-Simpson Method using modifier

def f(y):
F = -y**2
return F
x = [ 1 , 1.2 , 1.4 , 1.6 ]
y = [ 1 , 0.8333333 , 0.7142857 , 0.625 ]
h = 0.2 #

for i in range(0,2):
yp = y[(i)] + 4*h*( 2*f( y[(i+1)] ) - f( y[(i+2)] ) + 2*f( y[(i+3)] ) )/3
fp = f(yp) #
yc = y[( i+2)] + h*(f( y[(i+2)] ) + 4*f( y[(i+3)] ) + fp )/3
Etc = -(yc - yp)/29
y.append(yc + Etc)
print '\n y%dp = %f\n f%dp = %f \n y%dc = %f \n Modifier Etc = %f \n Modified y%dc = %f \n'%(i+4,yp,i+4,fp,i+4,yc,Etc,i+4,y[(i+4)])
err = exactanswer - y[5] #
print 'error = ',err

 y4p = 0.557351
f4p = -0.310640
y4c = 0.555396
Modifier Etc = 0.000067
Modified y4c = 0.555464

y5p = 0.500837
f5p = -0.250837
y5c = 0.499959
Modifier Etc = 0.000030
Modified y5c = 0.499989

error =  1.11332736336e-05


## Example No. 13_13 Pg No. 455¶

In [1]:
from __future__ import division

#System of differetial Equations

def f1(x,y1,y2):
F1 = x + y1 + y2
return F1
def f2(x,y1,y2):
F2 = 1 + y1 + y2
return F2

x0 = 0 #
y10 = 1 #
y20 = -1 #
h = 0.1 #
m1= [f1( x0 , y10 , y20 )]
m1.append(f2( x0 , y10 , y20 ))
m2=[f1( x0+h , y10 + h*m1[0] , y20 + h*m1[1] )]
m2.append(f2( x0+h , y10 + h*m1[0] , y20 + h*m1[1] ))
m= [(m1[0] + m2[0])/2 ]
m.append((m1[1] + m2[1])/2)

y1_0_1 = y10 + h*m[0]
y2_0_1 = y20 + h*m[1]

print 'm1(1) = %f\n m1(2) = %f\n m2(1) = %f\n m2(2) = %f\n m(1) = %f\n m(2) = %f\n y1(0.1) = %f\n y2(0.1) = %f\n'%(m1[0],m1[1],m2[0],m2[1],m[0],m[1],y1_0_1,y2_0_1)

m1(1) = 0.000000
m1(2) = 1.000000
m2(1) = 0.200000
m2(2) = 1.100000
m(1) = 0.100000
m(2) = 1.050000
y1(0.1) = 1.010000
y2(0.1) = -0.895000



## Example No. 13_14 Pg No. 457¶

In [1]:
from __future__ import division
#Higher Order Differential Equations

x0 = 0
y10 = 0
y20 = 1
h = 0.2
m1 = [y20] #
m1.append(6*x0 + 3*y10 - 2*y20)
m2= [y20 + h*m1[1]]
m2.append(6*(x0+h) + 3*(y10 + h*m1[0]) - 2*(y20 + h*m1[1]) )
m = [(m1[0] + m2[0])/2 ]
m.append((m1[1] + m2[1])/2)

y1_0_2 = y10 + h*m[0]
y2_0_2 = y20 + h*m[1]

print 'm1(1) = %f\n m1(2) = %f\n m2(1) = %f\n m2(2) = %f\n m(1) = %f\n m(2) = %f\n y1(0.1) = %f\n y2(0.1) = %f\n'%(m1[0],m1[1],m2[0],m2[1],m[0],m[1],y1_0_2,y2_0_2)

m1(1) = 1.000000
m1(2) = -2.000000
m2(1) = 0.600000
m2(2) = 0.600000
m(1) = 0.800000
m(2) = -0.700000
y1(0.1) = 0.160000
y2(0.1) = 0.860000