Numerical Methods of Solving BVP of differential
equations
Ayan Panja (M23MA2004)
m23ma2004@iitj.ac.in
November 2023
1
Contents
1 Introduction 2
2 Boundary value Problem 2
2.1 Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2 Corollary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
3 Solutions of BVP 3
3.1 SHOOTING METHOD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
3.1.1 EXAMPLE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
3.2 FINITE DIFFERENCE METHOD . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.2.1 EXAMPLE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.3 CUBIC SPLINE METHOD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.3.1 EXAMPLE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
3.4 GALERKIN’s METHOD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3.4.1 EXAMPLE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
4 CONCLUSION 12
5 REFERENCES 12
6 ACKNOWLEDGEMENT 13
1 Introduction
Ordinary Differential Equations are frequently occur in mathematical models that arise in many
branches of sciences,engineering and economics.Unfortunately it is seldom that these equations
have solutions,which can be expressed in closed forms ,so it is common to seek approximation
solutions by means of numerical methods.Nowadays this can be usually achieved very inexpen-
sively to high accuracy and with a reliable bound on the error between the analytical solutions
and its numerical approximation . In this section we shall be concerned with the construction
and the analysis of numerical methods of solving BVP differential equations.
2 Boundary value Problem
we have seen that we require m conditions to be specified in order to solve m-order differential
equation.It is not necessary to specify these conditions at one point of the independent variable.
They can be specified at different points in the interval (a,b) and, therefore such problems are
called the boundary value problems. A large no of problem fall into this category.
In case of boundary value problem, we seek solutions at specified points within the domain
of given boundaries, for instance,given:
d2 y
= f (x, y, y ′ ) (1)
dx2
y(a)=ya
y(b) = yb
Here we are interested in finding the values of y in the range x∈ (a, b).
2
2.1 Theorem
Suppose the function f in the boundary value problem
y ′′ = f (x, y, y ′ ) (2)
for x∈ (a, b)
y(a) = α
y(b) = β
is continuous on the set D=(x,y,y’) where x∈ (a, b)y, y ′ ∈ (−∞, ∞) and that the partial
derivatives with respect to y and y’ are also continuous on D if
1. f(x,y,y’)than0 for all (x,y,y’)∈ D
2. Let ∃M
M being a constant | f( y ′ )(x, y, y ′ ) |> 0f orall(x, y, y ′ ) ∈ D.
the BVP(2) has unique solution.
2.2 Corollary
Suppose the linear boundary-value problem
y ′′ = p(x)y ′ + q(x)y + r(x) (3)
for a≤ x ≤ b
y(a) = α
y(b) = β
satisfies
1. p(x),q(x),r(x) are continuous on [a,b].
2. q(x) is greater than 0 on [a,b]
Then the BVP(3) has unique solutions.
3 Solutions of BVP
Here I discuss about some methods to solve boundary value problems :
1. SHOOTING METHOD
2. FINITE DIFFERENCE METHOD
3. CUBIC SPLINE METHOD
4. GALERKIN’s METHOD
3
Figure 1: Illustration of SHOOTING METHOD
3.1 SHOOTING METHOD
This method is called the the shooting method because it resembles an artillery problem.In
this method,the given boundary value problem is first converted into an equivalent initial value
problem and then solve it using any numerical methods for solving IVP using numerical meth-
ods.
The approach is very simple.
Consider equation
y ′′ = f (x, y, y ′ ), y(a) = A, y(b) = B (4)
letting y’=z
we obtain the following set of two equations
y ′ = z; z ′ = f (x, y, z) (5)
In order to solve this set as an IVP, we need two conditions at x=a.We have one conditions
y(a)=A and y(b)=B and therefore,require another condition for z at x=a.Let us assume that
z(a)=M1 , where
M1
is a ”guess” . Note that it represents the slope y’(x) at x=a.
Thus the problem reduced to a system of two first order equations with the initial conditions
y ′ = z, y(a) = A; z ′ = f (x, y, z), z(a) = M1 (6)
Equation(6) can be solved for y and z using one step method using steps of n. Until the solution
at x=b is reached.
Let the estimated value of y(x) at x=b be B1
If B1 = B,
then we have obtained the required solution.
In practice, it is very unlikely that our initial guess z(a)=M1
is correct.
If B1 ̸= B
4
then we obtain the solution with another guess, say z(a)=M2 .
Let the new estimate of y(x) at x=b be B2 , If B2 ̸= B
then the process may be continued until we obtain the correct estimate of y(b).Let us assume
that z(a)=M3 ,
leads to the value y(b)=B . If we assume that the values of M and B are linearly related,
then
M3 − M2 M2 − M1
= (7)
B − B2 B2 − B1
B2 − B
M3 = M2 − (M2 − M1 ) (8)
B2 − B1
3.1.1 EXAMPLE
Using shooting method , to solve the equation
d2 y
= 6x; y(1) = 2, y(2) = 9. (9)
dx2
in the interval(1,2).Find y(1.5).
SOLUTION :
By the transformation we obtain the following dydz
=z
y(1) = 2
dz
dx
= 6x.
Let us assume that z(1)=y’(1)=2(M1 ).
Applying Heun’s method , we obtain the solution as follows :
ITERATION I
h = 0.5
x0 = 1, y(1) = y0 = 2, z(1) = z0 = 2
m1 (1) = z0 = 2
m1 (2) = 6x0 = 6
m2 (1) = z0 + hm1 (2) = 2 + 0.5(6) = 5
m2 (2) = 6(x0 + h) = 6(1.5) = 9
m(1) = m1 (1)+m 2
2 (1)
= 2+5
2
= 3.5
m1 (2)+m2 (2) 6+9
m(2) = 2
= 2 = 7.5
y(x1 ) = y(1.5) = y(1) + m(1)h = 2 + 3.5 · 0.5 = 3.75
z(x1 ) = z(1.5) = z(1) + m(2)h = 2 + 7.5 · 0.5 = 5.75.
ITERATION II
h = 0.5
x1 = 1.5, y1 = 3.75, z1 = 5.75
m1 (1) = z1 = 5.75
m1 (2) = 6x1 = 9
m2 (1) = z1 + hm1 (2) = 5.75 + 0.5(9) = 10.25
m2 (2) = 6(x1 + h) = 12
m(1) = 5.75+10.25
2
=8
m(2) = 9+122
= 10.5
y(x2 ) = y(2) = y1 + m(1)h = 3.75 + 8 · 0.5 = 7.75.
This gives B1 = 7.75.
which is less than B=9.
Now,let us assume z(1)=y’(1)=4(M2 ).
and again estimate y(2)
5
ITERATION I
h = 0.5, x0 = 1, y0 = 2, z0 = 2
m1 (1) = z0 = 4
m1 (2) = 6x0 = 6
m2 (1) = z0 + hm1 (2) = 4 + 0.5(6) = 7
m2 (2) = 6(x0 + h) = 6(1.5) = 9
m(1) = 4+72
= 5.5
6+9
m(2) = 2 = 7.5
y(x1 ) = y(1.5) = 2 + 5.5 · 0.5 = 4.75
z(x1 ) = z(1.5) = 4 + 7.5 · 0.5 = 7.75.
ITERATION II
h = 0.5, x1 = 1.5, y1 = 4.75, z1 = 7.75
m1 (1) = z1 = 7.75
m1 (2) = 6x1 = 9
m2 (1) = z1 + hm1 (2) = 7.75 + 0.5(9) = 12.25
m2 (2) = 6(x1 + h) = 12
m(1) = 7.75+12.25
2
= 10
9+12
m(2) = 2 = 10.5
y(x2 ) = y(2) = y1 + m(1)h = 4.75 + 10 · 0.5 = 9.75.
This gives B2 = 9.75.
which is greater than B=9
Now, let us have the third estimate of z(1)=M3 .
Using relationship :
−B
M3 = M2 − BB22−B 1
(M2 − M1 )
9.75−9
M3 = 4 − 9.75−9 (4 − 2) = 4 − 0.75 = 3.25
The new estimation for z(1)=y’(1)=3.25
ITERATION I
h = 0.5, x0 = 1, y0 = 2, z0 = 3.25
m1 (1) = z0 = 3.25
m1 (2) = 6x0 = 6
m2 (1) = z0 + hm1 (2) = 3.25 + 0.5(6) = 6.25
m2 (2) = 6(x0 + h) = 6(1.5) = 9
m(1) = 3.25+6.25
2
= 4.75
m(2) = 6+92
= 7.5
y(x1 ) = y(1.5) = 2 + 4.75 · 0.5 = 4.375
z(x1 ) = z(1.5) = 3.25 + 7.5 · 0.5 = 7.
ITERATION II
h = 0.5, x1 = 1.5, y1 = 3.75, z1 = 7
m1 (1) = z1 = 7
m1 (2) = 6x1 = 9
m2 (1) = z1 + hm1 (2) = 7. + 0.5(9) = 11.5
m2 (2) = 6(x1 + h) = 12
m(1) = 7.75+12.25
2
= 9.25
m(2) = 9+12
2
= 10.5
y(x2 ) = y(2) = 4.375 + 9.25 · 0.5 = 9.
The solution is y(1)=2,y(1.5)=4.375,y(2)=9.
The exact solution is y(x)=x3 + 1.
Hence therefore y(1.5)=4.375
6
Figure 2: Solution of DE by FINITE DIFFERENCE method
3.2 FINITE DIFFERENCE METHOD
In this method the derivatives are replaced by their finite equivalent, thus converting the
differential equation into a system of algebraic equations.For example we can use the following
central difference approximations:
y( i + 1) − y( i − 1)
yi′ = (10)
2h
y( i + 1) − 2yi − y( i − 1)
yi′′ = (11)
h2
These are second-order equations and the accuracy estimates can be improved by using higher
order equations.
The given interval (a,b) can be divided into n sub-intervals each of width h.Then
xi = x0 + ih = a + ih
yi = y(xi ) = y(a + ih)
y0 = y(a)
yn = y(a + nh) = y(b)
This is illustrated in Figure 2. The differential equation is written for each of the inter-
nal points i=1,2,...,n-1.If differential equation is linear,this will be result with n-1 unknowns
y1 , y2 , ..., yn − 1.
We can solve for these unknowns using any of the estimation methods.
3.2.1 EXAMPLE
Solve the boundary value problem
d2 y
−y =0 (12)
dx2
with y(0)=0 and y(2)=3.62686
The exact solution of this problem is y=sinh(x). Find y(1) ?
7
For this we divide the interval [0,2] into two sub-intervals with h=1.0 . The difference
equation at the single unknown point y1 is :
Y0 − 2Y1 + Y2 = Y1 .
Using the values of y0 andy2 .
we obtain y1 = 1.202895.
But the exact value of y=sinh(x) at x=1 is 1.17520. So the error is 0.03375 .So it is not a
good approximation
Now we subdivide the interval [0,2] into four equal parts so that h=0.5.Let the values of y
at the five points be y0 , y1 , y2 , y3 , y4 .
We are given that y0 = 0, y4 = 3.62686
Writing the difference equations at three interval point, we obtain
4(y0 − 2y1 + y2 ) = y1 , 4(y1 − 2y2 + y3 ) = y2 , 4(y2 − 2y3 + y4 ) = y3 (13)
respectively.Substituting for y0 andy4 .
After rearranging,we get the system
−9y1 + 4y2 = 0, 4y1 − 9y2 + 4y3 = 0, 4y2 − 9y3 = −14.50744 (14)
The solution of (13) is y1 = y(0.5) = 0.52635, y2 = y(1.0) = 1.18428, y3 = y(1.5) = 2.13829
Which is a better approximation since the error is now reduced to 0.00908
3.3 CUBIC SPLINE METHOD
Here we consider the boundary value problem
y ′′ (x) + f (x)y ′ (x) + g(x)y(x) = r(x) (15)
with the boundary conditions
y(x0 ) = A, y(xn ) = B (16)
Let s(x) be the cubic spline approximating the function y(x) and let s”(xi ) = Mi .At, x = xi
equation(15) becomes
Mi + fi s′ (xi ) + gi yi = ri (17)
But
h 1
s′ (xi −) = (2Mi + M( i − 1)) + (yi − y( i − 1)) (18)
3! h
and
h 1
s′ (xi +) = − (2Mi + M( i − 1)) + (y( i + 1) − yi ) (19)
3! h
substituting equation (18) and (19) successively in equation (17) , we obtain the equations
8
h 1
Mi + fi [ (2Mi + M( i − 1)) + (yi − y( i − 1))] + gi yi = ri (20)
3! h
and
h 1
Mi + fi [− (2Mi + M( i − 1)) + (y( i + 1) − yi )] + gi yi = ri (21)
3! h
since y0 , yn
are known ,equation (20) and (21) constitute a system of 2n equations in 2n unknowns,viz.,M0 , M1 ,
M2 , ..., Mn , y1 , y2 , ..., y( n − 1).
It is,however,possible to eliminate the Mi
and obtain a traditional system for yi .
3.3.1 EXAMPLE
Given the boundary value problem
x2 y ′′ + xy − y = 0; y(1) = 1, y(2) = 0.5 (22)
apply the cubic spline method to determine the value of y(1.5).
The given differential equation is
1 1
y ′′ = − y ′ + 2 y. (23)
x x
Setting x=xi , y ′′ (xi ) = Mi , Eq.(23)becomes.
1 ′ 1
Mi = − yi + 2 yi . (24)
xi xi
Using the expressions given in Eqs. (18) and (19) ,we obtain
1 h h y( i + 1) − yi 1
Mi = − [− Mi − M( i + 1)) + ] + 2 yi ; i = 0, 1, 2, ..., n − 1. (25)
xi 3 6 h xi
and
1 h h yi − y( i − 1) 1
Mi = − [− Mi − M( i − 1)) + ] + 2 yi ; i = 1, 2, ..., n. (26)
xi 3 6 h xi
If we divide [1,2] into two sub-intervals, we have h=0.5 and n=2.Then Eqs. (25) (26) give
10M0 − M1 + 24y1 = 36; 16M1 − M2 − 32y1 = −12; (27)
M0 + 20M1 + 16y1 = 24; M1 + 26M2 − 24y1 = −9 (28)
Eliminating M0 , M1 , M2 .
from these system of equations we obain
y1 = 0.65599.
Since the exact value of y1 = y(1.5) = 23 .
we can take y1 = 0.65599.
9
3.4 GALERKIN’s METHOD
This method uses trial functions (or approximating functions) which satisfy the boundary
conditions of the problem.The trial function is substituted in the given differential equation
problem and the result is called residual.The integral of the product of this residual and a
weighted function, taken over the domain,is then set to zero which yields a system of equations
for the unknown parameters in the trial functions.
Let the boundary value problem be defined by
y ′′ + p(x)y ′ + q(x)y = f (x), a < x < b (29)
with the boundary conditions
p0 y(a) + q0 y ′ (a) = r0 ; p1 y(b) + q1 y ′ (b) = rb (30)
Let the approximate solution given by
n
X
t(x) = αi ϕi (x)
i=1
(31)
where ϕi (x)
are called base functions.Substituting t(x) in equation (29), we obtain a residual.Denoting
this residual by R(t),
we obtain
R(t) = t′′ + p(x)t′ + q(x)t − f (x) (32)
Usually the base functions ϕi (x)
are chosen as weight functions.We,therefore,have
Z b
I= ϕi (x)R(t) dx = 0 (33)
a
which yields a system of equations for the parameters αi ,
which are known,t(x) can be calculated from equation (31)
NOTE: This method also called the WEIGHTED RESIDUAL Method.
3.4.1 EXAMPLE
Solve the boundary value problem defined by
y ′′ + y ′ + x = 0; y(0) = y(1) = 0. (34)
where 0 ∈ (0, 1)
Find y(0.5).
10
Let t(x)=α1 ϕ1 (x).
Since the both boundary conditions must be satisfied by t(x) , we can choose
ϕ1 (x) = x(1 − x).Substitutingf ort(x)inthegivendif f erentialequation, weobtain
R(t)=t”+t+x.
Hence we have
Z 1
I= (t′′ + t + x)α1 x(1 − x) dx = 0 (35)
0
Z 1
(t′′ + t + x)x(1 − x) dx = 0 (36)
0
Now,
R 1 ′′
doing integration by parts
0 (t + t + x)x(1 − x) dx
1 R1 ′
= [t′ x(1 − x)] − 0 t (1 − 2x) dx
0
− 1 t′ (1 − 2x) dx
R
= 0
= −[[t(1 − 2x)])10 − 01 t(−2)dx]
R
= −2 01 tdx,
R
since t=0 at x=0 and x=1.
Hence (36) Rsimplifies to
− 0 tdx + 01 tx(1 − x)dx + 01 x2 (1 − x) = 0
R1 R
putting the value of t,
3 4 1
- 0 α1 x(1 − x)dx + 01 α1 x2 (1 − x)d x + [ x3 − x4 ] = 0.
R1 R
0
5
Hence α1 = 18 = 0.2777778, simplied
Then a first approximation to the solution is
5
y(0.5) = 18 (0.5)(0.5) = 0.06944
The exact solution to the given boundary value problem is
y(x) = sinx
sin1
− x,
which means that our solution has an error of 0.0003.
11
4 CONCLUSION
5 REFERENCES
1. Introductory Methods of Numerical Analysis,5th ed. S.S. Sastry
2. Numerical Methods, E Balagurusamy(2011)
3. An Introduction to Numerical Analysis, Endre Süli and David Mayers(2008)
4. Numerical Analysis, 9th ed. Richard L. Burden and J.Douglas Faires
5. https://youtu.be/1QbYfKJHrJk?si=z6UkQ6H02-ZiZQEY
6. https://youtu.be/ud09P64-5fM?si=9iUOG4BeLCdpMeoO
12
6 ACKNOWLEDGEMENT
I would like to express my sincere gratitude to all those who have contributed to the success
for this report.I am thankful for the support and guidance of Kirankumar R. Hiremath sir for
sharing his knowledge and expertise in the subject matter,which helped me to shape my ideas
and concepts.I am thankful to my classmates for their suggestions that helped me to improve
my work.
Thank you all for your guidance,support,and encouragement.
13
