Chapter 27 : Numerical Solution Of Ordinary Differential Equations

Example 27.1, page no. 686

In [1]:
import numpy,sympy,math

x = sympy.Symbol('x')
print "Solution through picards method "
n = int(raw_input("The no of iterations required")) 
print "y(0)=1 and y(x)=x+y "
yo = 1
yn = 1
for i in range(1,n+1):
    yn = yo+sympy.integrate(yn+x,(x,0,x))
print "Y = ",yn

Solution through picards method 
The no of iterations required3
y(0)=1 and y(x)=x+y 
Y =  x**4/24 + x**3/3 + x**2 + x + 1

Example 27.2, page no. 687

In [5]:
import numpy,sympy,math

x = sympy.Symbol('x')
print "Solution through picards method "
n = int(raw_input("The no of iterations required : ")) 
print "y(0)=1 and y(x)=x+y "
yo = 1
y = 1
for i in range(1,n+1):
    f = (y-x)/(y+x) 
    y = yo+sympy.integrate(f,(x,0,x))
x = 0.1
print "Y = ",y.evalf()

Solution through picards method 
The no of iterations required : 2
y(0)=1 and y(x)=x+y 
Y =  -Integral(2*x/(2*log(x + 1) + 1), x) + Integral(2*x/(2*log(x + 1) + 1), (x, 0)) - Integral(-2*log(x + 1)/(2*log(x + 1) + 1), x) + Integral(-2*log(x + 1)/(2*log(x + 1) + 1), (x, 0)) - Integral(-1/(2*log(x + 1) + 1), x) + Integral(-1/(2*log(x + 1) + 1), (x, 0)) + 1.0

Example 27.5, page no. 690

In [10]:
print "Solution using Eulers Method "
print x
print y
n = int(raw_input("Input the number of iteration :− "))
x = 0
y = 1
for i in range(1,n+1):
    y1 = x+y
    y = y+0.1*y1
    x = x+0.1
print "The value of y is :− ",y

Solution using Eulers Method 
0.3
1.362
Input the number of iteration :− 4
The value of y is :−  1.5282

Example 27.6, page no. 691

In [11]:
print "Solution using Eulers Method "
print x
print y
n = int(raw_input("Input the number of iteration :− "))
x = 0
y = 1
for i in range(1,n+1):
    y1 = (y-x)/(y+x)
    y = y+0.02*y1
    x = x+0.1
    print y
print "The value of y is :− ",y

Solution using Eulers Method 
0.4
1.5282
Input the number of iteration :− 2
1.02
1.03642857143
The value of y is :−  1.03642857143

Example 27.7, page no. 692

In [12]:
print "Solution using Eulers Method "
print x
print y
n = int(raw_input("Input the number of iteration :− "))
x = 0.1
m = 1
y = 1
yn = 1
y1 = 1
k = 1
for i in range(1,n+1):
    yn = y 
    for i in range(1,5):
        m = (k+y1)/2
        yn = y+0.1*m
        y1 = (yn+x)
        print yn
    print "−−−−−−−−−−−−−−−−−−−−−−− "
    y = yn 
    m = y1
    yn = yn+0.1*m 
    print yn
    x = x+0.1
    yn = y 
    k = m
print "The value of y is :− ",y

Solution using Eulers Method 
0.2
1.03642857143
Input the number of iteration :− 3
1.1
1.11
1.1105
1.110525
−−−−−−−−−−−−−−−−−−−−−−− 
1.2315775
1.2315775
1.242630125
1.24318275625
1.24321038781
−−−−−−−−−−−−−−−−−−−−−−− 
1.38753142659
1.38753142659
1.39974747853
1.40035828113
1.40038882126
−−−−−−−−−−−−−−−−−−−−−−− 
1.57042770339
The value of y is :−  1.40038882126

Example 27.8, page no. 694

In [14]:
import math

