Chapter08:Numerical Solution of Ordinary Differential Equations

Ex8.1:pg-304

In [2]:
#example 8.1
#taylor's method
#page 304
import math
f=1    #value of function at 0
def f1(x):
    return x-f**2
def f2(x):
    return 1-2*f*f1(x)
def f3(x):
    return -2*f*f2(x)-2*f2(x)**2
def f4(x):
    return -2*f*f3(x)-6*f1(x)*f2(x)
def f5(x):
    return -2*f*f4(x)-8*f1(x)*f3(x)-6*f2(x)**2
h=0.1   #value at 0.1
k=f 
for j in range(1,5):
    if j==1:
        k=k+h*f1(0);
    elif j==2:
        k=k+(h**j)*f2(0)/math.factorial(j)
    elif j ==3:
        k=k+(h**j)*f3(0)/math.factorial(j)
    elif j ==4:
        k=k+(h**j)*f4(0)/math.factorial(j)
    elif j==5:
        k=k+(h**j)*f5(0)/math.factorial(j)
print "the value of the function at %.2f is :%0.4f" %(h,k)

the value of the function at 0.10 is :0.9113

Ex8.2:pg-304

In [4]:
#taylor's method
#example 8.2
#page 304
import math
f=1   #value of function at 0
f1=0   #value of first derivatie at 0
def f2(x):
    return x*f1+f
def f3(x):
    return x*f2(x)+2*f1
def f4(x):
    return x*f3(x)+3*f2(x)
def f5(x):
    return x*f4(x)+4*f3(x)
def f6(x):
    return x*f5(x)+5*f4(x)
h=0.1   #value at 0.1
k=f
for j in range(1,6):
    if j==1:
        k=k+h*f1
    elif j==2:
        k=k+(h**j)*f2(0)/math.factorial(j)
    elif j ==3:
        k=k+(h**j)*f3(0)/math.factorial(j)
    elif j ==4:
        k=k+(h**j)*f4(0)/math.factorial(j)
    elif j==5:
        k=k+(h**j)*f5(0)/math.factorial(j)
    else:
        k=k+(h**j)*f6(0)/math.factorial (j)
print "the value of the function at %.2f is :%0.7f" %(h,k)

the value of the function at 0.10 is :1.0050125

Ex8.3:pg-306

In [22]:
#example 8.3
#picard's method
#page 306
from scipy import integrate
from __future__ import division
def f(x,y):
    return x+y**2
y=[0,0,0,0]
y[1]=1
for i in range(1,3):
    a=integrate.quad(lambda x:x+y[i]**2,0,i/10)
    y[i+1]=a[0]+y[1]
    print "\n y (%g) = %g\n" %(i/10,y[i+1])

 y (0.1) = 1.105


 y (0.2) = 1.26421

Ex8.4:pg-306

In [16]:
#example 8.4
#picard's method
#page 306
from scipy import integrate
y=[0,0,0,0]   #value at 0
c=0.25
for i in range(0,3):
    a=integrate.quad(lambda x:(x**2/(y[i]**2+1)),0,c)
    y[i+1]=y[0]+a[0]
    print "\n y(%0.2f) = %g\n" %(c,y[i+1])
    c=c*2

 y(0.25) = 0.00520833


 y(0.50) = 0.0416655


 y(1.00) = 0.332756

Ex8.5:pg-308

In [6]:
#example 8.5
#euler's method
#page 308
def f(y):
    return -1*y
y=[0,0,0,0,0]
y[0]=1   #value at 0
h=0.01
c=0.01
for i in range(0,4):
    y[i+1]=y[i]+h*f(y[i])
    print "\ny(%g)=%g\n" %(c,y[i+1])
    c=c+0.01

y(0.01)=0.99


y(0.02)=0.9801


y(0.03)=0.970299


y(0.04)=0.960596

Ex8.6:pg-308

In [10]:
#example 8.6
#error estimates in euler's 
#page 308
import math
from __future__ import division
def f(y):
    return -1*y
