3 Shooting Methods for Boundary Value Problems

3.1 Shooting methods for boundary value problems with linear ODEs

Shooting methods are developed to transform boundary value problems (BVPs) for ordinary differential equations to an equivalent initial value problem (IVP). The term "shooting method" is inspired by the problem illustrated in Figure 30, where the problem is to "shoot" a ballistic object in the field of gravity, aiming to hit a target at a given length $ L $. The initial angle $ \alpha $ is then changed in an repeated fashion, based on observed lengths, until the target at $ L $ is hit.

Figure 30: The trajectory of a ballistic object launched with an initial angle $ \alpha $.

The problem is a boundary value problem with one condition given at $ x= 0 $ and another at $ x=L $. It might seems obvious that the problem of varying $ \alpha $ until the condition at $ x=L $ is satisfied has a solution. We will use this approach on problems which has nothing to do with ballistic trajectories.

Let us consider the following example: $$ \begin{equation} y'' = y(x) \tag{3.1} \end{equation} $$ which is a second order, linear ODE with initial conditions: $$ \begin{equation} y(0)=0, \ y'(0) =s \tag{3.2} \end{equation} $$

and consequently an initial value problem which can be shown to have the following analytical solution: $$ \begin{equation} y(x)=s \cdot \sinh(x) \tag{3.3} \end{equation} $$

For each choice of the initial value $ y'(0)=s $ we will get a new solution $ y(x) $, and in Figure 31 we see the solutions for $ s= 0.2 $ and $ 0.7 $.

The problem we really is (3.1) with the following boundary conditions: $$ \begin{equation} y(0)=0, \ y(1) =1 \tag{3.4} \end{equation} $$ which is what we call a boundary value problem.

From (3.3) we realize that the boundary problem may be solve by selecting $ s=s^* $ such that $ y(1)=s^*\cdot \sinh(1) $ or $$ \begin{equation} s^*= \frac{1}{\sinh(1)} \tag{3.5} \end{equation} $$

For this particular case we are able to find the analytical solution of both the initial problem and the boundary value problem. When this is not the case our method of approach is a shooting method where we select values of $ s $ until the condition $ y(1)=1 $ is fulfilled. For arbitrary values of $ s $ the boundary value $ y(1) $ becomes a function of $ s $, which is linear if the ODE is linear and nonlinear if the ODE is nonlinear.

To illustrate the shooting method, consider a somewhat general boundary problem for a second order linear ODE: $$ \begin{equation} y''(x) = p(x) \cdot y'(x) + q(x) \cdot y(x) + r(x) \tag{3.6} \end{equation} $$

where $ p(x) $, $ q(x) $ are arbitrary functions and $ r(x) $ may be considered as a source term. The boundary values may be prescribed as: $$ \begin{equation} y(a)=\alpha,\ y(b)=\beta \tag{3.7} \end{equation} $$

The second order ODE (3.6) may be reduced to a system of two ODEs (see 2.4 Reduction of Higher order Equations): $$ \begin{align} y'(x) & = g(x) \nonumber \\ g'(x) & = p(x) \cdot g(x) + q(x) \cdot y(x) + r(x) \tag{3.8} \end{align} $$ with the boundary conditions given in (3.7).

We start the shooting method by choosing the given boundary value $ y(a)=\alpha $ as an initial condition. As we need an initial condition for $ y'(\alpha) \equiv g(\alpha) $ to solve (3.8) as an initial value problem, we have to guess and initial value in a way such that the boundary value $ y(b)=\beta $ is satisfied. As (3.6) is a linear ODE it suffices to guess two values of $ s=g(\alpha) \equiv y'(a) $. The correct value for the initial value of $ s $ which gives the correct value of the solution $ y(b)=\beta $, may then be found by linear interpolation. Note that $ y(x) $ is always proportional to $ s $ when the ODE is linear.

To quantify the goodness of how well the boundary value resulting from our initial guess $ s=g(a) \equiv y'(a) $, we introduce a boundary value error function $ \phi $: $$ \begin{equation} \phi (s) =y(b;s)-\beta \tag{3.9} \end{equation} $$

The correct initial guess $ s= s^* $ is found when the boundary value error function is zero: $$ \begin{equation} \phi (s^*)=0 \tag{3.10} \end{equation} $$

The procedure may be outlined as follows:

The shooting method for linear boundary value problems.

Guess two initial values $ s^0 $ and $ s^1 $
Compute the corresponding values for the error function $ \phi ^0 $ and $ \phi ^1 $ by solving the initial value problem in (3.8)
Find the correct initial value $ s^* $ by linear interpolation.
Compute the solution to the boundary value problem by solving the initial value problem in (3.8) with $ g(a)=y'(a)=s^* $.