print "Solution using Eulers Method "
print x
print y
n = int(raw_input("Input the number of iteration :− "))
x = 0.2
m = 0.301
y = 2
yn = 2
y1 = math.log10(2)
k = 0.301
for i in range(1,n+1):
    yn = y
    for i in range(1,5):
        m = (k+y1)/2
        yn = y+0.2*m
        y1 = math.log10(yn+x)
        print yn
    print "−−−−−−−−−−−−−−−−−−−−−−−"
    y = yn
    m = y1
    yn = yn+0.2*m 
    print yn
    x = x+0.2
    yn = y
    k = m
print "The value of y is :− ",y

Solution using Eulers Method 
0.2
2
Input the number of iteration :− 3
2.06020299957
2.06551474469
2.06561668931
2.06561864352
−−−−−−−−−−−−−−−−−−−−−−−
2.13665600548
2.13665600548
2.14156348217
2.14164742068
2.14164885497
−−−−−−−−−−−−−−−−−−−−−−−
2.22267196493
2.22267196493
2.22722645093
2.22729646948
2.22729754504
−−−−−−−−−−−−−−−−−−−−−−−
2.31757184825
The value of y is :−  2.22729754504

Example 27.9, page no. 696

In [15]:
import math

print "Solution using Eulers Method "
print x
print y
n = int(raw_input("Input the number of iteration :− "))
x = 0.2 
m = 1
y = 1
yn = 1
y1 = 1
k = 1
for i in range(1,n+1):
    yn = y
    for i in range(1,5):
        m = (k+y1)/2
        yn = y+0.2*m 
        y1 = math.sqrt(yn)+x
        print yn
    print "−−−−−−−−−−−−−−−−−−−−−−− "
    y = yn
    m = y1
    yn = yn+0.2*m
    print yn
    x = x+0.2
    yn = y
    k = m
print "The value of y is :− ",y

Solution using Eulers Method 
0.8
2.22729754504
Input the number of iteration :− 2
1.2
1.2295445115
1.23088482816
1.23094524903
−−−−−−−−−−−−−−−−−−−−−−− 
1.49284119302
1.49284119302
1.52407510156
1.52534665765
1.52539814634
−−−−−−−−−−−−−−−−−−−−−−− 
1.85241216588
The value of y is :−  1.52539814634

Example 27.10, page no. 697

In [1]:
print "Runges method"
def f(x,y):
    y = x+y
    return y
x = 0
y = 1
h = 0.2
k1 = h*f(x,y)
k2 = h*f(x+1/2*h,y+1/2*k1)
kk = h*f(x+h,y+k1)
k3 = h*f(x+h,y+kk)
k = 1/6*(k1+4*k2+k3)
print "The required approximate value is :− "
y = y+k
print y

Runges method
The required approximate value is :− 
1.0

Example 27.11, page no. 697

In [2]:
print "Runges method"
def f(x,y):
    y = x+y
    return y
x = 0
y = 1
h = 0.2
k1 = h*f(x,y)
k2 = h*f(x+1/2*h,y+1/2*k1)
k3 = h*f(x+1/2*h,y+1/2*k2)
k4 = h*f(x+h,y+k3)
k = 1/6*(k1+2*k2+2*k3+k4)
print "The required approximate value is :− "
y = y+k
print y

Runges method
The required approximate value is :− 
1.0

Example 27.12, page no. 698

In [3]:
print "Runga kutta method"
def f(x,y):
    y = (y**2-x**2)/(x**2+y**2)
    return y
x = 0
y = 1
h = 0.2
k1 = h*f(x,y)
k2 = h*f(x+1/2*h,y+1/2*k1)
k3 = h*f(x+1/2*h,y+1/2*k2)
k4 = h*f(x+h,y+k3)
k = 1/6*(k1+2*k2+2*k3+k4)
print "The required approximate value is :− "
y = y+k
print y
print "To find y(0.4) put x=0.2 y=above value ie 1.196 h=0.2 "
x = 0.2
h = 0.2
k1 = h*f(x,y)
k2 = h*f(x+1/2*h,y+1/2*k1)
k3 = h*f(x+1/2*h,y+1/2*k2)
k4 = h*f(x+h,y+k3)
k = 1/6*(k1+2*k2+2*k3+k4)
print "The required approximate value is :− "
y = y+k
print y