y=[0,0,0,0,0]
L=[0,0,0,0,0]
e=[0,0,0,0,0]
y[0]=1    #value at 0
h=0.01
c=0.01;
for i in range(0,4):
    y[i+1]=y[i]+h*f(y[i])
    print "\ny(%g)=%g\n" %(c,y[i+1])
    c=c+0.01
for i in range(0,4):
    L[i]=abs(-(1/2)*(h**2)*y[i+1])
    print "L(%d) =%f\n\n" %(i,L[i])
e[0]=0
for i in range(0,4):
    e[i+1]=abs(y[1]*e[i]+L[0])
    print "e(%d)=%f\n\n" %(i,e[i])
Actual_value=math.exp(-0.04)
Estimated_value=y[4]
err=abs(Actual_value-Estimated_value)
if err<e[4]:
    print "VERIFIED"

y(0.01)=0.99


y(0.02)=0.9801


y(0.03)=0.970299


y(0.04)=0.960596

L(0) =0.000050


L(1) =0.000049


L(2) =0.000049


L(3) =0.000048


e(0)=0.000000


e(1)=0.000050


e(2)=0.000099


e(3)=0.000147


VERIFIED

Ex8.7:pg-310

In [25]:
#example 8.7
#modified euler's method
#page 310
h=0.05
f=1
def f1(x,y):
    return x**2+y
x=[0,0.05,0.1]
y1=[0,0,0,0]
y2=[0,0,0,0]
y1[0]=f+h*f1(x[0],f);
y1[1]=f+h*(f1(x[0],f)+f1(x[1],y1[0]))/2
y1[2]=f+h*(f1(x[0],f)+f1(x[2],y1[1]))/2
y2[0]=y1[1]+h*f1(x[1],y1[1])
y2[1]=y1[1]+h*(f1(x[1],y1[1])+f1(x[2],y2[0]))/2
y2[2]=y1[1]+h*(f1(x[1],y1[1])+f1(x[2],y2[1]))/2
print "y1(0)\t         y1(1)\t      y1(2)\t       y2(0)\t        y2(1)\t        y3(2)\n\n"
print " %f\t    %f\t    %f\t    %f\t     %f\t    %f\n" %(y1[0],y1[1],y1[2],y2[0],y2[1],y2[2])
print "\n\n the value of y at %0.2f is : %0.4f" %(x[2],y2[2])

y1(0)	         y1(1)	      y1(2)	       y2(0)	        y2(1)	        y3(2)


 1.050000	    1.051313	    1.051533	    1.104003	     1.105508	    1.105546



 the value of y at 0.10 is : 1.1055

Ex8.8:pg-313

In [29]:
#example 8.8
#runge-kutta formula
#page 313
from __future__ import division
def f(x,y):
    return y-x
y=2
x=0
h=0.1
K1=h*f(x,y)
K2=h*f(x+h,y+K1)
y1=y+( K1+K2)/2
print "\n y(0.1) by second order runge kutta method:%0.4f" %(y1)
y=y1
x=0.1
h=0.1
K1=h*f(x,y)
K2=h*f(x+h,y+K1)
y1=y+( K1+K2)/2
print "\n y(0.2) by second order runge kutta method:%0.4f" %(y1)
y=2
x=0
h=0.1
K1=h*f(x,y)
K2=h*f(x+h/2,y+K1/2)
K3=h*f(x+h/2,y+K2/2)
K4=h*f(x+h,y+K3)
y1=y+(K1+2*K2+2*K3+K4)/6
print "\n y(0.1) by fourth order runge kutta method:%0.4f" %(y1)
y=y1
x=0.1
h=0.1
K1=h*f(x,y)
K2=h*f(x+h/2,y+K1/2)
K3=h*f(x+h/2,y+K2/2)
K4=h*f(x+h,y+K3)
y1=y+(K1+2*K2+2*K3+K4)/6
print "\n y(0.1) by fourth order runge kutta method:%0.4f " %(y1)

 y(0.1) by second order runge kutta method:2.2050

 y(0.2) by second order runge kutta method:2.4210

 y(0.1) by fourth order runge kutta method:2.2052

 y(0.1) by fourth order runge kutta method:2.4214

