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 twopoint 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\text{\hspace{0.17em}}.$$
To specify the boundary conditions for a given BVP, you must:
Write a function of the form
res = bcfun(ya,yb)
, or use the formres = bcfun(ya,yb,p)
if there are unknown parameters involved. You supply this function to the solver as the second input argument. The function returnsres
, 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 isfunction 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 xcoordinates 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 fourthorder formula, while
bvp5c
implements a fifthorder 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 

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 is a finite difference code that
implements the fourstage Lobatto IIIa formula. This is a
collocation formula and the collocation polynomial provides a
C^{1}continuous solution that is
fifthorder accurate uniformly in [a,b]. The formula is implemented
as an implicit RungeKutta 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^{1}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 i
th 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
edit exampleFileName.m
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 



 Emden's equation, a singular BVP  

 —  FalknerSkan BVP on an infinite interval  

 —  Fourth eigenfunction of Mathieu's equation  


 Example comparing the errors controlled by
 Compare bvp4c and bvp5c Solvers ( Compare bvp4c and bvp5c Solvers ( 


 Solution with a shock layer near x = 0  

 —  BVP with exactly two solutions  

 —  Threepoint boundary value problem 
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(12), 2008, pp. 27–41.