5.2. The shooting method

Consider the two-point boundary value problem

\[ y'' = f(t,y), \qquad t \in [t_{\min}, t_{\max}], \qquad y(t_{\min}) = \alpha ,\qquad y(t_{\max}) = \beta. \]

Since the ODE solvers we use can only be applied to a first-order ODE we need to rewrite the second-order ODE as a system of two first-order ODEs

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

Since we have let \(y_1 = y\) we know that \(y_1(a) = \alpha\) from the definition of the boundary value problem but we do not know the value of \(y_2(a)\). To overcome this problem we simply guess this value and compute the solutions to the initial value problem and compare the solution of \(y_1(b)\) to \(\beta\) and then adjust the guess value accordingly. This method is called the shooting method because someone shooting at a target will adjust their next shot based where their previous shot landed.

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

\[ y'' - y' - y = 0, \qquad y(0) = 0, \qquad y(1) = 2. \]

Solution (click to show)

Rewriting this as a system of first-order ODEs gives

\[\begin{split} \begin{align*} y_1' &= y_2, & y_1 (0) &= 1, &&\\ 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 \(y'(0) = 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(5) = 1.585600\) which is less than the boundary value of \(y(1) = 2\). Lets try again but with an increased guess value of \(y'(0) = 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(5) = 3.171200\) which is greater than 2. Lets try again but with a decreased guess value of \(y'(0) = 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.

5.2.1. Improving the guess value using the secant method

In Example 5.2 we saw that by adjusting the guess value of \(y'(0)\) we obtained solutions at the upper boundary that were closer to the target value at the upper boundary \(\beta\). We could continue in this way by increasing or decreasing our guess value halfway towards the previous guess value depending on whether the solution at the upper boundary was higher or lower than the target (this method is known as the bisection method), however a more efficient approach is to use the Secant method. The Secant method is a root finding algorithm that calculates the value of \(s\) where \(g(s)=0\) for some function \(g\)

(5.3)\[s_{i+1} = s_i - g(s_i)\frac{s_i - s_{i-1}}{g(s_i) - g(s_{i-1})}.\]

This expression is iterated until two successive values of \(s\) are less than some small number, i.e., \(|s_i - s_{i-1}| < tol\). Since we want the solution of \(y'(a)\) to be equal to the upper boundary value \(\beta\) we can define the function

\[ g(s) = \beta - y_n. \]

Using the Secant method to find the root \(g(s) = 0\) will give the value of \(y'(0)\).

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 (click to show)

From Example 5.2 we have \(s_1 = 2\) and \(s_2 = 1.5\) which resulted in solutions at the upper boundary of \(y_5(s = 1) = 1.5856\), \(y_5(s = 2) = 3.1712\) and \(y_5(s = 1.5) = 2.3784\). 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.3) 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 = 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_5 = 2\) which is equal to \(y(5)=2\) correct to (at least) 4 decimal places.

Algorithm 5.1 (Shooting method for solving boundary value problems)

Inputs A system of two first-order ODEs \(\vec{y}' = f(t, y_1, y_2)\), a domain \(t\in [t_{\min}, t_{\max}]\), the boundary values \(y(t_{\min}) = \alpha\) and \(y(t_{\max}) = \beta\) and an accuracy tolerance \(tol\).

Outputs \((t_0, t_1, \ldots, t_n)\) and \((\vec{y}_0, \vec{y}_1, \ldots, \vec{y}_n)\)

• \(s_0 \gets 1\)

• \(s_1 \gets 2\)

• \(g(s_0) \gets 1\)

• For \(i = 1, \ldots, 10\)

  • Solve ODE using initial values \(\vec{y}_0 = (\alpha, s_i)\)

  • \(g(s_{i}) \gets \beta - y_1(t_{\max})\)

  • \(s_{i+1} \gets s_i - g(s_i) \dfrac{s_i - s_{i-1}}{g(s_i) - g(s_{i-1})}\)

  • If \(|s_{i+1} - s_{i}| < tol\)

    • Break

• Return \((t_0, t_1, \ldots, t_n)\) and \((\vec{y}_0, \vec{y}_1, \ldots, \vec{y}_n)\)

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-6):
    s, s_old, g_old = 1, 2, 1
    for _ in range(10):
        t, y = solveIVP(f, tspan, [bvals[0], s], h, solver)
        g = bvals[1] - y[-1,0]
        snew = s - g * (s - s_old) / (g - g_old)
        s, s_old, g_old = snew, s, g
        if abs(s - s_old) < tol:
            break
        
    return t, y

MATLAB

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

s = 1;
s_old = 2;
g_old = 1;
for i = 1 : 10
    [t, y] = solveIVP(f, tspan, [bvals(1), s], h, solver);
    g = bvals(2) - y(end, 1);
    snew = s - g * (s - s_old) / (g - g_old);
    s_old = s;
    s = snew;
    g_old = g;
    if abs(s - s_old) < tol
        break
    end
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)

• tol - the convergence tolerance for the Secant method

The variables s, s_old, g and g_old are the current and previous values of the guess value \(s\) and function \(g(s)\) which are intialised to ensure that at least two iterations of the Secant method are performed. The function calculates the solution to the initial value problem with the initial solution [bvals[0], s] and uses the last solution for \(y_1\) to calculate the next guess value using the Secant method. The iterations cease when the difference between two successive guess values are less than tol.

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 second-order Runge-Kutta method have been tabulated in Table 5.1 and plotted in Fig. 5.1.

Table 5.1 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\).

\(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	0.00e+00

Fig. 5.1 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.