Numerical Solution of Boundary Value Problems (BVP)

"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 X((Subscript[t, 0]))==c so that G==0. Since you are varying the initial conditions, it makes sense to think of X==(Subscript[X, c]) as a function of them, so shooting can be thought of as finding c such that

After setting up the function for G, the problem is effectively passed to (RefLink[FindRoot, paclet:ref/FindRoot]) to find the initial conditions c 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 J(t) is a matrix and (Subscript[F, 0])(t) is a vector both possibly depending on t, (Subscript[B, 0]) is a constant vector, and (Subscript[B, 1]),(Subscript[B, 2]),…,(Subscript[B, n]) are constant matrices.

Let Y==((∂(Subscript[X, c])(t))/(∂c)). Then, differentiating both the IVP and boundary conditions with respect to c gives

Since G is linear, when thought of as a function of c, you have G(c)==G((Subscript[c, 0]))+(((∂G)/(∂c)))(c-(Subscript[c, 0])), so the value of c for which G(c)==0 satisfies

for any particular initial condition (Subscript[c, 0]).

For nonlinear problems, let J(t) be the Jacobian for the nonlinear ODE system, and let (Subscript[B, i]) be the Jacobian of the i(^th) boundary condition. Then computation of ((∂G)/(∂c)) 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 (RefLink[FindRoot, paclet:ref/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 x''+sin(x)==0 with  x(0)==x(10)==0:

sols = Map[First[NDSolve[{x''[t] + Sin[x[t]] == 0, x[0] == x[10] == 0},x,t, Method->{"Shooting", "StartingInitialConditions"->{x[0] == 0, x'[0] == #}}]]&,
{1.5, 1.75, 2}];
Plot[Evaluate[x[t] /. sols], {t,0,10}, PlotStyle->{Black, Blue, Green}]

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 t==0 becomes unstable because of the growing (e^(λ(t-1))) and (e^(2λ(t-1))) terms. Similarly, starting at t==1, instability arises from the (e^(-λt)) 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 gives the equation, boundary conditions, and exact solution as Wolfram Language input:

eqn = x'''[t] - 2λ x''[t] - (λ^2)x'[t] + 2(λ^3)x[t] == ((λ^2)+(π^2)) (2 λ Cos[π t]+π Sin[π t]);
bcs = {x[0] ==  1 + ((1+(E^(-2 λ))+(E^(-λ)))/(2+(E^(-λ)))), x[1] == 0, x'[1] == ((3 λ-(E^(-λ)) λ)/(2+(E^(-λ))))};
xsol[t_] =  (((E^(λ (t-1)))+(E^(2 λ (t-1)))+(E^(-λ t)))/(2 + (E^(-λ))))+Cos[π t];

This solves the system with λ==10 shooting from the default t==0:

Block[{λ = 10}, 
sol = First[NDSolve[{eqn, bcs}, x, t]];
Plot[{xsol[t],x[t] /. sol}, {t,0,1}]]

Shooting from t==0, 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 t==0 is amplified by (e^20)≃4×(10^8). 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: (e^10)≃2×(10^4), so the solution should be much better.

This computes the solution using shooting from t==1:

Block[{λ = 10}, 
sol = First[NDSolve[{eqn, bcs}, x, t, 
Method->{"Shooting", "StartingInitialConditions"->{x[1] == 0, x'[1] == ((3 λ-(E^(-λ)) λ)/(2+(E^(-λ)))), x''[1] == 0}}]];
Plot[{xsol[t],x[t] /. sol}, {t,0,1}]]

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

This computes the solution for λ==15 shooting from t==2/3.

Block[{λ = 15}, 
sol = First[NDSolve[{eqn, bcs}, x, t, 
Method->{"Shooting", "StartingInitialConditions"->{x[2/3] == 0, x'[2/3] == 0, x''[2/3] == 0}}]];
Plot[{xsol[t],x[t] /. sol}, {t,0,1}]]

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.

Consider the linear ODE

X'(t)==A(t)X(t)+(Subscript[A, 0])(t)

where X(t)==((Subscript[x, 1])(t),(Subscript[x, 2])(t),…,(Subscript[x, n])(t)), A(t) is the coefficient matrix, and (Subscript[A, 0])(t) is the inhomogeneous coefficient vector, with n linear boundary conditions

(Subscript[B, i]).X((Subscript[t, i]))==(Subscript[b, i0]),i==1,2,⋯,n

where (Subscript[B, i])==((Subscript[b, i1]),(Subscript[b, i2]),⋯,(Subscript[b, in])) is a coefficient vector. From this, construct the augmented homogeneous system

(OverBar[X])'(t)==(OverBar[A])(t)(OverBar[X])(t),         (Subscript[(OverBar[B]), i]).(OverBar[X])((Subscript[t, i]))==0

where

The chasing method amounts to finding a vector function (Subscript[Φ, i])(t) such that (Subscript[Φ, i])((Subscript[t, i]))==(OverBar[(Subscript[B, i])]) and (Subscript[Φ, i])(t)(OverBar[X])(t)==0. Once the function (Subscript[Φ, i])(t) is known, if there is a full set of boundary conditions, solving

(({{(Subscript[Φ, 1])((Subscript[t, 0]))}, {(Subscript[Φ, 2])((Subscript[t, 0]))}, {⋮}, {(Subscript[Φ, n])((Subscript[t, 0]))}}))(OverBar[X])((Subscript[t, 0]))==0

can be used to determine initial conditions ((Subscript[x, 1])((Subscript[t, 0])),(Subscript[x, 2])((Subscript[t, 0])),⋯,(Subscript[x, n])((Subscript[t, 0]))) 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 X is always 1.

Thus, solving the boundary value problem is reduced to solving the auxiliary problems for the (Subscript[Φ, i])(t). Differentiating the equation for (Subscript[Φ, i])(t) gives

substituting the differential equation,

and transposing

Since this should hold for all solutions (OverBar[X]), you have the initial value problem for (Subscript[Φ, i]),

(Subscript[Φ, i])'(t)+((OverBar[A])^T)(t)(Subscript[Φ, i])(t)=0  with initial condition (Subscript[Φ, i])((Subscript[t, i]))=(Subscript[B, i])

Given (Subscript[t, 0]) where you want to have solutions to all of the boundary value problems, the Wolfram Language just uses (RefLink[NDSolve, paclet:ref/NDSolve]) to solve the auxiliary problems for (Subscript[Φ, 1]),(Subscript[Φ, 2]),…,(Subscript[Φ, m]) by integrating them to (Subscript[t, 0]). The results are then combined into the matrix of (3) that is solved for X((Subscript[t, 0])) to obtain the initial value problem that (RefLink[NDSolve, paclet:ref/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 (Subscript[B, i]), there is at least one nonzero component. Without loss of generality, assume that (Subscript[b, ij])≠0. Now solve for (Subscript[Φ, ij]) in terms of the other components of (Subscript[Φ, i]), (Subscript[Φ, ij])==(OverTilde[(Subscript[B, i])]).(OverTilde[(Subscript[Φ, i])]), where (OverTilde[(Subscript[Φ, i])])==(1,(Subscript[Φ, i1]),⋯,(Subscript[Φ, ij-1]),⋯,(Subscript[Φ, ij+1]),⋯,(Subscript[Φ, in])) and (OverTilde[(Subscript[B, i])])==((Subscript[b, i0, ])(Subscript[b, i1]),⋯,(Subscript[b, ij-1]),⋯,(Subscript[b, ij+1]),⋯,(Subscript[b, in]))/-(Subscript[b, ij]). Using (5) and replacing (Subscript[Φ, ij]), and thinking of (Subscript[x, n])(t) in terms of the other components of x(t), you get the nonlinear equation

where (OverTilde[A]) is A with the j(^th)column removed, and (Subscript[A, j]) is the j(^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.

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

nsol1 = NDSolve[{y''[t] + y[t]/4 == 8, y[0] == 0, y[10] == 0},y,{t,0,10}]

This shows a plot of the solution:

Plot[First[y[t] /. nsol1], {t,0,10}]

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

bconds = {
x[0] + x'[0] + y[0] + y'[0] == 1,
x[1] + 2 x'[1] + 3 y[1] + 4y'[1] == 5,
y[2] + 2 y'[2] == 4,
x[3] - x'[3] == 7};
nsol2 = NDSolve[{
x''[t] + x[t] + y[t] == t,y''[t] + y[t] == Cos[t],
bconds},
{x,y},
{t,0,4}]

In general, you cannot expect the boundary value equations to be satisfied to the close tolerance of (RefLink[Equal, paclet:ref/Equal]).

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

bconds /. First[nsol2]

They are typically only satisfied at most tolerances that come from the (RefLink[AccuracyGoal, paclet:ref/AccuracyGoal]) and (RefLink[PrecisionGoal, paclet:ref/PrecisionGoal]) options of (RefLink[NDSolve, paclet:ref/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:

Apply[Subtract, bconds, 1] /. First[nsol2]

When you give (RefLink[NDSolve, paclet:ref/NDSolve]) a problem that has no solution, numerical error may make it appear to be a solvable problem. Typically, (RefLink[NDSolve, paclet:ref/NDSolve]) will issue a warning message.

This is a boundary value problem that has no solution:

NDSolve[{x''[t] + x[t] == 0, x[0] == 1, x[Pi] == 0},x,{t,0,Pi},Method->"Chasing"]

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 (RefLink[DSolve, paclet:ref/DSolve]):

dsol = First[x /. DSolve[{x''[t] + x[t] == t, x'[0] == 1, x[Pi/2] == Pi/2},x,t]]

(RefLink[NDSolve, paclet:ref/NDSolve]) issues a warning message because the matrix to solve for the initial conditions is singular, but has a solution:

onesol = First[x /. NDSolve[{x''[t] + x[t] == t, x'[0] == 1, x[Pi/2] == Pi/2},x,{t, 0, Pi/2},Method->"Chasing"]]

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:

ip = onesol@"Coordinates"[1];
points = Transpose[{ip, onesol[ip]}];
model = dsol[t] /. (C[1])->α;
fit = FindFit[points, model,α, t];
ListPlot[Transpose[{ip, onesol[ip] - ((model /. fit) /. t->ip)}]]

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

Options for the "Chasing" method of (RefLink[NDSolve, paclet:ref/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.

Here is a boundary value problem with a simple solution computed symbolically using (RefLink[DSolve, paclet:ref/DSolve]):

bvp = {x''[t] + 1000 x[t] == 0, x[0] == 0, x[1] == 1};
dsol = First[x /. DSolve[bvp,x,t]]

This shows the error in the solution computed using the chasing method in (RefLink[NDSolve, paclet:ref/NDSolve]):

sol = First[x /. NDSolve[{x''[t] + 1000 x[t] == 0, x[0] == 0, x[1] == 1},x,{t, 0, 1},Method->"Chasing"]];
Plot[sol[t] - dsol[t],{t,0,1}]

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

sol = First[x /. NDSolve[{x''[t] + 1000 x[t] == 0, x[0] == 0, x[1] == 1},x,{t, 0, 1}, Method->{"Chasing", "ExtraPrecision"->10}]];
Plot[sol[t] - dsol[t],{t,0,1}]

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 (RefLink[AccuracyGoal, paclet:ref/AccuracyGoal]) and (RefLink[PrecisionGoal, paclet:ref/PrecisionGoal]) options to (RefLink[NDSolve, paclet:ref/NDSolve]).

This uses extra precision to compute the initial conditions along with more stringent settings for the (RefLink[AccuracyGoal, paclet:ref/AccuracyGoal]) and (RefLink[PrecisionGoal, paclet:ref/PrecisionGoal]) options:

sol = First[x /. NDSolve[{x''[t] + 1000 x[t] == 0, x[0] == 0, x[1] == 1},x,{t, 0, 1}, Method->{"Chasing", "ExtraPrecision"->10}, AccuracyGoal->10, PrecisionGoal->10]];
Plot[sol[t] - dsol[t],{t,0,1}]

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 R (the Reynolds number) is given, but a is to be determined.

To find the solution f  and the value of a, just add the equation a'==0.

This solves the flow problem with R==1 for f and a, plots the solution f and returns the value of a:

Block[{R = 1},sol = NDSolve[{f'''[t] - R((f'[t])^2 - f[t] f''[t]) + R a[t] == 0, a'[t] == 0, f[0] == f'[0]== f'[1] == 0, f[1] == 1},{f,a},t];
Column[{Plot[f[t] /. First[sol],{t,0,1}],
a[0] /. First[sol]}]]