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 first-order ODEs
, (1)
,
The vector contains the unknown functions of the independent variable . The unknown functions are ordered so that the first (0< < ) components of have first-kind boundary conditions at . The remaining := components of the solution have first-kind boundary conditions specified at a second point, . The Maple procedure shoot supports nonlinear boundary conditions at . In particular, can be replaced with , , where each is a differentiable function of n variables. The examples discussed in this paper all use := . (Note that if , then (1) is an initial value problem.)
>
The shooting method seeks to identify a vector of parameters in so that the solution, denoted by , to the initial value problem
,
,
agrees with the solution to (1). Note that (2) is simply (1) with the boundary conditions at replaced with unknown initial conditions at . To determine the correct initial values, consider the ``objective function'' with components
:= , .
Then [1, p. 476], (1) is solvable if and only if there exists in so that .
>
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 . 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 in , define a sequence of initial conditions { } by
:=
for all >= 0. To implement this, note that the vector is directly available from the solution of (2), but the Jacobian matrix requires the values of for all . These values can be obtained by solving the IVPs in (2) together with the sensitivity equations ([2, p. 226], [3, pp. 54--58]):
, ,
with corresponding initial conditions
, for all
,
, for all
>
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 . Even if (2) does have a solution on [a,b] the problem may be stiff. In such a case the solution at may be so inaccurate as to make the results of the Newton-Raphson step meaningless. When the accuracy of the solution at 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 increases, it seems appropriate to consider partitioning the interval into 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.
>