MATLAB®

BVPs

Function Summary

• BVP Solver

• BVP Helper Functions

• BVP Solver Options

BVP Solver

BVP Helper Functions

BVP Solver Options

An options structure contains named properties whose values are passed to bvp4c, and which affect problem solution. Use these functions to create, alter, or access an options structure.

Back to Top

Boundary Value Problems

The BVP solver is designed to handle systems of ordinary differential equations

where x is the independent variable, y is the dependent variable, and y′ represents the derivative of y with respect to x dy/dx.

See First Order ODEs for general information about ODEs.

Boundary Conditions

In a boundary value problem, the solution of interest satisfies certain boundary conditions. These conditions specify a relationship between the values of the solution at more than one x. In its basic syntax, bvp4c is designed to solve two-point BVPs, i.e., problems where the solution sought on an interval [a, b] must satisfy the boundary conditions

Unlike initial value problems, a boundary value problem may not have a solution, may have a finite number of solutions, or may have infinitely many solutions. As an integral part of the process of solving a BVP, you need to provide a guess for the required solution. The quality of this guess can be critical for the solver performance and even for a successful computation.

There may be other difficulties when solving BVPs, such as problems imposed on infinite intervals or problems that involve singular coefficients. Often BVPs involve unknown parameters p that have to be determined as part of solving the problem

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

Back to Top

BVP Solver

The BVP Solver

The function bvp4c solves two-point boundary value problems for ordinary differential equations (ODEs). It integrates a system of first-order ordinary differential equations

on the interval , subject to general two-point boundary conditions

It can also accommodate other types of BVP problems, such as those that have any of the following:

• Unknown parameters

• Singularities in the solutions

• Multipoint conditions

In this case, the number of boundary conditions must be sufficient to determine the solution and the unknown parameters.

bvp4c produces a solution that is continuous on and has a continuous first derivative there. You can use the function deval and the output of bvp4c to evaluate the solution at specific points on the interval of integration.

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 C¹-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, the solver adapts the mesh and repeats the process. The user must provide the points of the initial mesh as well as an initial approximation of the solution at the mesh points.

BVP Solver Syntax

The basic syntax of the BVP solver is

sol = bvp4c(odefun,bcfun,solinit)

The input arguments are

dydx = odefun(x,y)

res = bcfun(ya,yb)

The output argument sol is a structure created by the solver. In the basic case the structure has fields x, y, yp, and solver.

The structure sol returned by bvp4c contains an additional field if the problem involves unknown parameters:

The function deval uses the output structure sol to evaluate the numerical solution at any point from [a,b].

BVP Solver Options

For more advanced applications, you can specify solver options by passing an input argument options.

sol = bvp4c(odefun,bcfun,solinit,options)

Back to Top

Integrator Options

The default integration properties in the BVP solver bvp4c are selected to handle common problems. In some cases, you can improve solver performance by overriding these defaults. You do this by supplying bvp4c with an options structure that specifies one or more property values.

For example, to change the value of the relative error tolerance of bvp4c from the default value of 1e-3 to 1e-4,

1. Create an options structure using the function bvpset by entering

   options = bvpset('RelTol', 1e-4);

2. Pass the options structure to bvp4c as follows:

   sol = bvp4c(odefun,bcfun,solinit,options)

For a complete description of the available options, see the reference page for bvpset.

Back to Top

Examples

• Mathieu's Equation

• Continuation

• Singular BVPs

• Multipoint BVPs

• Additional Examples

Mathieu's Equation

• Solving the Problem

• Finding Unknown Parameters

• Evaluating the Solution

Solving the Problem.   This example determines the fourth eigenvalue of Mathieu's Equation. It illustrates how to write second-order differential equations as a system of two first-order ODEs and how to use bvp4c to determine an unknown parameter .

The task is to compute the fourth ( ) eigenvalue lambda of Mathieu's equation

Because the unknown parameter is present, this second-order differential equation is subject to three boundary conditions

