Chapter 15 - Solution of partial differential equations

Example No. 15_01 Pg No. 488

In [11]:
from numpy import zeros,mat
from numpy.linalg import solve
#Elliptic Equations

l = 15
h = 5
n = 1 + 15/5
f=zeros([4,4])
f[0,0:3]=100
f[0:3,0]=100

#At point 1 :  f2 + f3 - 4f1 + 100 + 100 = 0
#At point 2 :  f1 + f4 - 4f2 + 100 +   0 = 0
#At point 3 :  f1 + f4 - 4f3 + 100 +   0 = 0
#At point 4 :  f2 + f3 - 4f4 +   0 +   0 = 0
#
#Final Equations are
#         -4f1 + f2 + f3 + 0 = -200
#           f1 - 4f2 + 0 + f4 = -100
#           f1 + 0 - 4f3 + f4 = -100
#           0  + f2 + f3 - 4f4 = 0
A = mat([[ -4, 1 ,1 ,0],[1, -4, 0 ,1],[1, 0 ,-4 ,1],[0, 1, 1, -4 ]])
B = mat([[-200],[-100],[-100],[0]])
C = solve(A,B)
print '\n The solution of the system is \n f1 = %d \n f2 = %d \n f3 = %d \n f4 = %d \n '%(C[0],C[1],C[2],C[3])

 The solution of the system is 
 f1 = 75 
 f2 = 50 
 f3 = 50 
 f4 = 24

Example No. 15_02 Pg No. 489

In [1]:
from numpy import zeros,mat
from numpy.linalg import solve

#Liebmann's Iterative method

f=zeros([4,4])
f[0,0:3]=100
f[0:3,0]=100
A=zeros([5,5])
A[0,0:5] = range(0,5)
for n in range(1,6):
    for i in range(2,4):
        for j in range(2,4):
            if n == 1 and i == 2 and j == 2 :
                f[i-1,j-1] = ( f[i,j] + f[i-2,j-2] + f[i-2,j] + f[i,j-2] )/4
            else:
                f[i,j] = ( f[i,j-1] + f[i-2,j-1] + f[i-1,j]+ f[i-1,j-2] )/4
    A[1:5,n-1] = mat([f[1,1],f[1,2],f[2,1],f[2,2] ])


print 'First row of below matrix represents iteration number:\n',A

First row of below matrix represents iteration number:
[[  0.   1.   2.   3.   4.]
 [ 75.  75.  75.  75.  75.]
 [  0.   0.   0.   0.   0.]
 [  0.   0.   0.   0.   0.]
 [  0.  50.  50.  50.  50.]]

Example No. 15_03 Pg No. 490

In [5]:
from numpy import zeros,mat
from numpy.linalg import solve

#Poisson's Equation

#D2f = 2*x**2 * y**2
# f = 0
# h = 1 
#Point 1 : 0 + 0 + f2 + f3 - 4f1 = 2(1)**2 * 2**2
#          f2 + f3 - 4f1 = 8
#Point 2 : 0 + 0 + f1 + f4 -4f2 = 2*(2)**2*2**2
#           f1 - 4f2 = f4 = 32
#Point 3 : 0 + 0 + f1 + f4 - 4f4 = 2*(1**2)*1**2
#           f1 -4f3 + f4 = 2
#Point 4 : 0 + 0 + f2 + f3 - 4f4 = 2* 2**2 * 1**2
#            f2 + f3 - 4f4 = 8
#Rearranging the equations
#            -4f1 + f2 + f3 = 8
#             f1 - 4f2 + f4 = 32
#             f1 - 4f3 + f4 = 2
#              f2 + f3 - 4f4 = 8
A = mat([[ -4, 1, 1, 0],[1, -4 ,0 ,1],[1, 0, -4, 1],[0, 1, 1, -4]])
B = mat([[ 8],[32],[2],[8 ]])
C = solve(A,B)
print 'The solution is \n f1 = %f \n f2 = %f \n f3 = %f \n f4 = %f \n '%( C[0],C[1],C[2],C[3])

The solution is 
 f1 = -5.500000 
 f2 = -10.750000 
 f3 = -3.250000 
 f4 = -5.500000

Example No. 15_04 Pg No. 491

In [1]:
#Gauss-Seidel Iteration

