Chapter 5 - Boundary Value Problems

Exa 5.1 Page 87

In [1]:
from __future__ import division

#finite difference method
#given d2y_by_dx2-2=0 hence it is dirichlet's problem

x_1=0 ; y_1=0               #initial boundary conditions
x_3=1 ; y_3=0

delta_x=.5                #since we have to find y_2 at x=.5 so x_2=.5 
#in the central difference method substituting i=2 we have
#(y_3-2*y_2+y_1)/(delta_x)**2=2
#since y_1 & y_3=0 as per B.C.
y_2=(y_3+y_1-2*delta_x**2)/2
print "the value of y at x=.5 from finite difference method is",y_2
x=.5
y_exact=x**2-x
print "the value of y from exact soln at x=.5 is",y_exact

the value of y at x=.5 from finite difference method is -0.25
the value of y from exact soln at x=.5 is -0.25

Exa 5.2 Page 89

In [5]:
from __future__ import division
#boundary conditions are: x=0 at y=0# dy/dx=1 at x=1
print "to solve this problem we will take delta x=.5 since we have to find the value at x=.5"
delta_x=.5
y_1=0
#using central difference eqn
dy_by_dx=1                 #at x=1, i=3
#y_4=dy/dx*2*delta_x+y_2             sincefrom B.C. at node 3

#y_2=delta_x**2+y_3-delta_x           on substituting the value of y_4
y_3=-(2*delta_x**2+2*(delta_x**2-delta_x)-y_1) #on substituting the value of y_2
y_2=delta_x**2+y_3-delta_x
print "the value of y at x=.5 is",y_2

to solve this problem we will take delta x=.5 since we have to find the value at x=.5
the value of y at x=.5 is -0.25

Exa 5.3 Page 89

In [6]:
from __future__ import division
from numpy import arange
print "the solution of eg 5.3 -->Discretization in 1-D space"
#given the source term f(x)=4x**2-2x-4
#given eqn d2y/dx2-2y=f(x)
y_1=0
y_4=-1
delta_x=1/3           #since given 3 parts and length=1
for i in range(0,4):
    j=arange(0,1+delta_x,delta_x)

#given to divide the line in 3 parts
#at node 2
#y_3-2*y_2
def fx3(x):
    d=(4*x**2-2*x-4)
    return d
f2=fx3(j[1])
f3=fx3(j[2])
y_3=((f2)*delta_x**2+(2+2*delta_x**2)*((f3)*delta_x**2-y_4)-y_1)/(1-(2+2*delta_x**2)**2)
y_2=(f3+2*y_3)*delta_x**2+2*y_3-y_4
print "the value of y at x=",j[1],j[2],"is respectively",y_2,y_3

the solution of eg 5.3 -->Discretization in 1-D space
the value of y at x= 0.333333333333 0.666666666667 is respectively 0.111111111111 -0.222222222222

Exa 5.4 Page 90

In [7]:
from __future__ import division
from numpy import arange,zeros
y_1=0
dy_by_dx=-3            #at x=1
delta_x=1/3           #since given 3 parts and length=1
for i in range(0,4):
    j=arange(0,1+delta_x,delta_x)

#given to divide the line in 3 parts
def fx3(x):
    d=(4*x**2-2*x-4)
    return d

f2=fx3(j[1])
f3=fx3(j[2])
f4=fx3(j[3])
print "the exact analytical soln are"
y=[]
for i in range(1,4):
    y.append(-2*j[i]**2+j[i])
    print y[i-1]

#using B.C.-2 at node 4 we get
#(y_5-y_3)/(2*delta_x)=-3
#eqn at node 2
#-20*y_2+9*y_3=f2
#at node 3
#9*y_2-20*y_3+9*y_4=f3
#at node 4
#18*y_3-20*y_4=16
print "the value of f(x) at node 2 3 and 4 are",f2,f3,f4
#now solving the equations using TDMA method 

a=[0,9,18]                   #sub diagonal assignment

b=[]
for j in range(1,4):
    b.append(-20)            #main diagonal assignment
c=[]
for k in range(1,3):
    c.append(9)            #super diagonal assignment

d=[f2,f3,16]                     #given values assignment