Runga kutta method
The required approximate value is :− 
1.0
To find y(0.4) put x=0.2 y=above value ie 1.196 h=0.2 
The required approximate value is :− 
1.0

Example 27.13, page no. 699

In [4]:
print "Runga kutta method"
def f(x,y):
    yy = x+y**2
    return yy
x = 0
y = 1
h = 0.1
k1 = h*f(x,y)
k2 = h*f(x+1/2*h,y+1/2*k1)
k3 = h*f(x+1/2*h,y+1/2*k2) 
k4 = h*f(x+h,y+k3)
k = 1/6*(k1+2*k2+2*k3+k4)
print "The required approximate value is :− "
y = y+k
print y
print "To find y(0.4) put x=0.2 y=above value ie 1.196 h=0.2 "
x = 0.1
h = 0.1
k1 = h*f(x,y)
k2 = h*f(x+1/2*h,y+1/2*k1)
k3 = h*f(x+1/2*h,y+1/2*k2)
k4 = h*f(x+h,y+k3)
k = 1/6*(k1+2*k2+2*k3+k4)
print "The required approximate value is :− "
y = y+k
print y

Runga kutta method
The required approximate value is :− 
1.0
To find y(0.4) put x=0.2 y=above value ie 1.196 h=0.2 
The required approximate value is :− 
1.0

Example 27.14, page no. 700

In [30]:
import sympy

x = sympy.Symbol('x')
yo = 0
y = 0
h = 0.2
f = x-y**2
y = sympy.integrate(f,(x,0,x))
y1 = yo+y
print "y1= ",y1
f = x-y**2
y = sympy.integrate(f,(x,0,x))
y2 = yo+y
print "y2= ",y2
y = x-y**2
y = sympy.integrate(f,(x,0,x))
y3 = yo+y
print "y3= ",y3
print "Determining the initial values for milnes method using y3 "
print "x=0.0 y0=0.0 f0=0 "
print "x=0.2 y1= ",
x = 0.2
print y1
#y1 = sympy.solve(y1)
print "f1=",
f1 = x-y1**2
print f1
print "x=0.4 y2= ",
x = 0.4
print y2
print "f2= ",
f2 = x-y2**2
print f2
print "x=0.6 y3= ",
x = 0.6
print y3
print "f3=",
f3 = x-y3**2
print f3
print "Using predictor method to find y4"
x = 0.8
y4 = yo+4/3*h*(2*f1-f2+2*f3)
print "y4=",y4
f4 = x-y**2
print "f4= ",f4
print "Using predictor method to find y5 "
x = 1.0
y5 = sympy.solve(y1+4/3*h*(2*f2-f3+2*f4))
print y5
f5 = sympy.solve(x-y**2)
print "f5 =",f5
print "Hence y(1)= ",y5

