Chapter 28: Numerical Solution Of Partial Differential Equations

Example 28.1, page no. 725

In [1]:
print "D=Bˆ2−4AC"
print "if D<0 then elliptic if D=0 then parabolic if D>0 then hyperboic"
print "(i) A=xˆ2, B1−yˆ2 D=4ˆ2−4∗1∗4=0 so the equation is PARABOLIC "
print "(ii) D=4xˆ2(yˆ2−1)"
print "for −inf<x<inf and −1<y<1 D<0 "
print "So the equation is ELLIPTIC "
print "(iii) A=1+xˆ2, B=5+2xˆ2, C=4+xˆ2"
print "D=9>0"
print "So the equation is HYPERBOLIC "

D=Bˆ2−4AC
if D<0 then elliptic if D=0 then parabolic if D>0 then hyperboic
(i) A=xˆ2, B1−yˆ2 D=4ˆ2−4∗1∗4=0 so the equation is PARABOLIC 
(ii) D=4xˆ2(yˆ2−1)
for −inf<x<inf and −1<y<1 D<0 
So the equation is ELLIPTIC 
(iii) A=1+xˆ2, B=5+2xˆ2, C=4+xˆ2
D=9>0
So the equation is HYPERBOLIC

Example 28.2, page no. 726

In [5]:
print "See figure in question"
print "From symmetry u7=u1, u8=u2, u9=u3, u3=u1, u6=u4, u9=u7"
print "u5=1/4∗(2000+2000+1000+1000)=1500"
u5 = 1500
print "u1=1/4(0=1500+1000+2000)=1125"
u1 = 1125
print "u2=1/4∗(1125+1125+1000+1500)=1188"
u2 = 1188
print "u4=1/4(2000+1500+1125+1125)=1438"
u4 = 1438
print u1,u2,u4,u5 
print "Iterations :"
for i in range(1,7):
    u11 = (1000+u2+500+u4)/4
    u22 = (u11+u1+1000+u5)/4 
    u44 = (2000+u5+u11+u1)/4
    u55 = (u44+u4+u22+u2)/4
    print ""
    print u55,u44,u22,u11
    u1 = u11 
    u2 = u22 
    u4 = u44 
    u5 = u55

See figure in question
From symmetry u7=u1, u8=u2, u9=u3, u3=u1, u6=u4, u9=u7
u5=1/4∗(2000+2000+1000+1000)=1500
u1=1/4(0=1500+1000+2000)=1125
u2=1/4∗(1125+1125+1000+1500)=1188
u4=1/4(2000+1500+1125+1125)=1438
1125 1188 1438 1500
Iterations :

1301 1414 1164 1031

1250 1337 1087 1019

1199 1312 1062 981

1174 1287 1037 968

1155 1274 1024 956

1144 1265 1015 949

Example 28.3, page no. 727

In [6]:
print "See figure in question"
print "To find the initial values of u1 u2 u3 u4 we assume u4=0 "
print "u1=1/4∗(1000+0+1000+2000)=1000"
u1 = 1000
print "u2=1/4(1000+500+1000+500)=625 "
u2 = 625
print "u3=1/4∗(2000+0+1000+500)=875"
u3 = 875
print "u4=1/4(875+0+625+0)=375"
u4 = 375
print u1,u2,u3,u4 
print "Iterations:"
for i in range(1,7):
    u11 = (2000+u2+1000+u3)/4 
    u22 = (u11+500+1000+u4)/4
    u33 = (2000+u4+u11+500)/4
    u44 = (u33+0+u22+0)/4
    print ""
    print u44,u33,u22,u11
    u1 = u11 
    u2 = u22 
    u4 = u44 
    u3 = u33

See figure in question
To find the initial values of u1 u2 u3 u4 we assume u4=0 
u1=1/4∗(1000+0+1000+2000)=1000
u2=1/4(1000+500+1000+500)=625 
u3=1/4∗(2000+0+1000+500)=875
u4=1/4(875+0+625+0)=375
1000 625 875 375
Iterations:

437 1000 750 1125

453 1031 781 1187

457 1039 789 1203