Ex8.9:pg-315

In [30]:
#example 8.9
#runge kutta method
#page  315
from __future__ import division
def f(x,y):
    return 1+y**2
y=0
x=0
h=0.2
K1=h*f(x,y)
K2=h*f(x+h/2,y+K1/2)
K3=h*f(x+h/2,y+K2/2)
K4=h*f(x+h,y+K3)
y1=y+(K1+2*K2+2*K3+K4)/6
print "\n y(0.2) by fourth order runge kutta method:%0.4f" %(y1)
y=y1
x=0.2
h=0.2
K1=h*f(x,y)
K2=h*f(x+h/2,y+K1/2)
K3=h*f(x+h/2,y+K2/2)
K4=h*f(x+h,y+K3)
y1=y+(K1+2*K2+2*K3+K4)/6
print "\n y(0.4) by fourth order runge kutta method:%0.4f" %(y1)
y=2
x=0
h=0.1
y=y1
x=0.4
h=0.2
K1=h*f(x,y)
K2=h*f(x+h/2,y+K1/2)
K3=h*f(x+h/2,y+K2/2)
K4=h*f(x+h,y+K3)
y1=y+(K1+2*K2+2*K3+K4)/6
print "\n y(0.6) by fourth order runge kutta method:%0.4f" %(y1)

 y(0.2) by fourth order runge kutta method:0.2027

 y(0.4) by fourth order runge kutta method:0.4228

 y(0.6) by fourth order runge kutta method:0.6841

Ex8.10:pg-315

In [35]:
#example 8.10
#initial value problems
#page 315
from __future__ import division
def f1(x,y):
    return 3*x+y/2
y=[1,0,0]
h=0.1
c=0
for i in range(0,2):
    y[i+1]=y[i]+h*f1(c,y[i])
    print "\ny(%g)=%g\n" %(c,y[i])
    c=c+0.1

y(0)=1


y(0.1)=1.05

Ex8.11:pg-316

In [1]:
#example 8.11
#adam's moulton method
#page  316
def f(x,y):
    return 1+y**2
y=0
x=0
h=0.2
f1=[0,0,0]
K1=h*f(x,y)
K2=h*f(x+h/2,y+K1/2)
K3=h*f(x+h/2,y+K2/2)
K4=h*f(x+h,y+K3)
y1=y+(K1+2*K2+2*K3+K4)/6
f1[0]=y1
print "\n y(0.2) by fourth order runge kutta method:%0.4f" %(y1)
y=y1
x=0.2
h=0.2
K1=h*f(x,y)
K2=h*f(x+h/2,y+K1/2)
K3=h*f(x+h/2,y+K2/2)
K4=h*f(x+h,y+K3)
y1=y+(K1+2*K2+2*K3+K4)/6
f1[1]=y1
print "\n y(0.4) by fourth order runge kutta method:%0.4f" %(y1)
y=2
x=0
h=0.1
y=y1
x=0.4
h=0.2
K1=h*f(x,y)
K2=h*f(x+h/2,y+K1/2)
K3=h*f(x+h/2,y+K2/2)
K4=h*f(x+h,y+K3)
y1=y+(K1+2*K2+2*K3+K4)/6
f1[2]=y1
print "\n y(0.6) by fourth order runge kutta method:%0.4f" %(y1)
y_p=y1+h*(55*(1+f1[2]**2)-59*(1+f1[1]**2)+37*(1+f1[0]**2)-9)/24
y_c=y1+h*(9*(1+(y_p-1)**2)+19*(1+f1[2]**2)-5*(1+f1[1]**2)+(1+f1[0]**2))/24
print "\nthe predicted value is:%0.4f:\n" %(y_p)
print " the computed value is:%0.4f:" %(y_c)

 y(0.2) by fourth order runge kutta method:0.2027

 y(0.4) by fourth order runge kutta method:0.4228

 y(0.6) by fourth order runge kutta method:0.6841