y1=  x**2/2
y2=  -x**5/20 + x**2/2
y3=  -x**5/20 + x**2/2
Determining the initial values for milnes method using y3 
x=0.0 y0=0.0 f0=0 
x=0.2 y1=  x**2/2
f1= -x**4/4 + 0.2
x=0.4 y2=  -x**5/20 + x**2/2
f2=  -(-x**5/20 + x**2/2)**2 + 0.4
x=0.6 y3=  -x**5/20 + x**2/2
f3= -(-x**5/20 + x**2/2)**2 + 0.6
Using predictor method to find y4
y4= -0.1*x**4 - 0.2*(-x**5/20 + x**2/2)**2 + 0.24
f4=  -(-x**5/20 + x**2/2)**2 + 0.8
Using predictor method to find y5 
[-1.53353211362738, 2.59188514369719, -1.51825864909117 - 1.96576535664013*I, -1.51825864909117 + 1.96576535664013*I, -0.611593620665791 - 2.14551655586226*I, -0.611593620665791 + 2.14551655586226*I, 0.00495892544779777 - 0.780181689871346*I, 0.00495892544779777 + 0.780181689871346*I, 1.59571682927425 - 0.8271211188914*I, 1.59571682927425 + 0.8271211188914*I]
f5 = [-1.28459666867223, 2.38248090853886, -1.33416685318757 - 1.68618279990398*I, -1.33416685318757 + 1.68618279990398*I, -1.00716479884721 - 2.10036734766602*I, -1.00716479884721 + 2.10036734766602*I, 0.142926398918135 - 1.33989714259572*I, 0.142926398918135 + 1.33989714259572*I, 1.64946313318333 - 0.385565643189579*I, 1.64946313318333 + 0.385565643189579*I]
Hence y(1)=  [-1.53353211362738, 2.59188514369719, -1.51825864909117 - 1.96576535664013*I, -1.51825864909117 + 1.96576535664013*I, -0.611593620665791 - 2.14551655586226*I, -0.611593620665791 + 2.14551655586226*I, 0.00495892544779777 - 0.780181689871346*I, 0.00495892544779777 + 0.780181689871346*I, 1.59571682927425 - 0.8271211188914*I, 1.59571682927425 + 0.8271211188914*I]

Example 27.15, page no. 704

In [31]:
print "Runga kutta method"
def f(x,y):
    yy = x*y+y**2
    return yy
y0 = 1
x = 0
y = 1
h = 0.1
k1 = h*f(x,y)
k2 = h*f(x+1/2*h,y+1/2*k1)
k3 = h*f(x+1/2*h,y+1/2*k2)
k4 = h*f(x+h,y+k3)
ka = 1/6*(k1+2*k2+2*k3+k4)
print "The required approximate value is :− "
y1 = y+ka 
y = y+ka 
print y
print "To find y(0.4) put x=0.2 y=above value ie 1.196 h=0.2 "
x = 0.1
h = 0.1
k1 = h*f(x,y)
k2 = h*f(x+1/2*h,y+1/2*k1)
k3 = h*f(x+1/2*h,y+1/2*k2)
k4 = h*f(x+h,y+k3)
kb = 1/6*(k1+2*k2+2*k3+k4)
print "The required approximate value is :− "
y2 = y+kb 
y = y+kb
print y
print "To find y(0.4) put x=0.2 y=above value ie 1.196 h=0.2 "
x = 0.2
h = 0.1
k1 = h*f(x,y)
k2 = h*f(x+1/2*h,y+1/2*k1)
k3 = h*f(x+1/2*h,y+1/2*k2)
k4 = h*f(x+h,y+k3)
kc = 1/6*(k1+2*k2+2*k3+k4)
print "The required approximate value is :− "
y3 = y+kc
y = y+kc 
print y
f0 = f(0,y0)
f1 = f(0.1,y1)
f2 = f(0.2,y2)
f3 = f(0.3,y3)
print "y0 y1 y2 y3 are respectively: "
print y3,y2,y1,y0
print "f0 f1 f2 f3 are respectively: "
print f3,f2,f1,f0
print "finding y4 using predictors milne method x=0.4"
h = 0.1
y4 = y0+4*h/3*(2*f1-f2+2*f3)
print "y4 = ",y4
print "f4 = ",f(0.4,y4)
print "Using corrector method : "
y4 = y2+h/3*(f2+4*f3+f4) 
print "y4 = ",y4
print "f4 = ",f(0.4,y4) 

