# Chapter09:Numerical Solution of Partial Differential Equations¶

## Ex9.1:pg-350¶

In [30]:
#standard five point formula
#example 9.1
#page 350

u2=5.0
u3=1.0
for i in range(0,3):
u1=(u2+u3+6.0)/4.0
u2=(u1/2.0)+(5.0/2.0)
u3=(u1/2.0)+(1.0/2.0)
print" the values are u1=%d\t u2=%d\t u3=%d\t\n\n" %(u1,u2,u3)


 the values are u1=3	 u2=4	 u3=2

the values are u1=3	 u2=4	 u3=2

the values are u1=3	 u2=4	 u3=2



## Ex9.2:pg-351¶

In [51]:
#solution of laplace equation by jacobi method,gauss-seidel method and SOR method
#example 9.2
#page 351
u1=0.25
u2=0.25
u3=0.5
u4=0.5 #initial values
print "jacobis iteration process\n\n"
print"u1\t     u2\t     u3\t     u4\t \n\n"
print "%f\t     %f\t     %f\t     %f\t  \n" %(u1,u2,u3,u4)
for i in range(0,7):
u11=(0+u2+0+u4)/4
u22=(u1+0+0+u3)/4
u33=(1+u2+0+u4)/4
u44=(1+0+u3+u1)/4
u1=u11
u2=u22
u3=u33
u4=u44
print "%f\t     %f\t     %f\t     %f\t  \n" %(u11,u22,u33,u44)
print " gauss seidel process\n\n"
u1=0.25
u2=0.3125
u3=0.5625
u4=0.46875  #initial values
print "u1\t     u2\t     u3\t     u4\t \n\n"
print "%f\t     %f\t     %f\t     %f\t  \n" %(u1,u2,u3,u4)
for i in range(0,4):

u1=(0.0+u2+0.0+u4)/4.0
u2=(u1+0.0+0.0+u3)/4.0
u3=(1.0+u2+0.0+u4)/4.0
u4=(1.0+0.0+u3+u1)/4.0
print "%f\t     %f\t     %f\t     %f\t  \n" %(u1,u2,u3,u4)

jacobis iteration process

u1	     u2	     u3	     u4

0.250000	     0.250000	     0.500000	     0.500000

0.187500	     0.187500	     0.437500	     0.437500

0.156250	     0.156250	     0.406250	     0.406250

0.140625	     0.140625	     0.390625	     0.390625

0.132812	     0.132812	     0.382812	     0.382812

0.128906	     0.128906	     0.378906	     0.378906

0.126953	     0.126953	     0.376953	     0.376953

0.125977	     0.125977	     0.375977	     0.375977

gauss seidel process

u1	     u2	     u3	     u4

0.250000	     0.312500	     0.562500	     0.468750

0.195312	     0.189453	     0.414551	     0.402466

0.147980	     0.140633	     0.385775	     0.383439

0.131018	     0.129198	     0.378159	     0.377294

0.126623	     0.126196	     0.375872	     0.375624



## Ex9.4:pg-354¶

In [59]:
#poisson equation
#exaample 9.4
#page 354
u2=0.0
u4=0.0
print " u1\t    u2\t   u3\t   u4\t\n\n"
for i in range(0,6):
u1=(u2/2.0)+30.0
u2=(u1+u4+150.0)/4.0
u4=(u2/2.0)+45.0
print "%0.2f\t   %0.2f\t   %0.2f\t   %0.2f\n" %(u1,u2,u2,u4)
print " from last two iterates we conclude u1=67   u2=75   u3=75   u4=83\n"

 u1	    u2	   u3	   u4

30.00	   45.00	   45.00	   67.50

52.50	   67.50	   67.50	   78.75

63.75	   73.12	   73.12	   81.56

66.56	   74.53	   74.53	   82.27

67.27	   74.88	   74.88	   82.44

67.44	   74.97	   74.97	   82.49

from last two iterates we conclude u1=67   u2=75   u3=75   u4=83



## Ex9.6:pg-362¶