Note   The demo mat4bvp contains the complete code for this example. The demo uses nested functions to place all functions required by bvp4c in a single M-file. To run this example type mat4bvp at the command line. See BVP Solver Syntax for more information.

1. Rewrite the problem as a first-order system. To use bvp4c, you must rewrite the equations as an equivalent system of first-order differential equations. Using a substitution and , the differential equation is written as a system of two first-order equations

   

   Note that the differential equations depend on the unknown parameter . The boundary conditions become

   

2. Code the system of first-order ODEs. Once you represent the equation as a first-order system, you can code it as a function that bvp4c can use. Because there is an unknown parameter, the function must be of the form

   dydx = odefun(x,y,parameters)
   

   The following code represents the system in the function, mat4ode. Variable q is shared with the outer function:

   function dydx = mat4ode(x,y,lambda)
   dydx = [ y(2)
            -(lambda - 2*q*cos(2*x))*y(1) ];
   end   % End nested function mat4ode
   

3. Code the boundary conditions function. You must also code the boundary conditions in a function. Because there is an unknown parameter, the function must be of the form

   res = bcfun(ya,yb,parameters)
   

   The code below represents the boundary conditions in the function, mat4bc.

   function res = mat4bc(ya,yb,lambda)
   res = [  ya(2) 
            yb(2) 
           ya(1)-1 ];
   

4. Create an initial guess. To form the guess structure solinit with bvpinit, you need to provide initial guesses for both the solution and the unknown parameter.

   The function mat4init provides an initial guess for the solution. mat4init uses because this function satisfies the boundary conditions and has the correct qualitative behavior (the correct number of sign changes).

   function yinit = mat4init(x)
   yinit = [  cos(4*x)
             -4*sin(4*x) ];
   

   In the call to bvpinit, the third argument, lambda, provides an initial guess for the unknown parameter .

   lambda = 15;
   solinit = bvpinit(linspace(0,pi,10),@mat4init,lambda);
   

   This example uses @ to pass mat4init as a function handle to bvpinit.

   Note   See the function_handle (@), func2str, and str2func reference pages, and see Function Handles in MATLAB Programming Fundamentals.

5. Apply the BVP solver. The mat4bvp example calls bvp4c with the functions mat4ode and mat4bc and the structure solinit created with bvpinit.

   sol = bvp4c(@mat4ode,@mat4bc,solinit);

6. View the results. Complete the example by displaying the results:

   1. Print the value of the unknown parameter found by bvp4c.

          fprintf('Fourth eigenvalue is approximately %7.3f.\n',...
                  sol.parameters)

   2. Use deval to evaluate the numerical solution at 100 equally spaced points in the interval , and plot its first component. This component approximates .

          xint = linspace(0,pi);
          Sxint = deval(sol,xint);
          plot(xint,Sxint(1,:))
          axis([0 pi -1 1.1])
          title('Eigenfunction of Mathieu''s equation.') 
          xlabel('x')
          ylabel('solution y')

      The following plot shows the eigenfunction associated with the final eigenvalue = 17.097.

      

Finding Unknown Parameters.   The bvp4c solver can find unknown parameters for problems of the form

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

solinit = bvpinit(x,v,parameters)

The bvp4c function arguments odefun and bcfun must each have a third argument.

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

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

Evaluating the Solution.   The collocation method implemented in bvp4c produces a C¹-continuous solution over the whole interval of integration . You can evaluate the approximate solution, , at any point in using the helper function deval and the structure sol returned by bvp4c.

Sxint = deval(sol,xint)

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

Continuation

• Introduction

• Using Continuation to Solve a BVP

• Using Continuation to Verify Consistency

Introduction.   To solve a boundary value problem, you need to provide an initial guess for the solution. The quality of your initial guess can be critical to the solver performance, and to being able to solve the problem at all. However, coming up with a sufficiently good guess can be the most challenging part of solving a boundary value problem. Certainly, you should apply the knowledge of the problem's physical origin. Often a problem can be solved as a sequence of relatively simpler problems, i.e., a continuation.