458 1041 791 1207

458 1041 791 1208

458 1041 791 1208

Example 28.5, page no. 729

In [10]:
import numpy
print "Here cˆ2=4, h=1, k=1/8, therefore alpha=(cˆ2)∗k/(hˆ2)"
print "Using bendre−schmidits recurrence relation ie u(i)(j +1)=t∗u(i−1)(j)+t∗u(i+1)(j)+(1−2t)∗u(i,j)"
print "Now since u(0,t)=0=u(8,t) therefore u(0,i)=0 and u(8,j)=0 and u(x,0)=4x−1/2xˆ2 "
c = 2
h = 1
k = 1/8
t = (c**2)*k/(h**2)
A = numpy.ones((9,9))
for i in range(0,9):
    for j in range (0,9):
        A[0,i] = 0
        A[8,i] = 0
        A[i,0] = 4*(i-1)-1/2*(i-1)**2
for i in range(1,8):
    for j in range(1,7):
        A[i,j] = t*A[i-1,j-1]+t*A[i+1,j-1]+(1-2*t)*A[i-1,j-1]
for i in range(1,8):
    j = 2
    print A[i,j]

Here cˆ2=4, h=1, k=1/8, therefore alpha=(cˆ2)∗k/(hˆ2)
Using bendre−schmidits recurrence relation ie u(i)(j +1)=t∗u(i−1)(j)+t∗u(i+1)(j)+(1−2t)∗u(i,j)
Now since u(0,t)=0=u(8,t) therefore u(0,i)=0 and u(8,j)=0 and u(x,0)=4x−1/2xˆ2 
0.0
-4.0
0.0
4.0
8.0
12.0
16.0

Example 28.6, page no. 730

In [3]:
import numpy,math

print "Here cˆ2=1, h=1/3, k=1/36, therefore t=(cˆ2)∗k/(hˆ2)=1/4 "
print "So bendre−schmidits recurrence relation ie u(i)(j+1)=1/4(u(i−1)(j)+u(i+1)(j)+2u(i,j)"
print "Now since u(0,t)=0=u(1,t) therefore u(0,i)=0 and u(1,j)=0 and u(x,0)=sin(%pi)x"
c = 1.
h = 1./3
k = 1./36
t = (c**2)*k/(h**2)
A = numpy.ones((9,9))
for i in range(0,9):
    for j in range(0,9):
        A[0,i] = 0
        A[1,i] = 0
        A[i,0] = math.sin(math.pi/3*(i-1)) 
for i in range(1,8):
    for j in range(1,8):
        A[i,j] = t*A[i-1,j-1]+t*A[i+1,j-1]+(1-2*t)*A[i-1,j-1]
for i in range(1,8):
    j = 1
    print A[i,j]

