Numerical Solution of Boundary

Value Problem by Finite Difference

Method

July 26, 2024

2

Consider 2nd order linear ODE:

d2 y dy
2
+ p(x) + q(x)y = r(x) (1)
dx dx
with given boundary conditions: y(x=a)=α, y(x=b)=β
We are to solve y(x) in the range a≤x≤b by finite difference method.
We do uniform discretion of x within a≤x≤b by N intervals =⇒
(N+1) points.



 x0 = a

x1 = a + h





x = a + 2h
2 (b − a)
(N + 1)points h= (2)

 ................... N

xj = a + jh → f orj = 0, 1, . . . , N





xN = a + N h


The Centre Difference Formula gives;

dy 1
( ) xi = (yi+1 − yi−1 ) (3)
dx 2h
d2 y 1
( 2
)xi = 2 [yi+1 − 2yi + yi−1 ] (4)
dx h
Therefore, discretized version of equation (1) at i-th point would be,
1 1
2
(yi+1 − 2yi + yi−1 ) + pi (yi+1 − yi−1 ) + qi yi = ri (5)
h 2h
h h
=⇒ (1− pi )yi−1 +(h2 qi −2)yi +(1+ pi )yi+1 = h2 ri → For i=1,2,3,. . . ,(N-1)
2 2
(6)
 3

So, we have (N-1) number of simultaneous linear algebraic equations for

(N-1) unknowns: y1 , y2 , y3 , . . . , yN −1
From boundary conditions, we know that, y0 = α & yN = β.
The simultaneous equations:

[ci yi−1 + di yi + ei yi+1 = ri h2 ] (7)


h
ci ≡ (1 − 2 pi )

where, di ≡ (h2 qi − 2) (8)

ei ≡ (1 + h2 pi )


From equation (7):



 i = 1 → c1 y0 + d1 y1 + e1 y2 = h2 r1 ; y0 is known.

i = 2 → c2 y1 + d2 y2 + e2 y3 = h2 r2




i = 3 → c y + d y + e y = h2 r

3 2 3 3 3 4 3
F or(N −1)equations


 .......................................................
i = (N − 2) → cN −2 yN −3 + dN −2 yN −2 + eN −2 yN −1 = h2 rN −2





i = (N − 1) → cN −1 yN −2 + dN −1 yN −1 + eN −1 yN = h2 rN −1 ; yN is known.

(9)
In the form AY=B:
  
d1 e1 0 0 0 ... ... 0 y1  
 c2 d 2 e 2 h2 r1 − c1 y0
0 0 ... ... 0   y2  
  

 0 c3 d 3 e 3 h2 r2 
0 ... ... 0   y3  
2

    h r 
 0 0 c4 d4 e4
 ... ... 0   ...  = 
   
..
3 

 .. .. .. .. .. .. .. ..     . 
. . . . . . . .  yN −3   2


 0 0 . . . . . . . . . cN −2 dN −2 eN −2  yN −2 
   h r N −2 
2
h rN −1 − eN −1 yN
0 0 . . . . . . . . . . . . cN −1 dN −1 yN −1 | {z }
| {z } | {z } B(N −1)×1
A(N −1)×(N −1) Y
(10)

=⇒ a00 = d1 

=⇒ a11 = d2

(N − 1) diagonal elements.
.................... 

=⇒ aN −2,N −2 = dN −1


=⇒ a01 = e1 

=⇒ a12 = e2

(N −2) elements right side of diagonal elements.
.................... 

=⇒ aN −3,N −2 = eN −2

4

And,

=⇒ a10 = c2 

=⇒ a21 = c3

(N −2) elements lef t side of diagonal elements.
.................... 

=⇒ aN −2,N −3 = cN −1


Therefore in compact form:


aii = di+1 f or i = 0 to N − 2

ai,i+1 = ei+1 f or i = 0 to N − 3 All other elements are zero. (11)

ai,i−1 = ci+1 f or i = 1 to N − 2


Problem set on Finite Difference Method

1. Consider the ODE: y” + 2y ′ + y = 0 with given boundary
conditions: y(−0.5) = 0 and y(4) = 9e−4
(a) Solve the ODE in the range −0.5 ≤ x ≤ 4 using finite difference
method.
(b) This ODE has the analytic solution: y(x) = (2x + 1)e−x . Plot
the numerical solution along with analytic solution.
(c) Let, represent the numerical solution points by ỹi (where i
runs over the discrete points). The deviation of numerical
solution ỹ from analytic solution y for N points can be defined by,
v
u
u1 X N
σ(N ) = t (ỹi − yi )2 (12)
N i=1

Here, plot σ(N ) w.r.t. N taking different values of N.

Input:

import matplotlib.pyplot as plt

import numpy as np

