Solve boundary value problem — fifthorder method
integrates a system of differential equations of the form y′ =
f(x,y) specified by sol
= bvp5c(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
= bvp5c(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
.
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 bvp5c
with the derivative function, boundary condition function, and initial guess to solve the problem.
sol = bvp5c(@bvpfcn, @bcfcn, solinit);
Plot Solution
plot(sol.x, sol.y, 'o')
Local Functions
Listed here are the local functions that bvp5c
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 %
bvp4c
and bvp5c
SolversSolve 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 %
odefun
— Functions to solveFunctions 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
bvp5c
also can solve problems with singularities in the
solution or multipoint boundary
conditions.
Example: sol = bvp5c(@odefun, @bcfun, solinit)
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 conditionsBoundary 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 = bvp5c(@odefun, @bcfun, solinit)
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 solutionInitial 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 = bvp5c(@odefun, @bcfun, solinit)
Data Types: struct
options
— Option structureOption 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
sol
— Solution structureSolution 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

For multipoint boundary value problems, the boundary conditions are enforced at several points in the interval of integration.
bvp5c
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. bvp5c
enumerates the regions from left to right (from a to
b), with indices starting from 1. In region k,
[a_{k−1},a_{k}],
bvp5c
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.
bvp5c
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.
bvp5c
is a finite difference code that implements the fourstage
Lobatto IIIa formula [1]. 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. Some of the differences between bvp5c
and
bvp4c
are:
bvp5c
solves the algebraic equations directly.
bvp4c
uses analytical condensation.
bvp4c
handles unknown parameters directly.
bvp5c
augments the system with trivial differential equations for
the unknown parameters.
[1] 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.
A modified version of this example exists on your system. Do you want to open this version instead?
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.
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: .
Select web siteYou can also select a web site from the following list:
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.