bvp4c
Solve boundary value problem — fourthorder method
Description
integrates a system of differential equations of the form y′ =
f(x,y) specified by sol
= bvp4c(odefun
,bcfun
,solinit
)odefun
, subject to the boundary conditions
described by bcfun
and the initial solution guess
solinit
. Use the bvpinit
function to create the initial guess solinit
, which
also defines the points at which the boundary conditions in bcfun
are
enforced.
also uses the integration settings defined by sol
= bvp4c(odefun
,bcfun
,solinit
,options
)options
, which is an
argument created using the bvpset
function. For example, use the
AbsTol
and RelTol
options to specify absolute and
relative error tolerances, or the FJacobian
option to provide the
analytical partial derivatives of odefun
.
Examples
Solve SecondOrder BVP
Solve a secondorder BVP in MATLAB® using functions. For this example, use the secondorder equation
${{\mathit{y}}^{\prime}}^{\prime}+\mathit{y}=0$.
The equation is defined on the interval $\left[0,\pi /2\right]$ subject to the boundary conditions
$\mathit{y}\left(0\right)=0$,
$\mathit{y}\left(\pi /2\right)=2$.
To solve this equation in MATLAB, you need to write a function that represents the equation as a system of firstorder equations, a function for the boundary conditions, and a function for the initial guess. Then the BVP solver uses these three inputs to solve the equation.
Code Equation
Write a function that codes the equation. Use the substitutions ${\mathit{y}}_{1}=\mathit{y}$ and ${\mathit{y}}_{2}={\mathit{y}}^{\prime}$ to rewrite the equation as a system of firstorder equations.
${{\mathit{y}}_{1}}^{\prime}={\mathit{y}}_{2}$,
${{\mathit{y}}_{2}}^{\prime}={\mathit{y}}_{1}$.
The corresponding function is
function dydx = bvpfcn(x,y) dydx = zeros(2,1); dydx = [y(2) y(1)]; end
Note: All functions are included at the end of the example as local functions.
Code Boundary Conditions
Write a function that codes the boundary conditions in the form $\mathit{g}\left(\mathit{y}\left(\mathit{a}\right),\mathit{y}\left(\mathit{b}\right)\right)=0$. In this form the boundary conditions are
$\mathit{y}\left(0\right)=0$,
$\mathit{y}\left(\pi /2\right)2=0$.
The corresponding function is
function res = bcfcn(ya,yb) res = [ya(1) yb(1)2]; end
Create Initial Guess
Use the bvpinit
function to create an initial guess for the solution of the equation. Since the equation relates ${{\mathit{y}}^{\prime}}^{\prime}$ to $\mathit{y}$, a reasonable guess is that the solution involves trigonometric functions. Use a mesh of five points in the interval of integration. The first and last values in the mesh are where the solver applies the boundary conditions.
The function for the initial guess accepts $\mathit{x}$ as an input and returns a guess for the value of ${\mathit{y}}_{1}$ and ${\mathit{y}}_{2}$. The function is
function g = guess(x) g = [sin(x) cos(x)]; end
xmesh = linspace(0,pi/2,5); solinit = bvpinit(xmesh, @guess);
Solve Equation
Use bvp4c
with the derivative function, boundary condition function, and initial guess to solve the problem.
sol = bvp4c(@bvpfcn, @bcfcn, solinit);
Plot Solution
plot(sol.x, sol.y, 'o')
Local Functions
Listed here are the local functions that bvp4c
uses to solve the equation.
function dydx = bvpfcn(x,y) % equation to solve dydx = zeros(2,1); dydx = [y(2) y(1)]; end % function res = bcfcn(ya,yb) % boundary conditions res = [ya(1) yb(1)2]; end % function g = guess(x) % initial guess for y and y' g = [sin(x) cos(x)]; end %
Compare bvp4c
and bvp5c
Solvers
Solve a BVP at a crude error tolerance with two different solvers and compare the results.
Consider the secondorder ODE
${{\mathit{y}}^{\prime}}^{\prime}+\frac{2}{\mathit{x}}{\mathit{y}}^{\prime}+\frac{1}{{\mathit{x}}^{4}}\mathit{y}=0$.
The equation is defined on the interval $\left[\frac{1}{3\pi},1\right]$ subject to the boundary conditions
$\mathit{y}\left(\frac{1}{3\pi}\right)=0$,
$\mathit{y}\left(1\right)=\mathrm{sin}\left(1\right)$.
To solve this equation in MATLAB®, you need to write a function that represents the equation as a system of firstorder equations, write a function for the boundary conditions, set some option values, and create an initial guess. Then the BVP solver uses these four inputs to solve the equation.
Code Equation
With the substitutions ${\mathit{y}}_{1}=\mathit{y}$ and ${\mathit{y}}_{2}={\mathit{y}}^{\prime}$, you can rewrite the ODE as a system of firstorder equations
${{\mathit{y}}_{1}}^{\prime}={\mathit{y}}_{2}$,
${{\mathit{y}}_{2}}^{\prime}=\frac{2}{\mathit{x}}{\mathit{y}}_{2}\frac{1}{{\mathit{x}}^{4}}{\mathit{y}}_{1}$.
The corresponding function is
function dydx = bvpfcn(x,y) dydx = [y(2) 2*y(2)/x  y(1)/x^4]; end
Note: All functions are included at the end of the example as local functions.
Code Boundary Conditions
The boundary condition function requires that the boundary conditions are in the form $\mathit{g}\left(\mathit{y}\left(\mathit{a}\right),\mathit{y}\left(\mathit{b}\right)\right)=0$. In this form, the boundary conditions are
$\mathit{y}\left(\frac{1}{3\pi}\right)=0$,
$\mathit{y}\left(1\right)\mathrm{sin}\left(1\right)=0$.
The corresponding function is
function res = bcfcn(ya,yb) res = [ya(1) yb(1)sin(1)]; end
Set Options
Use bvpset
to turn on the display of solver statistics, and specify crude error tolerances to highlight the difference in error control between the solvers. Also, for efficiency, specify the analytical Jacobian
$\mathit{J}=\frac{\partial {\mathit{f}}_{\mathit{i}}}{\partial \mathit{y}}=\left[\begin{array}{cc}\frac{\partial {\mathit{f}}_{1}}{\partial {\mathit{y}}_{1}}& \frac{\partial {\mathit{f}}_{1}}{\partial {\mathit{y}}_{2}}\\ \frac{\partial {\mathit{f}}_{2}}{\partial {\mathit{y}}_{1}}& \frac{\partial {\mathit{f}}_{2}}{\partial {\mathit{y}}_{2}}\end{array}\right]=\left[\begin{array}{cc}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}0& \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}1\\ \frac{1}{{\mathit{x}}^{4}}& \frac{2}{\mathit{x}}\end{array}\right]$.
The corresponding function that returns the value of the Jacobian is
function dfdy = jac(x,y) dfdy = [0 1 1/x^4 2/x]; end
opts = bvpset('FJacobian',@jac,'RelTol',0.1,'AbsTol',0.1,'Stats','on');
Create Initial Guess
Use bvpinit
to create an initial guess of the solution. Specify a constant function as the initial guess with an initial mesh of 10 points in the interval $\left[1/3\pi ,1\right]$.
xmesh = linspace(1/(3*pi), 1, 10); solinit = bvpinit(xmesh, [1; 1]);
Solve Equation
Solve the equation with both bvp4c
and bvp5c
.
sol4c = bvp4c(@bvpfcn, @bcfcn, solinit, opts);
The solution was obtained on a mesh of 9 points. The maximum residual is 9.794e02. There were 157 calls to the ODE function. There were 28 calls to the BC function.
sol5c = bvp5c(@bvpfcn, @bcfcn, solinit, opts);
The solution was obtained on a mesh of 11 points. The maximum error is 6.742e02. There were 244 calls to the ODE function. There were 29 calls to the BC function.
Plot Results
Plot the results of the two calculations for ${\mathit{y}}_{1}$ with the analytic solution for comparison. The analytic solution is
${\mathit{y}}_{1}=\mathrm{sin}\left(\frac{1}{\mathit{x}}\right)$,
${\mathit{y}}_{2}=\frac{1}{{\mathit{x}}^{2}}\mathrm{cos}\left(\frac{1}{\mathit{x}}\right)$.
xplot = linspace(1/(3*pi),1,200); yplot = [sin(1./xplot); cos(1./xplot)./xplot.^2]; plot(xplot,yplot(1,:),'k',sol4c.x,sol4c.y(1,:),'r*',sol5c.x,sol5c.y(1,:),'bo') title('Comparison of BVP Solvers with Crude Error Tolerance') legend('True','BVP4C','BVP5C') xlabel('x') ylabel('solution y')
The plot confirms that bvp5c
directly controls the true error in the calculation, while bvp4c
controls it only indirectly. At more stringent error tolerances, this difference between the solvers is not as apparent.
Local Functions
Listed here are the local functions that the BVP solvers use to solve the problem.
function dydx = bvpfcn(x,y) % equation to solve dydx = [y(2) 2*y(2)/x  y(1)/x^4]; end % function res = bcfcn(ya,yb) % boundary conditions res = [ya(1) yb(1)sin(1)]; end % function dfdy = jac(x,y) % analytical jacobian for f dfdy = [0 1 1/x^4 2/x]; end %
Input Arguments
odefun
— Functions to solve
function handle
Functions to solve, specified as a function handle that defines the functions to be
integrated. odefun
and bcfun
must accept the same
number of input arguments.
To code odefun
, use the functional signature dydx =
odefun(x,y)
for a scalar x
and column vector
y
. The return value dydt
is a column vector of
data type single
or double
that corresponds to f(x,y). odefun
must accept both input arguments
x
and y
, even if one of the arguments is not
used in the function.
For example, to solve $$y\text{'}=5y3$$, use the function:
function dydt = odefun(t,y) dydt = 5*y3; end
For a system of equations, the output of odefun
is a vector. Each
element in the vector is the solution to one equation. For example, to solve
$$\begin{array}{l}y{\text{'}}_{1}={y}_{1}+2{y}_{2}\\ y{\text{'}}_{2}=3{y}_{1}+2{y}_{2}\end{array}$$
use the function:
function dydt = odefun(t,y) dydt = zeros(2,1); dydt(1) = y(1)+2*y(2); dydt(2) = 3*y(1)+2*y(2); end
bvp4c
also can solve problems with singularities in the
solution or multipoint boundary
conditions.
Example: sol = bvp4c(@odefun, @bcfun, solinit)
Unknown Parameters
If the BVP being solved includes unknown parameters, you instead can use the
functional signature dydx = odefun(x,y,p)
, where
p
is a vector of parameter values. When you use this functional
signature the BVP solver calculates the values of the unknown parameters starting from
the initial guess for the parameter values provided in
solinit
.
Data Types: function_handle
bcfun
— Boundary conditions
function handle
Boundary conditions, specified as a function handle that computes the residual error
in the boundary conditions. odefun
and bcfun
must
accept the same number of input arguments.
To code bcfun
, use the functional signature res =
bcfun(ya,yb)
for column vectors ya
and
yb
. The return value res
is a column vector of
data type single
or double
that corresponds to the
residual value of the boundary conditions at the boundary points.
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
ya(1)1
should
be 0
at the point x = a. Similarly, since
y(b) = 0, the residual value of yb(1)
should be
0
at the point x = b.
The boundary points x = a and x = b where the
boundary conditions are enforced are defined in the initial guess structure
solinit
. For twopoint boundary value problems, a =
solinit.x(1)
and b = solinit.x(end)
.
Example: sol = bvp4c(@odefun, @bcfun, solinit)
Unknown Parameters
If the BVP being solved includes unknown parameters, you instead can use the
functional signature res = bcfun(ya,yb,p)
, where
p
is a vector of parameter values. When you use this functional
signature the BVP solver calculates the values of the unknown parameters starting from
the initial guess for the parameter values provided in
solinit
.
Data Types: function_handle
solinit
— Initial guess of solution
structure
Initial guess of solution, specified as a structure. Use bvpinit
to create solinit
.
Unlike initial value problems, a boundary value problem can have no solution, a finite number of solutions, or 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. For some guidelines on creating a good initial guess, see Initial Guess of Solution.
Example: sol = bvp4c(@odefun, @bcfun, solinit)
Data Types: struct
options
— Option structure
structure
Option structure. Use the bvpset
function to create or modify the
options structure.
Example: options = bvpset('RelTol',1e5,'Stats','on')
specifies a
relative error tolerance of 1e5
and turns on the display of solver
statistics.
Data Types: struct
Output Arguments
sol
— Solution structure
structure
Solution structure. You can access the fields in sol
with
dotindexing, such as sol.field1
. The solution
(sol.x
,sol.y
) is continuous on the interval of
integration defined in the initial mesh solinit.x
and has a
continuous first derivative there. You can use sol
with the deval
function to evaluate the solution at other points in the
interval.
The structure sol
has these fields.
Field  Description 

 Mesh selected by 
 Approximation to y(x) at the mesh points of

 Approximation to y′(x) at the mesh points of

 Final values for the unknown parameters specified in



 Computational cost statistics related to the solution: the number
of mesh points, residual error, and number of calls to

More About
Multipoint Boundary Value Problems
For multipoint boundary value problems, the boundary conditions are enforced at several points in the interval of integration.
bvp4c
can solve multipoint boundary value problems where a = a_{0} < a_{1} < a_{2} < ...< a_{n} = b are boundary points in the interval
[a,b]. The points a_{1},a_{2},...,a_{n−1} represent interfaces that divide
[a,b] into regions. bvp4c
enumerates the regions from left to right (from a to
b), with indices starting from 1. In region k,
[a_{k−1},a_{k}],
bvp4c
evaluates the derivative as
yp = odefun(x,y,k)
In the boundary conditions function bcfun(yleft,yright)
,
yleft(:,k)
is the solution at the left boundary of
[a_{k−1},a_{k}].
Similarly, yright(:,k)
is the solution at the right boundary of region
k. In particular, yleft(:,1) = y(a)
and
yright(:,end) = y(b)
.
When you create an initial guess with bvpinit
, use double entries
in xinit
for each interface point. See the reference page for bvpinit
for more information.
If yinit
is a function, bvpinit
calls y
= yinit(x,k)
to get an initial guess for the solution at x
in
region k
. In the solution structure sol
returned by
bpv4c
, sol.x
has double entries for each interface
point. The corresponding columns of sol.y
contain the left and right
solution at the interface, respectively.
See Solve BVP with Multiple Boundary Conditions for an example that solves a threepoint boundary value problem.
Singular Boundary Value Problems
bvp4c
solves a class of singular boundary value
problems, including problems with unknown parameters p
, of the
form
$$\begin{array}{l}y\text{'}=S\frac{y}{x}+f\left(x,y,p\right),\\ 0=bc\left(y\left(0\right),y\left(b\right),p\right).\end{array}$$
The interval is required to be [0, b] with b >
0. Often such problems arise when computing a smooth solution of ODEs that result from
partial differential equations (PDEs) due to cylindrical or spherical symmetry. For singular
problems, you specify the (constant) matrix S
as the value of the
'SingularTerm'
option of bvpset
, and odefun
evaluates only
f(x,y,p). The
boundary conditions and initial guess must be consistent with the necessary condition for
smoothness S·y(0) = 0.
See Solve BVP with Singular Term for an example that solves a singular boundary value problem.
Algorithms
bvp4c
is a finite difference code that implements the threestage
Lobatto IIIa formula [1], [2]. This is a collocation
formula and the collocation polynomial provides a C^{1}continuous
solution that is fourthorder 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, 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.
References
[1] 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.
[2] 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.
Extended Capabilities
ThreadBased Environment
Run code in the background using MATLAB® backgroundPool
or accelerate code with Parallel Computing Toolbox™ ThreadPool
.
This function fully supports threadbased environments. For more information, see Run MATLAB Functions in ThreadBased Environment.
Version History
Open Example
You have a modified version of this example. Do you want to open this example with your edits?
MATLAB Command
You clicked a link that corresponds to this MATLAB command:
Run the command by entering it in the MATLAB Command Window. Web browsers do not support MATLAB commands.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list:
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
 América Latina (Español)
 Canada (English)
 United States (English)
Europe
 Belgium (English)
 Denmark (English)
 Deutschland (Deutsch)
 España (Español)
 Finland (English)
 France (Français)
 Ireland (English)
 Italia (Italiano)
 Luxembourg (English)
 Netherlands (English)
 Norway (English)
 Österreich (Deutsch)
 Portugal (English)
 Sweden (English)
 Switzerland
 United Kingdom (English)