Main Content

## 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\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 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.

SolverDescription

`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 C1-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`
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 UsedOptions Specified

Description

Example Link

`emdenbvp`

`bvp4c` or `bvp5c`

• `'SingularTerm'`

Emden's equation, a singular BVP

`fsbvp`

`bvp4c` or `bvp5c`

Falkner-Skan BVP on an infinite interval

`mat4bvp`

`bvp4c` or `bvp5c`

Fourth eigenfunction of Mathieu's equation

`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

`twobvp`

`bvp4c`

BVP with exactly two solutions

`threebvp`

`bvp4c` or `bvp5c`

Three-point 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(1-2), 2008, pp. 27–41.

Download ebook