#Number of Points=N
#Number of discrete uniform intervals=(N-1)
N=22

#(a)
A=[]
 5

A=np.zeros((N-2,N-2))

xi,xf=-0.5,4
yi=0
yf=9*(np.exp(-4))

h=(xf-xi)/(N-1)

xx=[]
pp,qq,rr=[],[],[]
cc,dd,ee=[],[],[]
for i in range(1,N-2+1):
x=xi+i*h
p=2
c=1-((h/2)*p)
e=1+((h/2)*p)
q=1
d=((h**2)*q)-2
r=0
cc.append(c)
dd.append(d)
ee.append(e)
pp.append(p)
qq.append(q)
rr.append(r)
xx.append(x)

for i in range(0,N-3):
A[i+1][i]=cc[i+1] #A[1][0]=C2=cc[1]
for i in range(0,N-2):
A[i][i]=dd[i] #A[0][0]=d1=dd[0]
for i in range(0,N-3):
A[i][i+1]=ee[i] #A[0][1]=e1=ee[0]

B=[]
B=np.zeros((N-2,1))
B[0][0]=((h**2)*rr[0])-(cc[0]*yi)
for i in range(1,N-3):
B[i][0]=(h**2)*rr[i]
B[N-3][0]=((h**2)*rr[N-3])-(ee[N-3]*yf)
6

Y=[]
Y=np.linalg.solve(A,B)
print(Y)

#(b)
xx1=[]
yy=[]
for i in range(N):
x1=xi+i*h
y=((2*x1)+1)*(np.exp(-x1))
xx1.append(x1)
yy.append(y)

plt.plot(xx,Y,’ok’,label=’Numerical Solution’)
plt.plot(xx1,yy,’-k’,label=’Analytic Solution’)
plt.title(’1.(b) Numerical vs Analytic Solution’)
plt.legend()
plt.grid()
plt.show()

#(c)
n=N-2
zz=[]
z=0
for i in range(n):
z=z+(Y[i]-yy[i+1])**2
zz.append(z)
ss=[]
nn=[]
for i in range(n):
s=((1/(i+1))*zz[i])**0.5
ss.append(s)
nn.append(i)
plt.plot(nn,ss,’.-k’)
plt.xlabel(’No. of Points(n)’)
plt.ylabel(’Sigma(n)’)
plt.title(’1.(c) Deviation of Numerical Solution vs No. of Points’)
plt.grid()
plt.show()
 7

Output:
8

2. Consider the ODE: y” + ex y ′ + y = 2(x − 1)e−x − x + 1 with given

boundary conditions: y(0) = 0 and y(4) = 4e−4
(a) Solve the ODE in the range 0 ≤ x ≤ 4 using finite difference
method.
(b) This ODE has the analytic solution: y(x) = xe−x . Plot the
numerical solution along with analytic solution.
Input:

import matplotlib.pyplot as plt

import numpy as np

#Number of Points=N
#Number of discrete uniform intervals=(N-1)
N=17

#(a)
A=[]
A=np.zeros((N-2,N-2))

xi,xf=0,4
yi=0
yf=4*(np.exp(-4))

h=(xf-xi)/(N-1)

xx=[]
pp,qq,rr=[],[],[]
cc,dd,ee=[],[],[]
for i in range(1,N-2+1):
x=xi+i*h
p=np.exp(x)
c=1-((h/2)*p)
e=1+((h/2)*p)
q=1
d=((h**2)*q)-2
r=(2*(x-1)*(np.exp(-x)))-x+1
cc.append(c)
dd.append(d)
ee.append(e)
pp.append(p)
qq.append(q)
 9

rr.append(r)
xx.append(x)

for i in range(0,N-3):
A[i+1][i]=cc[i+1] #A[1][0]=C2=cc[1]
for i in range(0,N-2):
A[i][i]=dd[i] #A[0][0]=d1=dd[0]
for i in range(0,N-3):
A[i][i+1]=ee[i] #A[0][1]=e1=ee[0]

B=[]
B=np.zeros((N-2,1))
B[0][0]=((h**2)*rr[0])-(cc[0]*yi)
for i in range(1,N-3):
B[i][0]=(h**2)*rr[i]
B[N-3][0]=((h**2)*rr[N-3])-(ee[N-3]*yf)

Y=[]
Y=np.linalg.solve(A,B)
print(Y)

#(b)
xx1=[]
yy=[]
for i in range(N):
x1=xi+i*h
y=x1*(np.exp(-x1))
xx1.append(x1)
yy.append(y)

plt.plot(xx,Y,’ok’,label=’Numerical Solution’)
plt.plot(xx1,yy,’-k’,label=’Analytic Solution’)
plt.title(’2.(b) Numerical & Analytic Solution’)
plt.legend()
plt.grid()
plt.show()
10

