Methods

Chasing

The method of chasing came from a manuscript of Gel`fand 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.

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 homogenous 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 , we have the initial value problem for ,

Given where we want to have solutions to all of the boundary value problems, Mathematica just uses NDSolve to solve the auxiliary problems for , , ... by integrating them to . The results are then combined into the matrix of Equation (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 Gel'fand 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 = (1, ,,,,,,) and = . Using (5) and replacing , and thinking of in terms of the other components of we get the nonlinear equation

where is with the column removed and is the 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 error 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 Mathematica is to use the linear method.

This solves a two-point boundary value problem for a second-order equation.

In[1]:=

Out[1]=

This shows a plot of the solution.

In[2]:=

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 t.)

In[3]:=

Out[4]=

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

This checks to see if the boundary conditions are "satisfied".

In[5]:=

Out[5]=

They are typically only be 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.

This checks the residual error at each of the boundary conditions.

In[6]:=

Out[6]=

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.

This is a boundary value problem that has no solution.

In[7]:=

Out[7]=

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.

Here is a boundary value problem that does not have a unique solution. Its general solution is shown as computed symbolically with DSolve.

In[8]:=

Out[8]=

NDSolve issues a warning message because the matrix to solve for the initial conditions is singular, but has a solution.

In[9]:=

Out[9]=

You can identify which solution it found by fitting it to the interpolating points. This makes a plot of the error relative to the actual best fit solution.

In[10]:=

Typically the default values Mathematica 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.

The method ChasingType->NonlinearChasing itself has two options.

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

Here is a boundary value problem with a simple solution computed symbolically using DSolve.

In[15]:=

Out[16]=

This shows the error in the solution computed by NDSolve.

In[17]:=

Using extra precision to solve for the initial conditions reduces the error substantially.

In[19]:=

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.

This uses extra precision to compute the initial conditions along with more stringent settings for the AccuracyGoal and PrecisionGoal options.

In[21]:=