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