This example shows how to use continuation to:

• Solve a difficult BVP

• Verify a solution's consistent behavior

Using Continuation to Solve a BVP.   This example solves the differential equation

for , on the interval [-1 1], with boundary conditions and . For , the solution has a transition layer at . Because of this rapid change in the solution for small values of , the problem becomes difficult to solve numerically.

The example solves the problem as a sequence of relatively simpler problems, i.e., a continuation. The solution of one problem is used as the initial guess for solving the next problem.

Note   The demo shockbvp contains the complete code for this example. The demo uses nested functions to place all required functions in a single M-file. To run this example type shockbvp at the command line.

Note   This problem appears in [1] to illustrate the mesh selection capability of a well established BVP code COLSYS.

1. Code the ODE and boundary condition functions. Code the differential equation and the boundary conditions as functions that bvp4c can use:

   The code below represents the differential equation and the boundary conditions in the functions shockODE and shockBC. Note that shockODE is vectorized to improve solver performance. The additional parameter is represented by e and is shared with the outer function.

   function dydx = shockODE(x,y)
   pix = pi*x;
   dydx = [  y(2,:)
            -x/e.*y(2,:) - pi^2*cos(pix) - pix/e.*sin(pix) ];
   end %   End nested function shockODE
   
   function res = shockBC(ya,yb)
   res = [ ya(1)+2 
           yb(1)   ];
   end %   End nested function shockBC
   

2. Provide analytical partial derivatives. For this problem, the solver benefits from using analytical partial derivatives. The code below represents the derivatives in functions shockJac and shockBCJac.

   function jac = shockJac(x,y)
   jac = [ 0   1
           0 -x/e ];
   end %   End nested function shockJac
   
   function [dBCdya,dBCdyb] = shockBCJac(ya,yb)
   dBCdya = [ 1 0
              0 0 ];
   dBCdyb = [ 0 0
              1 0 ];
   end %   End nested function shockBCJac
   

   shockJac shares e with the outer function.

   Tell bvp4c to use these functions to evaluate the partial derivatives by setting the options FJacobian and BCJacobian. Also set 'Vectorized' to 'on' to indicate that the differential equation function shockODE is vectorized.

   options = bvpset('FJacobian',@shockJac,...
                    'BCJacobian',@shockBCJac,...
                    'Vectorized','on');
   

3. Create an initial guess. You must provide bvp4c with a guess structure that contains an initial mesh and a guess for values of the solution at the mesh points. A constant guess of and , and a mesh of five equally spaced points on [-1 1] suffice to solve the problem for . Use bvpinit to form the guess structure.

   sol = bvpinit([-1 -0.5 0 0.5 1],[1 0]);
   

4. Use continuation to solve the problem. To obtain the solution for the parameter , the example uses continuation by solving a sequence of problems for . The solver bvp4c does not perform continuation automatically, but the code's user interface has been designed to make continuation easy. The code uses the output sol that bvp4c produces for one value of e as the guess in the next iteration.

   e = 0.1; 
   for i=2:4 
       e = e/10; 
       sol = bvp4c(@shockODE,@shockBC,sol,options); 
   end

5. View the results. Complete the example by displaying the final solution

   plot(sol.x,sol.y(1,:))
   axis([-1 1 -2.2 2.2])
   title(['There is a shock at x = 0 when \epsilon = '... 
         sprintf('%.e',e) '.']) 
   xlabel('x')
   ylabel('solution y')

   

Using Continuation to Verify Consistency.   Falkner-Skan BVPs arise from similarity solutions of viscous, incompressible, laminar flow over a flat plate. An example is

for on the interval with boundary conditions , , and .

The BVP cannot be solved on an infinite interval, and it would be impractical to solve it for even a very large finite interval. So, the example tries to solve a sequence of problems posed on increasingly larger intervals to verify the solution's consistent behavior as the boundary approaches .