f2 = 0
f3 = 0
for i in range(1,5):
    f1 = (f2 + f3 - 8)/4
    f4 = f1 
    f2 = (f1 + f4 -32)/4
    f3 = (f1 + f4 - 2)/4
    print '\n Iteration %d\n f1 = %f,    f2 = %f,    f3 = %f,    f4 = %f\n'%(i,f1,f2,f3,f4)

 Iteration 1
 f1 = -2.000000,    f2 = -9.000000,    f3 = -2.000000,    f4 = -2.000000


 Iteration 2
 f1 = -5.000000,    f2 = -11.000000,    f3 = -3.000000,    f4 = -5.000000


 Iteration 3
 f1 = -6.000000,    f2 = -11.000000,    f3 = -4.000000,    f4 = -6.000000


 Iteration 4
 f1 = -6.000000,    f2 = -11.000000,    f3 = -4.000000,    f4 = -6.000000

Example No. 15_05 Pg No. 494

In [1]:
from numpy import zeros,mat

#Initial Value Problems

h = 1 #
k = 2 #
tau = h**2/(2*k)
f=zeros([7,5])
for i in range(1,5):
    f[0,i] = 50*( 4 - (i) )

f[0:7,0] = 0
f[0:7,4] = 0 
for j in range(0,6):
    for i in range(1,4):
        f[j+1,i] = ( f[j,i-1] + f[j,i+1] )/2 
print 'The final results are \n',f

The final results are 
[[   0.   150.   100.    50.     0. ]
 [   0.    50.   100.    50.     0. ]
 [   0.    50.    50.    50.     0. ]
 [   0.    25.    50.    25.     0. ]
 [   0.    25.    25.    25.     0. ]
 [   0.    12.5   25.    12.5    0. ]
 [   0.    12.5   12.5   12.5    0. ]]

Example No. 15_06 Pg No. 497

In [6]:
from numpy import zeros,mat,linalg
#Crank-Nicholson Implicit Method

h = 1 #
k = 2 #
tau = h**2/(2*k)
f=zeros([5,5])
for i in range(1,4):
    f[0,i] = 50*( 4 - (i) )
f[0:5,0] = 0 
f[0:5,4] = 0 
A = mat([[4,  -1,  0 ],[-1 , 4 , -1 ],[0 , -1 , 4]])
B=zeros([3,1])
for j in range(0,4):
    for i in range(1,4):
        B[i-1,0] = f[j,i-1] + f[j,i+1]
    
    C = linalg.solve(A,B)
    f[j+1,1] = C[0]
    f[j+1,2] = C[1]
    f[j+1,3] = C[2]

print 'The final solution using crank nicholson implicit method is :\n',f

The final solution using crank nicholson implicit method is :
[[   0.          150.          100.           50.            0.        ]
 [   0.           42.85714286   71.42857143   42.85714286    0.        ]
 [   0.           26.53061224   34.69387755   26.53061224    0.        ]
 [   0.           13.70262391   20.11661808   13.70262391    0.        ]
 [   0.            7.70512287   10.70387339    7.70512287    0.        ]]

Example No. 15_07 Pg No. 500

In [11]:
from __future__ import division
from numpy import sqrt
#Hyperbolic Equations

h = 1
Tbyp = 4
tau = sqrt(h**2 /4)
r = int(1+(2.5 - 0)/tau)
c = int(1+(5 - 0)/h)

f=zeros([6,6])
for i in range(1,int(c)-1):
    f[0,i] = (i)*(5 - (i) )

f[0:r-1,0] = 0
f[0:r-1,c-1] = 0
g=[]
for i in range(1,c-1):
    g.append(0)
    f[1,i] = (f[0,i+1] + f[0,i-1])/2 + tau*g[i-1] 

for j in range(1,r-1):
    for i in range(1,c-1):
        f[j+1,i] = -f[j-1,i] + f[j,i+1] + f[j,i-1]
    

print 'The values estimated are :\n',f

The values estimated are :
[[ 0.  4.  6.  6.  4.  0.]
 [ 0.  3.  5.  5.  3.  0.]
 [ 0.  1.  2.  2.  1.  0.]
 [ 0. -1. -2. -2. -1.  0.]
 [ 0. -3. -5. -5. -3.  0.]
 [ 0. -4. -6. -6. -4.  0.]]