Runga kutta method
The required approximate value is :− 
1.0
To find y(0.4) put x=0.2 y=above value ie 1.196 h=0.2 
The required approximate value is :− 
1.0
To find y(0.4) put x=0.2 y=above value ie 1.196 h=0.2 
The required approximate value is :− 
1.0
y0 y1 y2 y3 are respectively: 
1.0 1.0 1.0 1
f0 f1 f2 f3 are respectively: 
1.3 1.2 1.1 1
finding y4 using predictors milne method x=0.4
y4 =  1.48
f4 =  2.7824
Using corrector method : 
y4 =  -0.0333333333333333*(-x**5/20 + x**2/2)**2 + 1.24
f4 =  -0.0133333333333333*(-x**5/20 + x**2/2)**2 + (-0.0333333333333333*(-x**5/20 + x**2/2)**2 + 1.24)**2 + 0.496

Example 27.16, page no. 705

In [32]:
def f(x,y):
    yy = x**2*(1+y)
    return yy
y3 = 1
y2 = 1.233
y1 = 1.548
y0 = 1.979
f3 = f(1,y3)
f2 = f(1.1,y2)
f1 = f(1.2,y1)
f0 = f(1.3,y0)
print "Using predictor method "
h = 0.1
y11 = y0+h/24*(55*f0-59*f1+37*f2-9*f3)
print "y11 = ",y11
x = 1.4
f11 = f(1.4,y11)
print "Using corrector method "
y11 = y0+h/24*(9*f11+19*f0-5*f1+f2)
print "y11 = ",y11
f11 = f(1.4,y11)
print "f11 = ",f11

Using predictor method 
y11 =  2.57229741667
Using corrector method 
y11 =  2.57494727679
f11 =  7.00689666251

Example 27.17, page no. 706

In [33]:
print "Runga kutta method"
def f(x,y):
    yy = x-y**2
    return yy
y0 = 1
x = 0
y = 1
h = 0.1
k1 = h*f(x,y)
k2 = h*f(x+1/2*h,y+1/2*k1)
k3 = h*f(x+1/2*h,y+1/2*k2)
k4 = h*f(x+h,y+k3)
ka = 1/6*(k1+2*k2+2*k3+k4) 
print "The required approximate value is :− "
y1 = y+ka
y = y+ka
print y
print "To find y(0.4) put x=0.2 y=above value ie 1.196 h=0.2 "
x = 0.1
h = 0.1
k1 = h*f(x,y)
k2 = h*f(x+1/2*h,y+1/2*k1)
k3 = h*f(x+1/2*h,y+1/2*k2)
k4 = h*f(x+h,y+k3)
kb = 1/6*(k1+2*k2+2*k3+k4)
print "The required approximate value is :− "
y2 = y+kb
y = y+kb
print y
print "To find y(0.4) put x=0.2 y=above value ie 1.196 h=0.2 "
x = 0.2
h = 0.1
k1 = h*f(x,y)
k2 = h*f(x+1/2*h,y+1/2*k1)
k3 = h*f(x+1/2*h,y+1/2*k2)
k4 = h*f(x+h,y+k3)
kc = 1/6*(k1+2*k2+2*k3+k4) 
print "The required approximate value is :− "
y3 = y+kc
y = y+kc
print y
f0 = f(0,y0)
f1 = f(0.1,y1)
f2 = f(0.2,y2)
f3 = f(0.3,y3) 
print "y0 y1 y2 y3 are respectively :"
print y3,y2,y1,y0
print "f0 f1 f2 f3 are respectively :"
print f3,f2,f1,f0
print "Using adams method "
print "Using the predictor"
h = 0.1
y4 = y3+h/24*(55*f3-59*f2+37*f1-9*f0)
x = 0.4
f4 = f(0.4,y4)
print "y4 = ",y4
print "Using corrector method "
y4 = y3+h/24*(9*f4+19*f3-5*f2+f1) 
print "y4 = ",y4
f4 = f(0.4,y4)
print "f4 = ",f4