Output:

3. Consider the ODE: y” + xy ′ − 4y = 14x2 − 2 with given boundary

conditions: y(−1) = y(1) = 0.
(a) Solve the ODE in the range −1 ≤ x ≤ 1 using finite difference
method.
(b) This ODE has the analytic solution: y(x) = x4 − x2 . Plot the
numerical solution along with analytic solution.
Input:

import matplotlib.pyplot as plt

import numpy as np

#Number of Points=N
#Number of discrete uniform intervals=(N-1)
N=22

#(a)
A=[]
A=np.zeros((N-2,N-2))
 11

xi,xf=-1,1
yi=0
yf=0

h=(xf-xi)/(N-1)

xx=[]
pp,qq,rr=[],[],[]
cc,dd,ee=[],[],[]
for i in range(1,N-2+1):
x=xi+i*h
p=x
c=1-((h/2)*p)
e=1+((h/2)*p)
q=-4
d=((h**2)*q)-2
r=(14*(x**2))-2
cc.append(c)
dd.append(d)
ee.append(e)
pp.append(p)
qq.append(q)
rr.append(r)
xx.append(x)

for i in range(0,N-3):
A[i+1][i]=cc[i+1] #A[1][0]=C2=cc[1]
for i in range(0,N-2):
A[i][i]=dd[i] #A[0][0]=d1=dd[0]
for i in range(0,N-3):
A[i][i+1]=ee[i] #A[0][1]=e1=ee[0]

B=[]
B=np.zeros((N-2,1))
B[0][0]=((h**2)*rr[0])-(cc[0]*yi)
for i in range(1,N-3):
B[i][0]=(h**2)*rr[i]
B[N-3][0]=((h**2)*rr[N-3])-(ee[N-3]*yf)

Y=[]
12

Y=np.linalg.solve(A,B)
print(Y)

#(b)
xx1=[]
yy=[]
for i in range(N):
x1=xi+i*h
y=(x1**4)-(x1**2)
xx1.append(x1)
yy.append(y)

plt.plot(xx,Y,’ok’,label=’Numerical Solution’)
plt.plot(xx1,yy,’-k’,label=’Analytic Solution’)
plt.title(’3.(b) Numerical & Analytic Solution’)
plt.legend()
plt.grid()
plt.show()
Output:
 13

4. Consider the ODE: (1 − x2 )y” − 2xy ′ + 6y = 0 with given

boundary conditions: y(−0.95) = y(0.95) = 0.85375
(a) Solve the ODE in the range −0.95 ≤ x ≤ 0.95 using finite
difference method.
(b) This ODE has the analytic solution: y(x) = (3x2 − 1)/2. Plot
the numerical solution along with analytic solution.
Input:

import matplotlib.pyplot as plt

import numpy as np

#Number of Points=N
#Number of discrete uniform intervals=(N-1)
N=22

#(a)
A=[]
A=np.zeros((N-2,N-2))

xi,xf=-0.95,0.95
yi,yf=0.85375,0.85375

h=(xf-xi)/(N-1)

xx=[]
pp,qq,rr=[],[],[]
cc,dd,ee=[],[],[]
for i in range(1,N-2+1):
x=xi+i*h
p=(-2*x)/(1-(x**2))
c=1-((h/2)*p)
e=1+((h/2)*p)
q=6/(1-(x**2))
d=((h**2)*q)-2
r=0
cc.append(c)
dd.append(d)
ee.append(e)
pp.append(p)
qq.append(q)
rr.append(r)
14

xx.append(x)

for i in range(0,N-3):
A[i+1][i]=cc[i+1] #A[1][0]=C2=cc[1]
for i in range(0,N-2):
A[i][i]=dd[i] #A[0][0]=d1=dd[0]
for i in range(0,N-3):
A[i][i+1]=ee[i] #A[0][1]=e1=ee[0]

B=[]
B=np.zeros((N-2,1))
B[0][0]=((h**2)*rr[0])-(cc[0]*yi)
for i in range(1,N-3):
B[i][0]=(h**2)*rr[i]
B[N-3][0]=((h**2)*rr[N-3])-(ee[N-3]*yf)

Y=[]
Y=np.linalg.solve(A,B)
print(Y)

#(b)
xx1=[]
yy=[]
for i in range(N):
x1=xi+i*h
y=(3*(x1**2)-1)/2
xx1.append(x1)
yy.append(y)

plt.plot(xx,Y,’ok’,label=’Numerical Solution’)
plt.plot(xx1,yy,’-k’,label=’Analytic Solution’)
plt.title(’4.(b) Numerical & Analytic Solution’)
plt.legend()
plt.grid()
plt.show()
 15

Output:

5. Consider the ODE: (1 − x2 )y” − 2xy ′ + 12y = 0 with given