the predicted value is:1.0234:

 the computed value is:0.9512:

Ex8.12:pg-320

In [4]:
#example 8.12
#milne's method
#page 320
def f(x,y):
    return 1+y**2
y=0
f1=[0,0,0]
Y1=[0,0,0,0]
x=0
h=0.2
f1[0]=0
print "x                y              y1=1+y^2\n\n"
Y1[0]=1+y**2
print "%0.4f         %0.4f             %0.4f\n" %(x,y,(1+y**2))
K1=h*f(x,y)
K2=h*f(x+h/2,y+K1/2)
K3=h*f(x+h/2,y+K2/2)
K4=h*f(x+h,y+K3)
y1=y+(K1+2*K2+2*K3+K4)/6
f1[0]=y1
Y1[1]=1+y1**2
print "%0.4f         %0.4f             %0.4f\n" %(x+h,y1,(1+y1**2))
y=y1
x=0.2
h=0.2
K1=h*f(x,y)
K2=h*f(x+h/2,y+K1/2)
K3=h*f(x+h/2,y+K2/2)
K4=h*f(x+h,y+K3)
y1=y+(K1+2*K2+2*K3+K4)/6
f1[1]=y1
Y1[2]=1+y1**2
print "%0.4f         %0.4f             %0.4f\n" %(x+h,y1,(1+y1**2))
y=y1
x=0.4
h=0.2
K1=h*f(x,y)
K2=h*f(x+h/2,y+K1/2)
K3=h*f(x+h/2,y+K2/2)
K4=h*f(x+h,y+K3)
y1=y+(K1+2*K2+2*K3+K4)/6
f1[2]=y1
Y1[3]=1+y1**2;
print "%0.4f         %0.4f             %0.4f\n" %(x+h,y1,(1+y1**2))
Y_4=4*h*(2*Y1[1]-Y1[2]+2*Y1[3])/3
print "y(0.8)=%f\n" %(Y_4)
Y=1+Y_4**2
Y_4=f1[1]+h*(Y1[2]+4*Y1[3]+Y)/3    #more correct value
print "y(0.8)=%f\n" %(Y_4)

x                y              y1=1+y^2


0.0000         0.0000             1.0000

0.2000         0.2027             1.0411

0.4000         0.4228             1.1788

0.6000         0.6841             1.4680

y(0.8)=1.023869

y(0.8)=1.029403

Ex8.13:pg-320

In [8]:
#example 8.13
#milne's method
#page 320
def f1(x,y):
    return x**2+y**2-2
x=[-0.1, 0, 0.1, 0.2]
y=[1.0900, 1.0, 0.8900, 0.7605]
Y1=[0,0,0,0]
h=0.1
for i in range(0,4):
    Y1[i]=f1(x[i],y[i])
print " x                y                y1=x^2+y^2-2    \n\n"
for i in range(0,4):
    print " %0.2f           %f           %f           \n" %(x[i],y[i],Y1[i])
Y_3=y[0]+(4*h/3)*(2*Y1[1]-Y1[2]+2*Y1[3])
print "y(0.3)=%f\n" %(Y_3)
Y1_3=f1(0.3,Y_3)
Y_3=y[2]+h*(Y1[2]+4*Y1[3]+Y1_3)/3   #corrected value
print "corrected y(0.3)=%f" %(Y_3)

 x                y                y1=x^2+y^2-2    


 -0.10           1.090000           -0.801900           

 0.00           1.000000           -1.000000           

 0.10           0.890000           -1.197900           

 0.20           0.760500           -1.381640           

y(0.3)=0.614616

corrected y(0.3)=0.614776

Ex8.14:pg322

In [10]:
#example 8.14
#initial-value problem
#page 322
from __future__ import division
import math
def f(x):
    return 13*math.exp(x/2)-6*x-12
s1=1.691358
s3=3.430879
print "the erorr in the computed values are %0.7g  %0.7g" %(abs(f(0.5)-s1),abs(f(1)-s3))