i=1#
n=3#
beta1=[b[i-1]]#                #initial b is equal to beta since a1=0
gamma1=[d[i-1]/beta1[i-1]]       #since c7=0
m=i+1#
for j in range(m,n+1):
    beta1.append(b[j-1]-a[j-1]*c[j-2]/beta1[j-2])
    gamma1.append((d[j-1]-a[j-1]*gamma1[j-2])/beta1[j-1])
x=zeros(n)
x[-1]=gamma1[n-1]               #since c7=0
n1=n-i#
for k in range(1,n1+1):
    j=n-k# 
    x[j-1]=gamma1[j-1]-c[j-1]*x[j]/beta1[j-1]

print "the values of y2, y3, y4 by TDMA method are"
for i in range(1,4):
    print x[i-1]

the exact analytical soln are
0.111111111111
-0.222222222222
-1.0
the value of f(x) at node 2 3 and 4 are -4.22222222222 -3.55555555556 -2.0
the values of y2, y3, y4 by TDMA method are
0.111111111111
-0.222222222222
-1.0

Exa 5.5 Page 91

In [1]:
from __future__ import division
from math import cosh
from numpy import zeros
delta_x=.1
y_11=1
dy_by_dx_1=0            #dy/dx at x=0
# given eqn d2y/dx2=y
print "the analytical soln are"
y_an=[]
for i in range(1,11):
    y_an.append(cosh((i-1)/10)/cosh(1))
    print y_an[i-1]

#using central difference method we have 
#y(i-1) - (2+delat_x**2)y(i) + y(i+1)=0
#therefore the eqns can be solved using TDMA method
a=[None]
for i in range(2,11):
    a.append(1)            #sub diagonal assignment
b=[]
for j in range(1,11):
    b.append(-2.01)           #main diagonal assignment

c=[2]
for k in range(2,10):
    c.append(1)            #super diagonal assignment
d=[]
for l in range(1,10):
    d.append(0)
                 #given values assignment
d.append(-1)
i=1#
n=10#
beta1=[b[i-1]]                #initial b is equal to beta since a1=0
gamma1=[d[i-1]/beta1[i-1]]      #since c7=0
m=i+1#
for j in range(m,n+1):
    beta1.append(b[j-1]-a[j-1]*c[j-2]/beta1[j-2])
    gamma1.append((d[j-1]-a[j-1]*gamma1[j-2])/beta1[j-1])

x=zeros(10)
x[n-1]=gamma1[n-1]                #since c7=0
n1=n-i#
for k in range(1,n1+1):
    j=n-k
    x[j-1]=gamma1[j-1]-c[j-1]*x[j]/beta1[j-1]

print "the values of y from y1 to y10 by TDMA method are"
for i in range(1,11):
    print x[i-1]

the analytical soln are
0.648054273664
0.651297246159
0.661058620401
0.677436091507
0.70059357071
0.730762825846
0.768245800962
0.81341763827
0.8667304327
0.92871775662
the values of y from y1 to y10 by TDMA method are
0.648259699273
0.65150099777
0.661257306244
0.67762618778
0.700771331194
0.730924187921
0.768386286526
0.813532247997
0.866813531948
0.928762951218

Exa 5.6 Page 92

In [1]:
from __future__ import division
from numpy import zeros
T_1=100 ;  T_10=200
delta_x=(T_10-T_1)/9               #since 9 divisions are to be made
#using central difference method we get
#for node 2--> 2*T_2-T_3=100
a=[0]
for i in range(2,9):
    a.append(-1)            #sub diagonal assignment
b=[]
for j in range(1,9):
    b.append(2)           #main diagonal assignment
c=[]
for k in range(1,8):
    c.append(-1)            #super diagonal assignment

d=[100]
for l in range(2,8):
    d.append(0);         #given values assignment

d.append(200)


i=1;
n=8;
beta1=[b[i-1]]                #initial b is equal to beta since a1=0
gamma1=[d[i-1]/beta1[i-1]];      #since c7=0
m=i+1;
for j in range(m,n+1):
    beta1.append(b[j-1]-a[j-1]*c[j-2]/beta1[j-2])
    gamma1.append((d[j-1]-a[j-1]*gamma1[j-2])/beta1[j-1])