To see how we find $ s^* $ from linearization, consider first $ \phi $ as a linear function of $ s $: $$ \begin{equation} \phi=k_a\cdot s+k_b \tag{3.11} \end{equation} $$

where $ k_a $ is the slope and $ k_b $ is the offset when $ s=0 $. We easily see from (3.11) that $ \phi=0 $ for $ s^*=-\dfrac{k_b}{k_a} $. When we have two samples $ \phi ^0 $ and $ \phi ^1 $ of $ \phi $ at $ s^0 $ and $ s^1 $, both $ k_a $ and $ k_b $ may be found from (3.11) to be: $$ \begin{equation} k_a=\frac{\phi ^1-\phi ^0}{s^1-s^0}, \qquad k_b= \frac{s^1\cdot \phi ^0 - \phi ^1 \cdot s^0}{s^1-s^0} \tag{3.12} \end{equation} $$

and consequently we may find $ s^* $ as: $$ \begin{equation} s^*= \frac{\phi^1 s^0-\phi^0 s^1}{\phi^1-\phi^0} \tag{3.13} \end{equation} $$

which also may be presented on an incremental form as: $$ \begin{equation} \tag{3.14} s^*=s^1+\Delta s, \qquad \Delta s = -\phi ^1 \cdot \left(\frac{s^1-s^0}{\phi^1 - \phi^0}\right) \end{equation} $$

Let us go back to our model boundary value example (3.1) with boundary values as prescribed in (3.4) which may be presented as reduced system of first order ODEs: $$ \begin{equation} \tag{3.15} \left.\begin{matrix} y'(x)=g(x) \\ g'(x)=y(x) \end{matrix}\right.\ \end{equation} $$

For the given values of $ b $ and $ \beta $ and $ s=y'(0)=g(0) $ the boundary value error function in (3.9) takes the form: $$ \begin{equation} \phi (s)=y(1;s)-1 \tag{3.16} \end{equation} $$

In accordance with the shooting method we choose two initial values $ s^0=0.2 $ and $ s^1=0.7 $ and solve (3.15) as an initial value problem with e.g. and RK4-solver with $ \Delta x=0.1 $ and obtain the following results:.

$ m $	$ s^m $	$ \phi(s^m) $
0	0.2	-0.7650
1	0.7	-0.1774

Note that we use $ m $ as a superindex for iterations, which will be used in case of nonlinear equations. Substitution into (3.13) gives the $ s^*=0.8510 $. Subsequent use of the $ s^* $-value in the RK4-solver yields $ \phi (0.8510)=0.0001 $. We observe a good approximation to the correct value for the initial value may be found from the analytical solution as: $ y'(0)=\dfrac{1}{\sinh(1)}=0.8509 $

Figure 31: The solution $ y(x) $ of a boundary value problem resulting from two initial guesses $ y(x;s^0) $ and $ y(x;s^1) $.

The presentation above with the shooting method is chosen as it can easily be generalized to the solution of nonlinear ODEs.

Alternative approach for linear second order ODEs.

For linear second order ODEs the solution may be found in somewhat simpler fashion, by solving the following sub-problems:

Sub-problem 1 $$ \begin{align} y_0''(x)&=p(x) \cdot y_0'(x) + q(x) \cdot y_0(x) + r(x), & \qquad y_0(a)=\alpha , \quad y_0'(a)=0 \tag{3.17} \end{align} $$ and

Sub-problem 2 $$ \begin{align} y_1''(x)&=p(x) \cdot y_1'(x) + q(x) \cdot y_1(x), & \qquad y_1(a)=0 , \quad y_1'(a)=1 \tag{3.18} \end{align} $$

That is, the two sub-problems differ only by the source term $ r(x) $ and the boundary conditions. Notice that the condition $ y'(\alpha)=0 $ in (3.17) corresponds to $ s^0=0 $ and the condition $ y'(\alpha)=1 $ i (3.18) corresponds to $ s^1=1 $.

Let $ y_0(x) $ represent the solution to (3.17) and $ y_1(x) $ be the solution to (3.18). The complete solution of the boundary value problem in (3.6) with the boundary conditions (3.7) may then be shown to be: $$ \begin{equation} y(x)=y_0(x) + \left[\frac{\beta - y_0(b)}{y_1(b)}\right]\cdot y_1(x) = y_0(x) - \left[\frac{\phi^0}{\phi^1+\beta}\right]\cdot y_1(x) \tag{3.19} \end{equation} $$

with $ s^0=0 $ and $ s^1=1 $. Note that the solution of (3.17) corresponds to a particular solution of (3.1) and (3.4), and that the solution of (3.18) corresponds to the homogenous solution of (3.1) and (3.4) as the source term is missing. From general theory on ODEs we know that the full solution may be found by the sum of a particular solution and the homogenous solution, given that the boundary conditions are satisfied.