The example imposes the infinite boundary condition at a finite point called infinity. The example then uses continuation in this end point to get convergence for increasingly larger values of infinity. It uses bvpinit to extrapolate the solution sol for one value of infinity as an initial guess for the new value of infinity. The plot of each successive solution is superimposed over those of previous solutions so they can easily be compared for consistency.

Note   The demo fsbvp contains the complete code for this example. The demo uses nested functions to place all required functions in a single M-file. To run this example type fsbvp at the command line.

1. Code the ODE and boundary condition functions. Code the differential equation and the boundary conditions as functions that bvp4c can use. The problem parameter beta is shared with the outer function.

   function dfdeta = fsode(eta,f)
   dfdeta = [ f(2)
              f(3)
             -f(1)*f(3) - beta*(1 - f(2)^2) ];
   end %   End nested function fsode
   
   function res = fsbc(f0,finf)
   res = [f0(1)
          f0(2)
          finf(2) - 1];
   end %   End nested function fsbc
   

2. Create an initial guess. You must provide bvp4c with a guess structure that contains an initial mesh and a guess for values of the solution at the mesh points. A crude mesh of five points and a constant guess that satisfies the boundary conditions are good enough to get convergence when infinity = 3.

   infinity = 3;
   maxinfinity = 6;
   
   solinit = bvpinit(linspace(0,infinity,5),[0 0 1]);
   

