A Maple Implementation of the Simple Shooting Method

The basic idea of the shooting method for two-point boundary value problems is to reformulate the problem as a nonlinear parameter estimation problem. The new problem requires the solution of a related IVP with initial conditions chosen to approximate the boundary conditions at the other endpoint. If these boundary conditions are not satisfied to the desired accuracy, the process is repeated with a new set of initial conditions until the desired accuracy is achieved or an iteration limit is reached.

>

To be more specific, consider the two-point BVP for a coupled system of n first-order ODEs

diff(y(t),t) = f(t,y(t))

y[i](a) = alpha[i] , i = 1 .. m[1] (1)

y[m[1]+j](b) = beta[j] , j = 1 .. m[2]

The vector y contains the n unknown functions of the independent variable t . The unknown functions are ordered so that the first m[1] (0< m[1] < n ) components of y have first-kind boundary conditions at t = a . The remaining m[2] := n-m[1] components of the solution have first-kind boundary conditions specified at a second point, t = b . The Maple procedure shoot supports nonlinear boundary conditions at t = b . In particular, `(1)`[3] can be replaced with r[j](y(b)) = 0 , j = 1 .. m[2] , where each r[j] is a differentiable function of n variables. The examples discussed in this paper all use r[j](y(b)) := y[m[1]+j](b)-beta[j] . (Note that if m[2] = 0 , then (1) is an initial value problem.)

>

The shooting method seeks to identify a vector of parameters s in R^m[2] so that the solution, denoted by y(t,s) , to the initial value problem

diff(y,t) = f(t,y(t,s))

y[i](a,s) = alpha[i] , i = 1 .. m[1]

y[m[1]+j](a,s) = s[j] , j = 1 .. m[2]

agrees with the solution to (1). Note that (2) is simply (1) with the boundary conditions at t = b replaced with unknown initial conditions at t = a . To determine the correct initial values, consider the ``objective function'' F with components

F[j](s) := y[m[1]+j](b,s)-beta[j] , j = 1 .. m[2] .

Then [1, p. 476], (1) is solvable if and only if there exists s in R^m[2] so that F(s) = 0 .

>

The success of this process depends primarily on the iterative procedure used to construct a sequence of parameter vectors that converges to a zero of F . While any numerical root-finding algorithm could be employed for this step, one step of the Newton-Raphson method is most commonly used. That is, given an initial guess s^0 in R^m[2] , define a sequence of initial conditions { s^k } by

s^(k+1) := s^k-JF(s^k)^(-1)*F(s^k)

for all k >= 0. To implement this, note that the vector F(s^k) is directly available from the solution of (2), but the Jacobian matrix JF(s^k) requires the values of diff(y[m[1]+i](b,s^k),s[j]) for all i, j = 1 .. m[2] . These values can be obtained by solving the n IVPs in (2) together with the n*m[2] sensitivity equations ([2, p. 226], [3, pp. 54--58]):

diff(diff(y[i](t),s[j]),t) = diff(f,y[i])*diff(y[i]... , i = 1 .. n , j = 1 .. m[2]

with corresponding initial conditions

diff(y[i],s[j]) = 0 , for all i = 1 .. m[1] , j = 1 .. m[2]
diff(y[m[1]+j],s[j]) = delta[i*j] , for all i, j = 1 .. m[2]

>

The above algorithm is known as the simple , or single , shooting method . While this method is effective for many problems, three specific problems should be mentioned [4]. Assume (1) has a unique solution. There is no guarantee that the initial value problem (2) has a solution on the interval [a,b] for all s in R^m[2] . Even if (2) does have a solution on [a,b] the problem may be stiff. In such a case the solution at t = b may be so inaccurate as to make the results of the Newton-Raphson step meaningless. When the accuracy of the solution at t = b is known with sufficient accuracy, the local convergence of the Newton-Raphson step may prevent the iterations from converging to a solution of the original boundary value problem (1). This difficulty can be addressed by replacing the Newton-Raphson step with another iterative solver with improved convergence properties ( e.g. , modified Newton's method and Broyden's method). The multiple , or parallel , shooting method addresses the other difficulties. As these problems are generally more pronounced as b-a increases, it seems appropriate to consider partitioning the interval into N subintervals. Then, using compatibility conditions between the subintervals, a well-posed IVP is obtained on each subinterval.

>

The combination of Maple's facilities for symbolic manipulation and numerical solution of IVPs provide an excellent environment for the translation of the simple shooting method into the Maple programming language. The syntax for this procedure, called shoot , closely follows the syntax of dsolve . In particular, the data structures returned by shoot are identical to the ones returned by dsolve/numeric . Thus, all the techniques and tools for manipulating Maple's numerical solutions of differential equations, e.g. plots[odeplot] , can be used to interpret the results from shoot . (The source code for shoot , including on-line help documentation, is available upon request from the first author.)

>

Other numerical methods for the solution of two-point boundary value problems have also been developed on other platforms. One of the most well-known is the fortran subroutine COLSYS (available from Netlib). The numerical method used in this routine is collocation of B-splines at Gaussian points [5]. Finite difference methods can be easily implemented using Matlab 's fsolve command.

>

The three examples discussed in this paper provide samples of how simple shooting, in particular shoot , can be applied to the analysis of boundary value problems that are encountered in chemical engineering.

>