x=zeros(n)    
x[-1]=gamma1[n-1]               #since c7=0
n1=n-i;
for k in range(1,n1+1):
    j=n-k
    x[j-1]=gamma1[j-1]-c[j-1]*x[j]/beta1[j-1];

print "the values of T from T2 to T9 by TDMA method are"
for i in range(1,9):
    print x[i-1]

the values of T from T2 to T9 by TDMA method are
111.111111111
122.222222222
133.333333333
144.444444444
155.555555556
166.666666667
177.777777778
188.888888889

Exa 5.7 Page 93

In [2]:
from __future__ import division
from math import pi,sqrt
from numpy import zeros
dia=.02
l=.05
T_0=320
delta_x=l/4
k=50
h=100
T_surr=20
#B.C--> d(theta)/dx+h(theta)/k=0 at x=0.05
#let m=sqrt(hP/kA)
P=pi*dia
A=pi*dia**2/4
m=sqrt(h*P/(k*A))#
#using central difference method we get 
#-theta(i-1)+(2+(m*delta_x**2)*theta(i)+theta(i+1))=0
theta_0=T_0-T_surr
#using B.C. at node 4 we get--> theta(5)=theta(3)-2h*theta(4)*delta_x/k
#now the eqns can be solved using TDMA method
a=[None]
for i in range(2,4):
    a.append(-1 )           #sub diagonal assignment

a.append(-2)
b=[]
for j in range(1,4):
    b.append(2.0625)#           #main diagonal assignment

b.append(2.1125)
c=[]
for k in range(1,4):
    c.append(-1)            #super diagonal assignment
d=[300]
for l in range(2,5):
    d.append(0)
                  #given values assignment

i=1#
n=4#
beta1=[b[i-1]]                #initial b is equal to beta since a1=0
gamma1=[d[i-1]/beta1[i-1]]      #since c7=0
m=i+1
for j in range(m,n+1):
    beta1.append(b[j-1-1]-a[j-1]*c[j-2-1]/beta1[j-1-1])
    gamma1.append((d[j-1]-a[j-1]*gamma1[j-2])/beta1[j-1])

x=zeros(n)    
x[n-1]=gamma1[n-1]#               #since c7=0
n1=n-i#
for k in range(1,n1+1):
    j=n-k#
    x[j-1]=gamma1[j-1]-c[j-1]*x[j]/beta1[j-1]

print "the values of T from T1 to T4 in Celsius by TDMA method are"
for i in range(1,5):
    print x[i-1]-T_surr

the values of T from T1 to T4 in Celsius by TDMA method are
231.893171899
199.529667042
180.886766374
174.799288605

Exa 5.8 Page 94

In [1]:
from __future__ import division
from math import pi,sqrt
from numpy import zeros

lnght=.001
k_const=.001
l=.05
D=10**-9
delta_x=l/100
C=1                 #in mol/m3
#B.C. are C=1 at x=0
#         dC/dx=0 at x=10**-3    since at the end point conc. is const.
#using central difference method we get the following eqns which can be solved using TDMA method
a=[0]
for i in range(2,100):
    a.append(1)            #sub diagonal assignment

a.append(2)                  #since C99=C100 using B.C. 

b=[]
for j in range(1,101):
    b.append(-2.0001)           #main diagonal assignment

c=[]
for k in range(1,100):
    c.append(1)            #super diagonal assignment

d=[-1]
for l in range(2,101):
    d.append(0)
          #given values assignment
i=1#
n=100#
beta1 =[b[i-1]]                #initial b is equal to beta since a1=0
gamma1=[d[i-1]/beta1[i-1]]      #since c7=0
m=i+1#
for j in range(m,n+1):
    beta1.append(b[j-1]-a[j-1]*c[j-2]/beta1[j-2])
    gamma1.append((d[j-1]-a[j-1]*gamma1[j-2])/beta1[j-1])

x=zeros(n)    
x[n-1]=gamma1[n-1]               #since c7=0
n1=n-i#
for k in range(1,n1+1):
    j=n-k# 
    x[j-1]=gamma1[j-1]-c[j-1]*x[j]/beta1[j-1]

print "the values of conc. at x=.5mm or at the 50th node is",x[49]

the values of conc. at x=.5mm or at the 50th node is 0.730764441227