3. Solve on the initial interval. The example obtains the solution for infinity = 3. It then prints the computed value of for comparison with the value reported by Cebeci and Keller [2]:

   sol = bvp4c(@fsode,@fsbc,solinit);
   eta = sol.x;
   f = sol.y;
   
   fprintf('\n');
   fprintf('Cebeci & Keller report that f''''(0) = 0.92768.\n')
   fprintf('Value computed using infinity = %g is %7.5f.\n', ...
            infinity,f(3,1))
   

   The example prints

   Cebeci & Keller report that f''(0) = 0.92768.
   Value computed using infinity = 3 is 0.92915.
   

4. Setup the figure and plot the initial solution.

   figure
   plot(eta,f(2,:),eta(end),f(2,end),'o');
   axis([0 maxinfinity 0 1.4]);
   title('Falkner-Skan equation, positive wall shear, ...
         \beta = 0.5.')
   xlabel('\eta')
   ylabel('df/d\eta')
   hold on
   drawnow 
   shg 
   

   

5. Use continuation to solve the problem and plot subsequent solutions. The example then solves the problem for infinity = 4, 5, 6. It uses bvpinit to extrapolate the solution sol for one value of infinity as an initial guess for the next value of infinity. For each iteration, the example prints the computed value of and superimposes a plot of the solution in the existing figure.

   for Bnew = infinity+1:maxinfinity
     
     solinit = bvpinit(sol,[0 Bnew]); % Extend solution to Bnew.
     sol = bvp4c(@fsode,@fsbc,solinit);
     eta = sol.x;
     f = sol.y;
     
     fprintf('Value computed using infinity = %g is %7.5f.\n', ...
              Bnew,f(3,1))
     plot(eta,f(2,:),eta(end),f(2,end),'o');
     drawnow
     
   end
   hold off
   

   The example prints

   Value computed using infinity = 4 is 0.92774.
   Value computed using infinity = 5 is 0.92770.
   Value computed using infinity = 6 is 0.92770.
   

   Note that the values approach 0.92768 as reported by Cebeci and Keller. The superimposed plots confirm the consistency of the solution's behavior.

   

Singular BVPs

• Introduction

• Emden's equation

Introduction.   The function bvp4c solves a class of singular BVPs of the form

(6-1)

It can also accommodate unknown parameters for problems of the form

Singular problems must be posed on an interval with . Use bvpset to pass the constant matrix to bvp4c as the value of the 'SingularTerm' integration property. Boundary conditions at must be consistent with the necessary condition for a smooth solution, . An initial guess should also satisfy this necessary condition.

When you solve a singular BVP using

sol = bvp4c(@odefun,@bcfun,solinit,options)

bvp4c requires that your function odefun(x,y) return only the value of the term in Equation 5-2.

Emden's equation.   Emden's equation arises in modeling a spherical body of gas. The PDE of the model is reduced by symmetry to the ODE

on an interval . The coefficient is singular at , but symmetry implies the boundary condition . With this boundary condition, the term

is well-defined as approaches 0. For the boundary condition , this BVP has the analytical solution

Note   The demo emdenbvp contains the complete code for this example. The demo uses subfunctions to place all required functions in a single M-file. To run this example type emdenbvp at the command line.

1. Rewrite the problem as a first-order system and identify the singular term. Using a substitution and , write the differential equation as a system of two first-order equations

   

   The boundary conditions become

   

   Writing the ODE system in a vector-matrix form

   

   the terms of Equation 5-2 are identified as

   

   and

   

2. Code the ODE and boundary condition functions. Code the differential equation and the boundary conditions as functions that bvp4c can use.

   function dydx = emdenode(x,y)
   dydx = [  y(2) 
            -y(1)^5 ];
   function res = emdenbc(ya,yb)
   res = [ ya(2)
           yb(1) - sqrt(3)/2 ];
   

3. Setup integration properties. Use the matrix as the value of the 'SingularTerm' integration property.

   S = [0,0;0,-2];
   options = bvpset('SingularTerm',S);
   

4. Create an initial guess. This example starts with a mesh of five points and a constant guess for the solution.

   

   Use bvpinit to form the guess structure

   guess = [sqrt(3)/2;0];
   solinit = bvpinit(linspace(0,1,5),guess);
   

5. Solve the problem. Use the standard bvp4c syntax to solve the problem.

   sol = bvp4c(@emdenode,@emdenbc,solinit,options);
   

6. View the results. This problem has an analytical solution

   

   The example evaluates the analytical solution at 100 equally spaced points and plots it along with the numerical solution computed using bvp4c.

   x = linspace(0,1);
   truy = 1 ./ sqrt(1 + (x.^2)/3);
   plot(x,truy,sol.x,sol.y(1,:),'ro');
   title('Emden problem -- BVP with singular term.')
   legend('Analytical','Computed');
   xlabel('x');
   ylabel('solution y');
   

   

Multipoint BVPs

In multipoint boundary value problems, the solution of interest satisfies conditions at points inside the interval of integration. The bvp4c function is useful in solving such problems.

The following example shows how the multipoint capability in bvp4c can improve efficiency when you are solving a nonsmooth problem. The following equations are solved on for constant parameters n, , , and . These are subject to boundary conditions v(0) = 0 and :

v' = (C - 1)/n
C' = (v * C - min(x,1))/

The term min(x,1) is not smooth at xc = 1, and this can affect the solver's efficiency. By introducing an interface point at xc = 1, smooth solutions can be obtained on [0,1] and [1, ]. To get a continuous solution over the entire interval [0, ], the example imposes matching conditions at the interface.

Note   The demo threebvp contains the complete code for this example and solves the problem for = 2, n = 0.05, and several values of . The demo uses nested functions to place all functions required by bvp4c in a single M-file and to communicate problem parameters efficiently. To run this example, type threebvp at the MATLAB command prompt.

The demo takes you through the following steps:

1. Determine the interfaces and divide the interval of integration into regions. Introducing an interface point at xc = 1 divides the problem into two regions in which the solutions remain smooth. The differential equations for the two regions are

   Region 1:

   v' = (C - 1)/n 
   C' = (v * C - x)/
   
   

   Region 2:

   v' = (C - 1)/n 
   C' = (v * C - 1)/
   
   

   Note that the interface xc = 1 is included in both regions. At xc = 1, bvp4c produces a left and right solution. These solutions are denoted as v(1-), C(1-) and v(1+), C(1+) respectively.

2. Determine the boundary conditions. Solving two first-order differential equations in two regions requires imposing four boundary conditions. Two of these conditions come from the original formulation; the others enforce the continuity of the solution across the interface xc = 1:

   v(0) = 0
   C(
   ) - 1 = 0
   v(1-) - v(1+) = 0
   C(1-) - C(1+) = 0
   

   Here, v(1-), C(1-) and v(1+), C(1+) denote the left and right solution at the interface.

3. Code the derivative function. In the derivative function, y(1) corresponds to v(x), and y(2) corresponds to C(x). The additional input argument region identifies the region in which the derivative is evaluated. bvp4c enumerates regions from left to right, starting with 1. Note that the problem parameters n and are shared with the outer function:

   function dydx = f(x,y,region)
      dydx = zeros(2,1);
      dydx(1) = (y(2) - 1)/n;
   
      % The definition of C'(x) depends on the region.
      switch region
         case 1                                % x in [0 1]
            dydx(2) = (y(1)*y(2) - x)/
   ; 
         case 2                                % x in [1 
   ]
            dydx(2) = (y(1)*y(2) - 1)/
   ; 
       end
   end                % End nested function f
   

4. Code the boundary conditions function. For multipoint BVPs, the arguments of the boundary conditions function, YL and YR, become matrices. In particular, the kth column YL(:,k) represents the solution at the left boundary of the kth region. Similarly, YR(:,k) represents the solution at the right boundary of the kth region.

   In the example, y(0) is approximated by YL(:,1), while y( ) is approximated by YR(:,end). Continuity of the solution at the internal interface requires that YR(:,1) = YL(:,2). Nested function bc computes the residual in the boundary conditions:

   function res = bc(YL,YR)
      res = [YL(1,1)               % v(0) = 0
             YR(1,1) - YL(1,2)     % Continuity of v(x) at x=1
             YR(2,1) - YL(2,2)     % Continuity of C(x) at x=1
             YR(2,end) - 1];       % C(
   ) = 1
   end                             % End nested function bc
   

5. Create an initial guess. For multipoint BVPs, when creating an initial guess using bvpinit, use double entries in xinit for the interface point xc. This example uses a constant guess yinit = [1;1]:

   xc = 1;
   xinit = [0, 0.25, 0.5, 0.75, xc, xc, 1.25, 1.5, 1.75, 2];
   solinit = bvpinit(xinit,yinit)
   

   For multipoint BVPs, you can use different guesses in different regions. To do that, you specify the initial guess for y as a function using the following syntax:

   solinit = bvpinit(xinit,@yinitfcn)
   

   The initial guess function must have the following general form:

   function y = yinitfcn(x,region)
      switch region
      case 1            % x in [0, 1]
         y = [1;1];     % initial guess for y(x) 
   
      case 2 % x in [1, 
   ]
         y = [1;1];     % initial guess for y(x), 
   
   end
   

6. Apply the solver. The bvp4c function uses the same syntax for multipoint BVPs as it does for two-point BVPs:

   sol = bvp4c(@f,@bc,solinit);
   

   The mesh points returned in sol.x are adapted to the solution behavior, but the mesh still includes a double entry for the interface point xc = 1. Corresponding columns of sol.y represent the left and right solution at xc.

7. View the results. Using deval, the solution can be evaluated at any point in the interval of integration.

   Note that, with the left and right values computed at the interface, the solution is not uniquely defined at xc = 1. When evaluating the solution exactly at the interface, deval issues a warning and returns the average of the left and right solution values. Call deval at xc-eps(xc) and xc+eps(xc) to get the limit values at xc.

   The example plots the solution approximated at the mesh points selected by the solver:

   plot(sol.x,sol.y(1,:),sol.x,sol.y(2,:),'--')
   legend('v(x)','C(x)')
   title('A three-point BVP solved with BVP4C')
   xlabel(['\
    = ',num2str(
   ), ...
           ', \
    = ',num2str(
   ),'.'])
   ylabel('v and C') 
   

   

Additional Examples

The following additional examples are available. Type

edit examplename

to view the code and

examplename

to run the example.

Additional examples are provided by "Tutorial on Solving BVPs with BVP4C," available at http://www.mathworks.com/bvp_tutorial.

Back to Top

DDEs

PDEs