boundary conditions: y(−0.95) = −y(0.95) = −0.7184375
(a) Solve the ODE in the range −0.95 ≤ x ≤ 0.95 using finite
difference method.
(b) This ODE has the analytic solution: y(x) = (5x3 − 3x)/2. Plot
the numerical solution along with analytic solution.
Input:

import matplotlib.pyplot as plt

import numpy as np

#Number of Points=N
#Number of discrete uniform intervals=(N-1)
N=17

#(a)
A=[]
A=np.zeros((N-2,N-2))
16

xi,xf=-0.95,0.95
yi,yf=-0.7184375,0.7184375

h=(xf-xi)/(N-1)

xx=[]
pp,qq,rr=[],[],[]
cc,dd,ee=[],[],[]
for i in range(1,N-2+1):
x=xi+i*h
p=(-2*x)/(1-(x**2))
c=1-((h/2)*p)
e=1+((h/2)*p)
q=12/(1-(x**2))
d=((h**2)*q)-2
r=0
cc.append(c)
dd.append(d)
ee.append(e)
pp.append(p)
qq.append(q)
rr.append(r)
xx.append(x)

for i in range(0,N-3):
A[i+1][i]=cc[i+1] #A[1][0]=C2=cc[1]
for i in range(0,N-2):
A[i][i]=dd[i] #A[0][0]=d1=dd[0]
for i in range(0,N-3):
A[i][i+1]=ee[i] #A[0][1]=e1=ee[0]

B=[]
B=np.zeros((N-2,1))
B[0][0]=((h**2)*rr[0])-(cc[0]*yi)
for i in range(1,N-3):
B[i][0]=(h**2)*rr[i]
B[N-3][0]=((h**2)*rr[N-3])-(ee[N-3]*yf)

Y=[]
Y=np.linalg.solve(A,B)
 17

print(Y)

#(b)
xx1=[]
yy=[]
for i in range(N):
x1=xi+i*h
y=(5*(x1**3)-(3*x1))/2
xx1.append(x1)
yy.append(y)

plt.plot(xx,Y,’ok’,label=’Numerical Solution’)
plt.plot(xx1,yy,’-k’,label=’Analytic Solution’)
plt.title(’5.(b) Numerical & Analytic Solution’)
plt.legend()
plt.grid()
plt.show()

Output:
18

6. Consider the ODE: (1 − x2 )y” − 2xy ′ + 20y = 0 with given

boundary conditions: y(−0.95) = y(0.95) = 0.5540898438
(a) Solve the ODE in the range −0.95 ≤ x ≤ 0.95 using finite
difference method.
(b) This ODE has the analytic solution: y(x) = (35x4 − 30x2 + 3)/8.
Plot the numerical solution along with analytic solution.
Input:

import matplotlib.pyplot as plt

import numpy as np

#Number of Points=N
#Number of discrete uniform intervals=(N-1)
N=32

#(a)
A=[]
A=np.zeros((N-2,N-2))

xi,xf=-0.95,0.95
yi,yf=0.5540898438,0.5540898438

h=(xf-xi)/(N-1)

xx=[]
pp,qq,rr=[],[],[]
cc,dd,ee=[],[],[]
for i in range(1,N-2+1):
x=xi+i*h
p=(-2*x)/(1-(x**2))
c=1-((h/2)*p)
e=1+((h/2)*p)
q=20/(1-(x**2))
d=((h**2)*q)-2
r=0
cc.append(c)
dd.append(d)
ee.append(e)
pp.append(p)
qq.append(q)
rr.append(r)
 19

xx.append(x)

for i in range(0,N-3):
A[i+1][i]=cc[i+1] #A[1][0]=C2=cc[1]
for i in range(0,N-2):
A[i][i]=dd[i] #A[0][0]=d1=dd[0]
for i in range(0,N-3):
A[i][i+1]=ee[i] #A[0][1]=e1=ee[0]

B=[]
B=np.zeros((N-2,1))
B[0][0]=((h**2)*rr[0])-(cc[0]*yi)
for i in range(1,N-3):
B[i][0]=(h**2)*rr[i]
B[N-3][0]=((h**2)*rr[N-3])-(ee[N-3]*yf)

Y=[]
Y=np.linalg.solve(A,B)
print(Y)

#(b)
xx1=[]
yy=[]
for i in range(N):
x1=xi+i*h
y=(35*(x1**4)-(30*(x1**2))+3)/8
xx1.append(x1)
yy.append(y)

plt.plot(xx,Y,’ok’,label=’Numerical Solution’)
plt.plot(xx1,yy,’-k’,label=’Analytic Solution’)
plt.title(’6.(b) Numerical & Analytic Solution’)
plt.legend()
plt.grid()
plt.show()
20

Output:

The End