the erorr in the computed values are 0.0009724169  0.002497519

Ex8.15:pg-328

In [13]:
#boundary value problem using finite difference method
#example 8.15
#page 328
import math
from numpy import matrix
def f(x):
    return math.cos(x)+((1-math.cos(1))/math.sin(1))*math.sin(x)-1
h1=1/2
Y=f(0.5)
y0=0
y2=0
y1=4*(1/4+y0+y2)/7
print "computed value with h=%f of y(0.5) is %f\n" %(h1,y1)
print "error in the result with actual value %f\n" %(abs(Y-y1))
h2=1/4
y0=0
y4=0
#solving the approximated diffrential equation
A=matrix([[-31/16, 1, 0],[1, -31/16, 1],[0, 1, -31/16]])
X=matrix([[-1/16],[-1/16],[-1/16]])
C=A.I*X
print "computed value with h=%f of y(0.5) is %f\n" %(h2,C[1][0])
print "error in the result with actual value %f\n" %(abs(Y-C[1][0]))

computed value with h=0.500000 of y(0.5) is 0.142857

error in the result with actual value 0.003363

computed value with h=0.250000 of y(0.5) is 0.140312

error in the result with actual value 0.000818

Ex8.16:pg-329

In [14]:
#boundary value problem using finite difference method
#example 8.16
#page 329\
from numpy import matrix
import math
def f(x):
    return math.sinh(x)
y0=0    #y(0)=0
y4=3.62686         #y(2)=3.62686
h1=0.5
Y=f(0.5)
#arranging and calculating the values
A=matrix([[-9, 4, 0],[4, -9, 4],[0, 4, -9]])
C=matrix([[0],[0],[-14.50744]])
X=A.I*C
print "computed value with h=%f of y(0.5) is %f\n" %(h1,X[0][0])
print "error in the result with actual value %f\n" %(abs(Y-X[0][0]))
h2=1.0
y0=0    #y(0)=0
y2=3.62686   #y(2)=3.62686
y1=(y0+y2)/3
Y=(4*X[1][0]-y1)/3
print "with better approximation error is reduced to %f" %(abs(Y-f(1.0)))

computed value with h=0.500000 of y(0.5) is 0.526347

error in the result with actual value 0.005252

with better approximation error is reduced to 0.000855

Ex8.17:pg-330

In [15]:
#cubic spline method
#example 8.17
#page 330
import math
def f(x):
    return math.cos(x)+((1-math.cos(1))/math.sin(1))*math.sin(x)-1
y=[f(0), f(0.5), f(1)]
h=1/2
Y=f(0.5)
y0=0
y2=0
M0=-1
M1=-22/25
M2=-1
s1_0=47/88
s1_1=-47/88
s1_05=0
print "the cubic spline value is %f" %(Y)

the cubic spline value is 0.139494

Ex8.18:pg-331

In [17]:
#cubic spline method
#example 8.18
#page 331
from numpy import matrix
from __future__ import division
import math
#after arranging and forming equation 
A=matrix([[10, -1, 0, 24],[0, 16, -1, -32],[1, 20, 0, 16],[0, 1, 26, -24]])
C=matrix([[36],[-12],[24],[-9]])
X=A.I*C;
print " Y1=%f\n\n" %(X[3][0])
print " the error in the solution is :%f" %(abs((2/3)-X[3][0]))

 Y1=0.653890


 the error in the solution is :0.012777

Ex8.19:pg-331

In [18]:
#boundary value problem by cubisc spline nethod
#example 8.18
#page 331
from numpy import matrix
h=1/2
#arranging in two subintervals we get
A=matrix([[10, -1, 0,24],[0, 16, -1, -32],[1, 20, 0, 16],[0, 1, 26, -24]])
C=matrix([[36],[-12],[24],[-9]])
X=A.I*C
print "the computed value of y(1.5) is %f "%(X[3][0])

the computed value of y(1.5) is 0.653890