In [77]:
#bender-schmidt formula
#example 9.6
#page 362
def f(x):
return (4*x)-(x*x)
#u=[f(0),f(1),f(2),f(3),f(4)]
u1=f(0);u2=f(1);u3=f(2);u4=f(3);u5=f(4);
u11=(u1+u3)/2
u12=(u2+u4)/2
u13=(u3+u5)/2
print "u11=%0.2f\t  u12=%0.2f\t  u13=%0.2f\t \n" %(u11,u12,u13)
u21=(u1+u12)/2.0
u22=(u11+u13)/2.0
u23=(u12+0)/2.0
print "u21=%0.2f\t  u22=%0.2f\t  u23=%0.2f\t \n" %(u21,u22,u23)
u31=(u1+u22)/2.0
u32=(u21+u23)/2.0
u33=(u22+u1)/2.0
print "u31=%0.2f\t  u32=%0.2f\t  u33=%0.2f\t \n" % (u31,u32,u33)
u41=(u1+u32)/2.0
u42=(u31+u33)/2.0
u43=(u32+u1)/2.0
print "u41=%0.2f\t  u42=%0.2f\t  u43=%0.2f\t \n" % (u41,u42,u43)
u51=(u1+u42)/2.0
u52=(u41+u43)/2.0
u53=(u42+u1)/2.0
print "u51=%0.2f\t  u52=%0.2f\t  u53=%0.2f\t \n" % (u51,u52,u53)

u11=2.00	  u12=3.00	  u13=2.00

u21=1.50	  u22=2.00	  u23=1.50

u31=1.00	  u32=1.50	  u33=1.00

u41=0.75	  u42=1.00	  u43=0.75

u51=0.50	  u52=0.75	  u53=0.50



## Ex9.7:pg-363¶

In [28]:
#bender-schimdt's formula and crank-nicolson formula
#example 9.7
#page 363
#bender -schimdt's formula
import math
from numpy import matrix
z=math.pi
def f(x,t):
return math.exp(z*z*t*-1)*sin(z*x)
#u=[f(0,0),f(0.2,0),f(0.4,0),f(0.6,0),f(0.8,0),f(1,0)]
u1=f(0,0)
u2=f(0.2,0)
u3=f(0.4,0)
u4=f(0.6,0)
u5=f(0.8,0)
u6=f(1.0,0)
u11=u3/2
u12=(u2+u4)/2
u13=u12
u14=u11
print "u11=%f\t  u12=%f\t  u13=%f\t  u14=%f\n\n" % (u11,u12,u13,u14)
u21=u12/2
u22=(u12+u14)/2
u23=u22
u24=u21
print "u21=%f\t  u22=%f\t  u23=%f\t  u24=%f\n\n" % (u21,u22,u23,u24)
print "the error in the solution is: %f\n\n" % (math.fabs(u22-f(0.6,0.04)))
#crank-nicolson formula
#by putting i=1,2,3,4 we obtain four equation
A=matrix([[4, -1, 0, 0] ,[-1, 4, -1, 0],[0, -1, 4, -1],[0, 0, -1, 4]])
C=matrix([[0.9510],[1.5388],[1.5388],[0.9510]])
X=A.I*C
print "u00=%f\t  u10=%f\t  u20=%f\t  u30=%f\t\n\n" %(X[0][0],X[1][0],X[2][0],X[3][0])
print "the error in the solution is: %f\n\n" %(abs(X[1][0]-f(0.6,0.04)))

u11=0.475528	  u12=0.769421	  u13=0.769421	  u14=0.475528

u21=0.384710	  u22=0.622475	  u23=0.622475	  u24=0.384710

the error in the solution is: 0.018372

u00=0.399255	  u10=0.646018	  u20=0.646018	  u30=0.399255

the error in the solution is: 0.005172



## Ex9.8:pg-364¶

In [1]:
#heat equation using crank-nicolson method
#example 9.8
#page 364
from numpy import matrix
import math
z=0.01878
#h=1/2 l=1/8,i=1
u01=0.0
u21=1.0/8.0
u11=(u21+u01)/6.0
print " u11=%f\n\n" % (u11)
print "error is %f\n\n" % (math.fabs(u11-z))
#h=1/4,l=1/8,i=1,2,3
A=matrix([[-3.0 ,-1.0 ,0.0],[1.0,-3.0,1.0],[0.0,1.0,-3.0]])
C=matrix([[0.0],[0.0],[-0.125]])
#here we found inverese of A then we multipy it with C
X=A.I*C
print "u12=%f\n\n" % (X[1][0])
print "error is %f\n\n" %(math.fabs(X[1][0]-z))

 u11=0.020833

error is 0.002053

u12=0.013889

error is 0.004891