Runga kutta method
The required approximate value is :− 
1.0
To find y(0.4) put x=0.2 y=above value ie 1.196 h=0.2 
The required approximate value is :− 
1.0
To find y(0.4) put x=0.2 y=above value ie 1.196 h=0.2 
The required approximate value is :− 
1.0
y0 y1 y2 y3 are respectively :
1.0 1.0 1.0 1
f0 f1 f2 f3 are respectively :
-0.7 -0.8 -0.9 -1
Using adams method 
Using the predictor
y4 =  0.935
Using corrector method 
y4 =  0.9397165625
f4 =  -0.483067217837

Example 27.18, page no. 708

In [36]:
import sympy

print "Picards method"
x0 = 0
y0 = 2
z0 = 1
x = sympy.Symbol('x')
def f(x,y,z):
    yy = x+z
    return yy
def g(x,y,z):
    yy = x-y**2
    return yy
print "First approximation"
y1 = y0+sympy.integrate(f(x,y0,z0),(x,x0,x))
print "y1 = ",y1
z1 = z0+sympy.integrate(g(x,y0,z0),(x,x0,x))
print "z1 = ",z1
print "Second approximation"
y2 = y0+sympy.integrate(f(x,y1,z1),(x,x0,x))
print "y2 = ",y2
z2 = z0+sympy.integrate(g(x,y1,z1),(x,x0,x))
print "z2 = ",z2
print "Third approximation"
y3 = y0+sympy.integrate(f(x,y2,z2),(x,x0,x))
print "y3 = ",y3
z3 = z0+sympy.integrate(g(x,y2,z2),(x,x0,x))
print "z3 = ",z3
x = 0.1
print "y(0.1) = ",y3.evalf()
print "z(0.1) = ",z3.evalf()

 Picards method
First approximation
y1 =  x**2/2 + x + 2
z1 =  x**2/2 - 4*x + 1
Second approximation
y2 =  x**3/6 - 3*x**2/2 + x + 2
z2 =  -x**5/20 - x**4/4 - x**3 - 3*x**2/2 - 4*x + 1
Third approximation
y3 =  -x**6/120 - x**5/20 - x**4/4 - x**3/2 - 3*x**2/2 + x + 2
z3 =  -x**7/252 + x**6/12 - 31*x**5/60 + 7*x**4/12 + 5*x**3/3 - 3*x**2/2 - 4*x + 1
y(0.1) =  -0.00833333333333333*x**6 - 0.05*x**5 - 0.25*x**4 - 0.5*x**3 - 1.5*x**2 + x + 2.0
z(0.1) =  -0.00396825396825397*x**7 + 0.0833333333333333*x**6 - 0.516666666666667*x**5 + 0.583333333333333*x**4 + 1.66666666666667*x**3 - 1.5*x**2 - 4.0*x + 1.0

Example 27.19, page no. 710

In [37]:
import sympy

x = sympy.Symbol('x')
def f(x,y,z):
    yy = z 
    return yy
def g(x,y,z):
    yy = x*y**2-y**2
    return yy
x0 = 0
y0 = 1
z0 = 0
h = 0.2
print "Using k1 k2.. for f and l1 l2... for g runga kutta formula becomes "
h = 0.2
k1 = h*f(x0,y0,z0)
l1 = h*g(x0,y0,z0)
k2 = h*f(x0+1/2*h,y0+1/2*k1,z0+1/2*l1)
l2 = h*g(x0+1/2*h,y0+1/2*k1,z0+1/2*l1)
k3 = h*f(x0+1/2*h,y0+1/2*k2,z0+1/2*l2)
l3 = h*g(x0+1/2*h,y0+1/2*k2,z0+1/2*l2)
k4 = h*f(x0+h,y0+k3,z0+l3)
l4 = h*g(x0+h,y0+k3,z0+l3)
k = 1/6*(k1+2*k2+2*k3+k4)
l = 1/6*(l1+2*l2+2*l3+2*l4)
x = 0.2
y = y0+k
y1 = z0+l
print "y = ",y
print "y1 = ",y1

Using k1 k2.. for f and l1 l2... for g runga kutta formula becomes 
y =  1.0
y1 =  0.0