Solving Boundary Value Problems

In a boundary value problem (BVP), the goal is to find a solution to an ordinary differential equation (ODE) that also satisfies certain specified boundary conditions. The boundary conditions specify a relationship between the values of the solution at two or more locations in the interval of integration. In the simplest case, the boundary conditions apply at the beginning and end (or boundaries) of the interval.

The MATLAB^® BVP solvers bvp4c and bvp5c are designed to handle systems of ODEs of the form

$$y\text{'}=f\left(x,y\right)$$

where:

• x is the independent variable

• y is the dependent variable

• y′ represents the derivative of y with respect to x, also written as dy/dx

Boundary Conditions

In the simplest case of a two-point BVP, the solution to the ODE is sought on an interval [a, b], and must satisfy the boundary conditions

$$g\left(y\left(a\right),y\left(b\right)\right)=0.$$

To specify the boundary conditions for a given BVP, you must:

• Write a function of the form res = bcfun(ya,yb), or use the form res = bcfun(ya,yb,p) if there are unknown parameters involved. You supply this function to the solver as the second input argument. The function returns res, which is the residual value of the solution at the boundary point. For example, if y(a) = 1 and y(b) = 0, then the boundary condition function is

  function res = bcfun(ya,yb)
  res = [ya(1)-1
         yb(1)];
  end

• In the initial guess for the solution, the first and last points in the mesh specify the points at which the boundary conditions are enforced. For the above boundary conditions, you can specify bvpinit(linspace(a,b,5),yinit) to enforce the boundary conditions at a and b.

The BVP solvers in MATLAB also can accommodate other types of problems that have:

• Unknown parameters p

• Singularities in the solutions

• Multipoint conditions (internal boundaries that separate the integration interval into several regions)

In the case of multipoint boundary conditions, the boundary conditions apply at more than two points in the interval of integration. For example, the solution might be required to be zero at the beginning, middle, and end of the interval. See bvpinit for details on how to specify multiple boundary conditions.

Initial Guess of Solution

Unlike initial value problems, a boundary value problem can have:

• No solution

• A finite number of solutions

• Infinitely many solutions

An important part of the process of solving a BVP is providing a guess for the required solution. The quality of this guess can be critical for the solver performance and even for a successful computation.

Use the bvpinit function to create a structure for the initial guess of the solution. The solvers bvp4c and bvp5c accept this structure as the third input argument.

Creating a good initial guess for the solution is more an art than a science. However, some general guidelines include:

• Have the initial guess satisfy the boundary conditions, since the solution is required to satisfy them as well. If the problem contains unknown parameters, then the initial guess for the parameters also should satisfy the boundary conditions.

• Try to incorporate as much information about the physical problem or expected solution into the initial guess as possible. For example, if the solution is supposed to oscillate or have a certain number of sign changes, then the initial guess should as well.

• Consider the placement of the mesh points (the x-coordinates of the initial guess of the solution). The BVP solvers adapt these points during the solution process, so you do not need to specify too many mesh points. Best practice is to specify a few mesh points placed near where the solution changes rapidly.

• If there is a known, simpler solution on a smaller interval, then use it as an initial guess on a larger interval. Often you can solve a problem as a series of relatively simpler problems, a practice called continuation. With continuation, a series of simple problems are connected by using the solution of one problem as the initial guess to solve the next problem.

Finding Unknown Parameters

Often BVPs involve unknown parameters p that have to be determined as part of solving the problem. The ODE and boundary conditions become

$$\begin{array}{l}y\text{'}=f\left(x,y,p\right)\\ g\left(y\left(a\right),y\left(b\right),p\right)=0\end{array}$$

In this case, the boundary conditions must suffice to determine the values of the parameters p.

You must provide the solver with an initial guess for any unknown parameters. When you call bvpinit to create the structure solinit, specify the initial guess as a vector in the third input argument parameters.

solinit = bvpinit(x,v,parameters)

Additionally, the functions odefun and bcfun that encode the ODE equations and boundary conditions must each have a third argument.

dydx = odefun(x,y,parameters)
res = bcfun(ya,yb,parameters)

While solving the differential equations, the solver adjusts the value of the unknown parameters to satisfy the boundary conditions. The solver returns the final values of these unknown parameters in sol.parameters.

Singular BVPs

bvp4c and bvp5c can solve a class of singular BVPs of the form

$$\begin{array}{c}{y}^{\prime }=\frac{1}{x}Sy+f\left(x,y\right),\\ 0=g\left(y\left(0\right),y\left(b\right)\right).\end{array}$$

The solvers can also accommodate unknown parameters for problems of the form

$$\begin{array}{c}{y}^{\prime }=\frac{1}{x}Sy+f\left(x,y,p\right),\\ 0=g\left(y\left(0\right),y\left(b\right),p\right).\end{array}$$

Singular problems must be posed on an interval [0,b] with b > 0. Use bvpset to pass the constant matrix S to the solver as the value of the 'SingularTerm' option. Boundary conditions at x = 0 must be consistent with the necessary condition for a smooth solution, Sy(0) = 0. The initial guess of the solution also should satisfy this condition.

When you solve a singular BVP, the solver requires that your function odefun(x,y) return only the value of the f(x, y) term in the equation. The term involving S is handled by the solver separately using the 'SingularTerm' option.

BVP Solver Selection

MATLAB includes the solvers bvp4c and bvp5c to solve BVPs. In most cases you can use the solvers interchangeably. The main difference between the solvers is that bvp4c implements a fourth-order formula, while bvp5c implements a fifth-order formula.

The bvp5c function is used exactly like bvp4c, with the exception of the meaning of error tolerances between the two solvers. If S(x) approximates the solution y(x), bvp4c controls the residual |S′(x) – f(x,S(x))|. This approach indirectly controls the true error |y(x) – S(x)|. Use bvp5c to control the true error directly.

Solver	Description
bvp4c	bvp4c is a finite difference code that implements the 3-stage Lobatto IIIa formula. This is a collocation formula, and the collocation polynomial provides a C1-continuous solution that is fourth-order accurate uniformly in the interval of integration. Mesh selection and error control are based on the residual of the continuous solution. The collocation technique uses a mesh of points to divide the interval of integration into subintervals. The solver determines a numerical solution by solving a global system of algebraic equations resulting from the boundary conditions and the collocation conditions imposed on all the subintervals. The solver then estimates the error of the numerical solution on each subinterval. If the solution does not satisfy the tolerance criteria, then the solver adapts the mesh and repeats the process. You must provide the points of the initial mesh, as well as an initial approximation of the solution at the mesh points.
bvp5c	bvp5c is a finite difference code that implements the four-stage Lobatto IIIa formula. This is a collocation formula and the collocation polynomial provides a C1-continuous solution that is fifth-order accurate uniformly in [a,b]. The formula is implemented as an implicit Runge-Kutta formula. bvp5c solves the algebraic equations directly, whereas bvp4c uses analytical condensation. bvp4c handles unknown parameters directly, while bvp5c augments the system with trivial differential equations for unknown parameters.

Evaluating the Solution

The collocation methods implemented in bvp4c and bvp5c produce C¹-continuous solutions over the interval of integration [a,b]. You can evaluate the approximate solution, S(x), at any point in [a,b] using the helper function deval and the structure sol returned by the solver. For example, to evaluate the solution sol at the mesh of points xint, use the command

Sxint = deval(sol,xint)

The deval function is vectorized. For a vector xint, the ith column of Sxint approximates the solution y(xint(i)).

BVP Examples and Files

Several available example files serve as excellent starting points for most common BVP problems. To easily explore and run examples, simply use the Differential Equations Examples app. To run this app, type

odeexamples

To open an individual example file for editing, type

edit exampleFileName.m

To run an example, type

exampleFileName

This table contains a list of the available BVP example files, as well as the solvers and the options they use.

Example File	Solver Used	Options Specified	Description	Example Link
emdenbvp	bvp4c or bvp5c	'SingularTerm'	Emden's equation, a singular BVP	Solve BVP with Singular Term
fsbvp	bvp4c or bvp5c	—	Falkner-Skan BVP on an infinite interval	Verify BVP Consistency Using Continuation
mat4bvp	bvp4c or bvp5c	—	Fourth eigenfunction of Mathieu's equation	Solve BVP with Unknown Parameter
rcbvp	bvp4c and bvp5c	'FJacobian' 'AbsTol' 'RelTol' 'Stats'	Example comparing the errors controlled by bvp4c and bvp5c	Compare bvp4c and bvp5c Solvers (bvp4c) Compare bvp4c and bvp5c Solvers (bvp5c)
shockbvp	bvp4c or bvp5c	'FJacobian' 'BCJacobian' 'Vectorized'	Solution with a shock layer near x = 0	Solve BVP Using Continuation
twobvp	bvp4c	—	BVP with exactly two solutions	Solve BVP with Two Solutions
threebvp	bvp4c or bvp5c	—	Three-point boundary value problem	Solve BVP with Multiple Boundary Conditions

References

[1] Ascher, U., R. Mattheij, and R. Russell. “Numerical Solution of Boundary Value Problems for Ordinary Differential Equations.” Philadelphia, PA: SIAM, 1995, p. 372.

[2] Shampine, L.F., and J. Kierzenka. "A BVP Solver based on residual control and the MATLAB PSE." ACM Trans. Math. Softw. Vol. 27, Number 3, 2001, pp. 299–316.

[3] Shampine, L.F., M.W. Reichelt, and J. Kierzenka. "Solving Boundary Value Problems for Ordinary Differential Equations in MATLAB with bvp4c." MATLAB File Exchange, 2004.

[4] Shampine, L.F., and J. Kierzenka. "A BVP Solver that Controls Residual and Error." J. Numer. Anal. Ind. Appl. Math. Vol. 3(1-2), 2008, pp. 27–41.