5.2. The shooting method

Consider the two-point boundary value problem

\[ y'' = f(t,y), \qquad t \in [t_0, t_{\max}], \qquad y(t_0) = a ,\qquad y(t_{\max}) = b. \]

To solve a second-order boundary value problem numerically, we first rewrite it as a system of first-order ODEs. This allows us to apply the initial value problem solvers developed in earlier chapters.

Let

\[\begin{split} \begin{align*} y_1' &= y_2 ,\\ y_2' &= f(t, y_1, y_2). \end{align*} \end{split}\]

Then the boundary condition \(y(t_0) = a\) gives \(y_1(t_0) = a\), but the initial value \(y_2(t_0) = y'(t_0)\) is unknown. The shooting method treats this missing value as an unknown parameter \(s\). For a chosen value of \(s\), we solve the resulting initial values problem and compare the computed value of \(y_1(t_{\max})\) with the required boundary condition \(b\). The value of \(s\) is then adjusted until the boundary condition at \(t = t_{\max}\) is satisfied.

The method is called the shooting method because it resembles aiming at a target. If the computed solution misses the required boundary value, the initial slope is adjusted and another “shot” is taken until the target is hit.

Example 5.2

Use the Euler method with a step length of \(h = 0.2\) and the shooting method to solve the following boundary value problem

\[ \begin{align*} y'' - y' - y &= 0, & y(0) &= 0, & y(1) &= 2. \end{align*} \]

---

Solution

Rewriting this as a system of first-order ODEs gives

\[\begin{split} \begin{align*} y_1' &= y_2, & y_1 (0) &= 0, &&\\ y_2' &= y_1 +y_2, & y_2 (0) &= s, && \end{align*} \end{split}\]

where \(s\) is a guess of \(y'(0)\).

The Euler method solutions using a guess value of \(s = 1\) is tabulated below.

\(t\)	\(y_1\)	\(y_2\)
0.00	0.000000	1.000000
0.20	0.200000	1.200000
0.40	0.440000	1.480000
0.60	0.736000	1.864000
0.80	1.108800	2.384000
1.00	1.585600	3.082560

The numerical solution at the upper boundary is \(y_1(1) = 1.585600\). Lets try again but with an increased guess value of \(s = 2\).

\(t\)	\(y_1\)	\(y_2\)
0.00	0.000000	2.000000
0.20	0.400000	2.400000
0.40	0.880000	2.960000
0.60	1.472000	3.728000
0.80	2.217600	4.768000
1.00	3.171200	6.165120

Here we have \(y_1(1) = 3.171200\). Since the solution obtained with \(s\) gives \(y_1(1) < 2\) and the solution obtained with \(s=2\) gives \(y_1(1) > 2\), the correct value of \(s\) must lie between 1 and 2. Since the desired initial slope lies between 1 and 2, a natural next choice is \(s = 1.5\).

\(t\)	\(y_1\)	\(y_2\)
0.00	0.000000	1.500000
0.20	0.300000	1.800000
0.40	0.660000	2.220000
0.60	1.104000	2.796000
0.80	1.663200	3.576000
1.00	2.378400	4.623840

The solutions using the three guess values of \(s=1\), \(s=2\) and \(s=1.5\) are plotted below.

Fig. 5.1 Solutions of the BVP \(y'' - y' - y = 0\), \(y(0) = 0\), \(y(1) = 2\) using initial values for \(y'(0) = s\) where \(s = 1\), \(s = 2\) and \(s = 1.5\).

The approximation obtained using \(s=1.5\) is closer to satisfying the boundary condition at \(t=1\) than the previous guesses. In practice, the shooting method is usually combined with a root-finding algorithm such as the secant method to determine the value of \(s\) automatically.

5.2.1. Improving the guess value using the secant method

In Example 5.2, successive guesses for the unknown initial slope produced numerical solutions that moved closer to the required boundary value at \(t = t_{\max}\). One approach would be to repeatedly halve the interval containing the correct value of \(s\), a strategy known as the bisection method. A more efficient alternative is the secant method which uses previous approximations to predict and improved value of \(s\).

The goal of the shooting method is to determine the value of \(s = y'(t_0)\) that satisfies the boundary condition

\[ y(t_{\max}) = b. \]

Define the residual function

\[ g(s) = b - y(t_{\max}, s), \]

where \(y(t_{\max}, s)\) denotes the numerical solution obtained using the initial slope \(s\).

If \(g(s) = 0\), then the computed solution satisfies the boundary condition exact. Therefore, solving the boundary value problem reduces to finding the root of \(g(s)\).

The secant method generates a series of approximations

\[s_0, \; s_1, \; s_2, \; \ldots \]

that converge to a root of \(g(s) = 0\). Given two previous approximations \(s_{n - 1}\) and \(s_{n - 2}\), the next approximation is computed using

(5.4)\[ s_{n+1} = s_n - g(s_n)\frac{s_n - s_{n-1}}{g(s_n) - g(s_{n-1})}. \]

The iteration continues until successive approximations differ by less than a prescribed tolerance, that is,

\[ |s_n - s_{n-1}| < tol. \]

Example 5.3

Use the secant method to calculate the next guess value \(s\) for the solution of the boundary value problem Example 5.2 and calculate the solution to the BVP using this guess value.

---

Solution

From Example 5.2 we have solutions at the upper boundary of

\[\begin{split} \begin{align*} y_1(1) &= 3.1712, & s_1 &= 2, \\ y_1(1) &= 2.3784, & s_2 &= 1.5. \end{align*} \end{split}\]

The values of \(g(s_1)\) and \(g(s_2)\) are

\[\begin{split} \begin{align*} g(s_1) &= 2 - 3.1712 = -1.1712, \\ g(s_2) &= 2 - 2.3784 = -0.3784, \end{align*} \end{split}\]

therefore, using equation (5.4) the next guess value calculated using the secant method is

\[\begin{split} \begin{align*} s_3 &= s_2 - g(s_2) \left( \frac{s_2 - s_1}{g(s_2) - g(s_1)} \right) \\ &= 1.5 - (-0.3784) \left( \frac{1.5 - 2}{-0.3784 - (-1.1712)} \right) = 1.2614. \end{align*} \end{split}\]

The solution using a guess value of \(s_3 = 1.2614\) are tabulated below.

\(t\)	\(y_1\)	\(y_2\)
0.00	0.0000	1.2614
0.20	0.2523	1.5136
0.40	0.5550	1.8668
0.60	0.9284	2.3512
0.80	1.3986	3.0071
1.00	2.0000	3.8882

Here the solution at the upper boundary is \(y(t_{\max}) = 2\) which agrees to the required boundary condition \(y(1) = 2\) correct to 4 decimal places.

Fig. 5.2 Solution of the BVP \(y'' - y' - y = 0\), \(y(0) = 0\), \(y(1) = 2\) using the shooting method.

5.2.2. Code

The code below defines a function called shooting_method() that calculates the solution to a boundary value problem using the shooting method.

Python

def shooting_method(f, tspan, bvals, h, solver, tol=1e-8, maxiter=100):
    s0 = 1
    s1 = s0 + 1
    [t, y] = solver(f, tspan, [bvals[0], s0], h)
    g0 = bvals[1] - y[-1,0]

    for _ in range(maxiter):
        [t, y] = solver(f, tspan, [bvals[0], s1], h)
        g1 = bvals[1] - y[-1,0]
        s2 = s1 - g1 * (s1 - s0) / (g1 - g0)

        if abs(s2 - s1) < tol:
            return [t, y]

        s0 = s1
        g0 = g1
        s1 = s2

MATLAB

function [t, y] = shooting_method(f, tspan, bvals, h, solver)

s0 = 1;
s1 = s0 + 1;
[t, y] = solver(f, tspan, [bvals(1), s0], h);
g0 = bvals(2) - y(end, 1);

for i = 1 : 10
    [t, y] = solver(f, tspan, [bvals(1), s1], h);
    g1 = bvals(2) - y(end, 1);
    s2 = s1 - g1 * (s1 - s0) / (g1 - g0);

    if abs(s2 - s1) < 1e-8
        return
    end

    s0 = s1;
    g0 = g1;
    s1 = s2;
end

end

The inputs to the function are

• f - the name of the ODE function

• tspan - an array of two values defining the boundaries of the \(t\) domain

• bvals - an array of two values defining the upper and lower boundary values

• h - the step length used in the solver function

• solver - the name of the solver function, e.g., euler, rk4 etc. (this needs to be defined elsewhere)

5.2.3. A note about accuracy

In Example 5.3 we iterated the Secant method until convergence and got a very accurate solution for the value of \(y(b)\). We must be careful not to forget that this solution was obtained using the Euler method which being only first-order so expect it to be relatively inaccurate for the other values in the domain.

The exact solution to the boundary value problem used here is

\[ \begin{align*} y = \frac{2e^{(1 - \sqrt{5})(t - 1)/2} (e^{\sqrt{5}t} - 1)}{e^{\sqrt{5}} - 1}. \end{align*} \]

and the solutions using the Euler method and the RK4 method have been tabulated below and plotted in Fig. 5.3.

\(t\)	Euler	RK4	Exact	Euler error	RK4 error
0.00	0.000000	0.000000	0.000000	0.00e+00	0.00e+00
0.20	0.252270	0.221310	0.221296	3.10e-02	1.41e-05
0.40	0.554995	0.501444	0.501419	5.36e-02	2.50e-05
0.60	0.928355	0.865870	0.865840	6.25e-02	3.00e-05
0.80	1.398587	1.349436	1.349412	4.92e-02	2.42e-05
1.00	2.000000	2.000000	2.000000	0.00e+00	4.44e-16

Fig. 5.3 Solutions to the boundary value problem \(y'' - y' - y = 0\), \(t \in [0,1]\), \(y(0) = 0\), \(y(1) = 2\) using the Euler and RK4 methods with \(h=0.2\).

So despite the Secant method giving a guess value that gives an accurate solution at the upper boundary, the use of the Euler method does not give an accurate solution across the domain. The solution using the second-order Runge-Kutta method in comparison gives much more accurate solutions over the domain.