"Shooting" Method

The shooting method works by considering the boundary conditions as a multivariate function of initial conditions at some point, reducing the boundary value problem to finding the initial conditions that give a root. The advantage of the shooting method is that it takes advantage of the speed and adaptivity of methods for initial value problems. The disadvantage of the method is that it is not as robust as finite difference or collocation methods: some initial value problems with growing modes are inherently unstable even though the BVP itself may be quite well posed and stable.

Consider the BVP system

The shooting method looks for initial conditions so that . Since you are varying the initial conditions, it makes sense to think of as a function of them, so shooting can be thought of as finding such that

After setting up the function for , the problem is effectively passed to FindRoot to find the initial conditions giving the root. The default method is to use Newton's method, which involves computing the Jacobian. While the Jacobian can be computed using finite differences, the sensitivity of solutions of an initial value problem (IVP) to its initial conditions may be too much to get reasonably accurate derivative values, so it is advantageous to compute the Jacobian as a solution to ODEs.

Linearization and Newton's Method

Linear problems can be described by

where is a matrix and is a vector both possibly depending on , is a constant vector, and are constant matrices.

Let . Then, differentiating both the IVP and boundary conditions with respect to gives

Since is linear, when thought of as a function of , you have , so the value of for which satisfies

for any particular initial condition .

For nonlinear problems, let be the Jacobian for the nonlinear ODE system, and let be the Jacobian of the ^th boundary condition. Then computation of for the linearized system gives the Jacobian for the nonlinear system for a particular initial condition, leading to a Newton iteration,

"StartingInitialConditions"

For boundary value problems, there is no guarantee of uniqueness as there is in the initial value problem case. "Shooting" will find only one solution. Just as you can affect the particular solution FindRoot gets for a system of nonlinear algebraic equations by changing the starting values, you can change the solution that "Shooting" finds by giving different initial conditions to start the iterations from.

"StartingInitialConditions" is an option of the "Shooting" method that allows you to specify the values and position of the initial conditions to start the shooting process from.

The shooting method by default starts with zero initial conditions so that if there is a zero solution, it will be returned.

This computes a very simple solution to the boundary value problem with :

By default, "Shooting" starts from the left side of the interval and shoots forward in time. There are cases where it is advantageous to go backward, or even from a point somewhere in the middle of the interval.

Consider the linear boundary value problem

that has a solution

For moderate values of , the initial value problem starting at becomes unstable because of the growing and terms. Similarly, starting at , instability arises from the term, though this is not as large as the term in the forward direction. Beyond some value of , shooting will not be able to get a good solution because the sensitivity in either direction will be too great. However, up to that point, choosing a point in the interval that balances the growth in the two directions will give the best solution.

This solves the system with shooting from the default :

Shooting from , the "Shooting" method gives warning messages about an ill-conditioned matrix and that the boundary conditions are not satisfied as well as they should be. This is because a small error at is amplified by Since the reciprocal of this is of the same order of magnitude as the local truncation error, visible errors such as those seen in the plot are not surprising. In the reverse direction, the magnification will be much less: , so the solution should be much better.

This computes the solution using shooting from :

A good point to choose is one that will balance the sensitivity in each direction, which is about at . With this, the error with will still be under reasonable control.

This computes the solution for shooting from

Option Summary

"Shooting" method options.

"Chasing" Method

The method of chasing came from a manuscript of Gelfand and Lokutsiyevskii first published in English in [BZ65] and further described in [Na79]. The idea is to establish a set of auxiliary problems that can be solved to find initial conditions at one of the boundaries. Once the initial conditions are determined, the usual methods for solving initial value problems can be applied. The chasing method is, in effect, a shooting method that uses the linearity of the problem to good advantage. As such, the chasing method can be used for linear differential equations.

Consider the linear ODE

where , is the coefficient matrix, and is the inhomogeneous coefficient vector, with linear boundary conditions

where is a coefficient vector. From this, construct the augmented homogeneous system

where

The chasing method amounts to finding a vector function such that and . Once the function is known, if there is a full set of boundary conditions, solving

can be used to determine initial conditions that can be used with the usual initial value problem solvers. Note that the solution to system (3) is nontrivial because the first component of is always 1.

Thus, solving the boundary value problem is reduced to solving the auxiliary problems for the . Differentiating the equation for gives

substituting the differential equation,

and transposing

Since this should hold for all solutions , you have the initial value problem for ,

Given where you want to have solutions to all of the boundary value problems, the Wolfram Language just uses NDSolve to solve the auxiliary problems for by integrating them to . The results are then combined into the matrix of (3) that is solved for to obtain the initial value problem that NDSolve integrates to give the returned solution.

This variant of the method is further described in and used by the MathSource package [R98], which also allows you to solve linear eigenvalue problems.

There is an alternative, nonlinear way to set up the auxiliary problems that is closer to the original method proposed by Gelfand and Lokutsiyevskii. Assume that the boundary conditions are linearly independent (if not, then the problem is insufficiently specified). Then in each , there is at least one nonzero component. Without loss of generality, assume that . Now solve for in terms of the other components of , , where and . Using (5) and replacing , and thinking of in terms of the other components of , you get the nonlinear equation

where is with the ^thcolumn removed, and is the ^th column of A. The nonlinear method can be more numerically stable than the linear method, but it has the disadvantage that integration along the real line may lead to singularities. This problem can be eliminated by integrating on a contour in the complex plane. However, the integration in the complex plane typically has more numerical errors than a simple integration along the real line, so in practice, the nonlinear method does not typically give results better than the linear method. For this reason, and because it is also generally faster, the default for the Wolfram Language is to use the linear method.

The solver can solve multipoint boundary value problems of linear systems of equations. (Note that each boundary equation must be at one specific value of .)

In general, you cannot expect the boundary value equations to be satisfied to the close tolerance of Equal.

They are typically only satisfied at most tolerances that come from the AccuracyGoal and PrecisionGoal options of NDSolve. Usually, the actual accuracy and precision is less because these are used for local, not global, error control.

When you give NDSolve a problem that has no solution, numerical error may make it appear to be a solvable problem. Typically, NDSolve will issue a warning message.

In this case, it is not able to integrate over the entire interval because of nonexistence.

Another situation in which the equations can be ill-conditioned is when the boundary conditions do not give a unique solution.

Typically, the default values the Wolfram Language uses work fine, but you can control the chasing method by giving NDSolve the option Method->{"Chasing",chasing options}. The possible chasing options are shown in the following table.

Options for the "Chasing" method of NDSolve.

The method "ChasingType"->"NonlinearChasing" itself has two options.

Options for the "NonlinearChasing" option of the "Chasing" method.

These options, especially "ExtraPrecision", can be useful in cases where there is a strong sensitivity to computed initial conditions.

Increasing the extra precision beyond this really will not help because a significant part of the error results from computing the solution once the initial conditions are found. To reduce this, you need to give more stringent AccuracyGoal and PrecisionGoal options to NDSolve.

Boundary Value Problems with Parameters

In many of the applications where boundary value problems arise, there may be undetermined parameters, such as eigenvalues, in the problem itself that may be a part of the desired solution. By introducing the parameters as dependent variables, the problem can often be written as a boundary value problem in standard form.

For example, the flow in a channel can be modeled by

where (the Reynolds number) is given, but is to be determined.

To find the solution and the value of , just add the equation .

This solves the flow problem with for and , plots the solution and returns the value of :