Here cˆ2=1, h=1/3, k=1/36, therefore t=(cˆ2)∗k/(hˆ2)=1/4 
So bendre−schmidits recurrence relation ie u(i)(j+1)=1/4(u(i−1)(j)+u(i+1)(j)+2u(i,j)
Now since u(0,t)=0=u(1,t) therefore u(0,i)=0 and u(1,j)=0 and u(x,0)=sin(%pi)x
-0.433012701892
0.216506350946
0.649519052838
0.433012701892
-0.216506350946
-0.649519052838
-0.433012701892

Example 28.7, page no. 732

In [15]:
import numpy

print "Here cˆ2=16, taking h=1, finding k such that cˆ2tˆ2=1 "
print "So bendre−schmidits recurrence relation ie u(i)(j+1)=(16tˆ2(u(i−1)(j)+u(i+1)(j))+2(1−16∗tˆ2u(i,j)−u(i)(j−1)"
print "Now since u(0,t)=0=u(5,t) therefore u(0,i)=0 and u(5,j)=0 and u(x,0)=xˆ2(5−x)"
c = 4
h =1 
k = (h/c)
t = k/h
A = numpy.zeros((6,6))
print "Also from 1st derivative (u(i)(j+1)−u(i,j−1))/2k=g(x) and g(x)=0 in this case"
print "So if j=0 this gives u(i)(1)=1/2∗(u(i−1)(0)+u(i+1)(0))"
for i in range(0,6):
    for j in range(1,9):
        A[0,i] = 0
        A[5,i] = 0;
        A[i,1] = (i)**2*(5-i)
for i in range(0,4):
    A[i+1,2] =1/2*(A[i,1]+A[i+2,1])
for i in range(2,5):
    for j in range(2,5):
        A[i-1,j] = (c*t)**2*(A[i-2,j-1]+A[i,j-1])+2*(1-(c*t)**2)*A[i-1,j-1]-A[i-1,j-2]
for i in range(0,5):
    for j in range(0,5):
        print  "  ",A[i,j],
    print ""

Here cˆ2=16, taking h=1, finding k such that cˆ2tˆ2=1 
So bendre−schmidits recurrence relation ie u(i)(j+1)=(16tˆ2(u(i−1)(j)+u(i+1)(j))+2(1−16∗tˆ2u(i,j)−u(i)(j−1)
Now since u(0,t)=0=u(5,t) therefore u(0,i)=0 and u(5,j)=0 and u(x,0)=xˆ2(5−x)
Also from 1st derivative (u(i)(j+1)−u(i,j−1))/2k=g(x) and g(x)=0 in this case
So if j=0 this gives u(i)(1)=1/2∗(u(i−1)(0)+u(i+1)(0))
   0.0    0.0    0.0    0.0    0.0 
   0.0    4.0    8.0    12.0    16.0 
   0.0    12.0    24.0    36.0    48.0 
   0.0    18.0    36.0    54.0    72.0 
   0.0    16.0    0.0    0.0    0.0

Example 28.8, page no. 734

In [1]:
import numpy

print "Here cˆ2=4, taking h=1, finding k such that cˆ2tˆ2=1 "
print "So bendre−schmidits recurrence relation ie u(i)(j+1)=(16tˆ2(u(i−1)(j)+u(i+1)(j))+2(1−16∗tˆ2u(i,j)−u(i)(j−1)"
print "Now since u(0,t)=0=u(4,t) therefore u(0,i)=0 and u(4,j)=0 and u(x,0)=x(4−x) "
c = 2
h = 1
k = (h/c)
t = k/h
A = numpy.zeros((6,6))
print "Also from 1st derivative (u(i)(j+1)−u(i,j−1))/2 k=g(x) and g(x)=0 in this case"
print "So if j=0 this gives u(i)(1)=1/2∗(u(i−1)(0)+u(i+1)(0))"
for i in range(0,6):
    for j in range(1,9):
        A[0,i] = 0
        A[4,i] = 0
        A[i,0] = (i)*(4-i)
for i in range(0,4):
    A[i+1,2] = 1/2*(A[i,1]+A[i+2,1])
for i in range(2,5):
    for j in range(2,5):
        A[i-1,j] = (c*t)**2*(A[i-2,j-1]+A[i,j-1])+2*(1-(c*t)**2)*A[i-1,j-1]-A[i-1,j-2]
for i in range (0,5):
    for j in range(0,5):
        print "  ",A[i,j],
    print ""

Here cˆ2=4, taking h=1, finding k such that cˆ2tˆ2=1 
So bendre−schmidits recurrence relation ie u(i)(j+1)=(16tˆ2(u(i−1)(j)+u(i+1)(j))+2(1−16∗tˆ2u(i,j)−u(i)(j−1)
Now since u(0,t)=0=u(4,t) therefore u(0,i)=0 and u(4,j)=0 and u(x,0)=x(4−x) 
Also from 1st derivative (u(i)(j+1)−u(i,j−1))/2 k=g(x) and g(x)=0 in this case
So if j=0 this gives u(i)(1)=1/2∗(u(i−1)(0)+u(i+1)(0))
   0.0    0.0    0.0    0.0    0.0 
   3.0    0.0    -3.0    -6.0    -9.0 
   4.0    0.0    -4.0    -8.0    -12.0 
   3.0    0.0    -3.0    -6.0    -9.0 
   0.0    0.0    0.0    0.0    0.0