Boundary Value Problems for ODE

Basic Definitions

We wish to solve the differential equation

y'' + ay' + by = c
in the interval t_i £ t £ t_f. The prime denotes defferentiation of y with respect to t. We need to specify two conditions for a complete solution of this problem.

Initial Value Problem (IVP):

Both y and y' are given at the initial instant of time t_i, as follows

y( t_i ) = f, y'( t_i ) = g

Boundary Value Problem (BVP):

Both y and y' are not given at the initial instant, instead the conditions are mixed between the boundaries t_f and t_i of the range of t. Three such boundary conditions are shown below.

(a) y( t_i ) = f, y( t_f ) = g
(b) y( t_i ) = f, y'( t_f ) = g
(c) y'( t_i ) = f, y( t_f ) = g

Shooting Method

We will describe the method in detail for the type (a) problem, other situations are similar.

First, we convert the 2nd-order ODE into two 1st-order ODE's, as follows

y₁' = y₂
y₂' = c - ay₂ - by₁
with boundary conditions

y₁( t_i ) = f, y₁( t_f ) = g . . . Eq. (A)
where y₁ = y, and y₂ = y'.

We cannot use Runge-Kutta (RK) type marching solution technique for this problem, bacause RK-type methods require that the intital conditions are specified. Thus we convert the BVP into an IVP. As the inital value of y₂ is not given, we guess this value. Thus, instead of the conditions in (A), we use the conditions

y₁( t_i ) = f, y₂( t_i ) = h . . . Eq. (B)
With the initial conditions of (B), we can use the RK-type marching method to solve the set of ODE's. However, this solution will not satisfy the condition y₁( t_f ) = g. Instead, we will get y₁( t_f ) = K. Note that K depends on the choice of h. Thus, K is a function of h. What we want is to find a suitable h such that

K( h ) = g or K( h ) - g = 0 . . . Eq. (C)

The Eq.(C) is solved by using the Newton-Raphson method. Begin with a guess value h₀ for h, and update h using the formula

h_n+1 = h_n - [ K( h_n ) - g ] / K'( h_n )
The derivative K'( h_n ) can be obtained as follows

K'( h_n ) = [ K( h_n + D h ) - K( h_n ) ] / D h
A suitable D h can be 0.01 or 0.001.

Examples by using SIMULINK

Convective Heat Transfer Along a Rod ( D h = 0.001 )

Iteration	h_n	K ( h_n )	K( h_n + D h )	K'( h_n )	h_n+1
0	5	109.622	109.633	11.7	12.7246
1	12.7246	200.401	200.410	9.0	12.6800
2	12.6800	199.877	199.889	11.8	12.6904

Radiative Heat Transfer Along a Rod ( D h = 0.001 )

Iteration	h_n	K ( h_n )	K( h_n + D h )	K'( h_n )	h_n+1
0	5	98.8492	98.8637	14.5	11.9759
1	11.9759	506.841	507.187	346.0	11.0891
2	11.0891	326.674	326.796	122.1	10.0516
3	10.0516	238.397	238.457	60.6	9.4180
4	9.4180	205.343	205.388	44.5	9.2979
5	9.2979	200.1412	200.184	42.2	9.2945