]> 5.10 Using Matlab for solving ODEs: boundary value problems
Skip navigation.
This site will look best in a browser that supports web standards, but it is accessible to any browser or Internet device. We recommend the Mozilla Firefox browser.
University of OxfordSolving ODEs in MATLAB

Unit navigation

5.10 Using Matlab for solving ODEs: boundary value problems

Problem definition

Suppose we wish to solve the system of equations d y d x = f ( x , y ), with conditions applied at two different points x = a and x = b.

More commonly, problems of this sort will be written as a higher-order (that is, a second-order) ODE with derivative boundary conditions. We can reduce any such equation to a system of first-order equations, however, so we only need to consider how to solve systems of first-order equations.

Such problems are known as ‘Boundary Value Problems’ (BVPs).

For a system to be well defined, there should be as many conditions as there are first-order equations. For example, to solve two second-order ODEs you would need four conditions, as this system would equate to one with four first-order ODEs.

Matlab commands

Suppose we wish to solve the system of n equations, d y d x = f ( x , y ) , with conditions applied at two different points x = a and x = b .

In order to use the inbuilt MATLAB ODE solvers, you need to follow the steps below:

  1. Construct a function (here called deriv) which has input arguments x MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLnhiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=xfr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamiEaaaa@36EA@ and y = ( y 1 , , y n ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLnhiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=xfr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaaCyEaiabg2da9iaacIcacaWG5bWaaSbaaSqaaiaaigdaaeqaaOGaaiilaiabl+UimjaacYcacaWG5bWaaSbaaSqaaiaad6gaaeqaaOGaaiykaaaa@40B2@ and returns the value of the derivative d y d x MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLnhiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=xfr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaWaaSaaaeaacaWGKbGaaCyEaaqaaiaadsgacaWG4baaaaaa@39CE@ , that is f ( x , y ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLnhiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=xfr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaaCOzaiaacIcacaWG4bGaaiilaiaahMhacaGGPaaaaa@3AE4@ .
  2. Construct a function (here called bcs) which has input arguments y ( a ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLnhiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=xfr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaaCyEaiaacIcacaWGHbGaaiykaaaa@392E@ and y ( b ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLnhiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=xfr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaaCyEaiaacIcacaWGIbGaaiykaaaa@392F@ and returns the value of the residual for each specified boundary condition. For example to apply y 1 ( a ) = 1 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLnhiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=xfr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamyEamaaBaaaleaacaaIXaaabeaakiaacIcacaWGHbGaaiykaiabg2da9iaaigdaaaa@3BDC@ and y 1 ( b ) = 0 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLnhiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=xfr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamyEamaaBaaaleaacaaIXaaabeaakiaacIcacaWGIbGaaiykaiabg2da9iaaicdaaaa@3BDC@ use

    function res = bcs(ya,yb)
    res = [ ya(1) - 1
            yb(1)];
    

    To apply conditions to the other variables, y 2 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLnhiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=xfr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamyEamaaBaaaleaacaaIYaaabeaaaaa@37D3@ , etc. change the index to ya(2) or yb(2), for example.

  3. Define the solution domain and provide an initial guess for the solution on the solution domain. Use the command

    solinit = bvpinit([a,b],[0,0]);

    This defines the domain for solution as [a,b], and the initial guess for the solution at the points specified in the domain as [0,0]. (Note that we could use a more accurate initial guess, that is define the domain using linspace(a,b,100) and then define the solution on these points.)

  4. Call the ODE solver bvp4c, using the following command
    sol = bvp4c(@deriv,@bcs,solinit);

    The various parameters are:

    • @deriv, a handle to a function that given x MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLnhiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=xfr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamiEaaaa@36EA@ and y MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLnhiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=xfr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaaCyEaaaa@36EF@ returns the value of the derivative d y d x MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLnhiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=xfr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaWaaSaaaeaacaWGKbGaaCyEaaqaaiaadsgacaWG4baaaaaa@39CE@ ;
    • @bcs, a handle to a function that defines the boundary conditions;
    • solinit, the structure defining the solution’s domain and that initial guess at the solution; and
    • sol, a structure that contains the solution.
  5. Plot the results, which are now stored as sol.x and sol.y.
    plot(sol.x,sol.y(1,:),'b-x');

This is shown in the following walkthrough example.

Walkthrough

Suppose we wish to solve the following boundary value problem.

Consider the equation

d 2 y d x 2 + y = 0.

subject to y ' ( 0 ) = 1 and y ( π ) = 0 .

The exact solution is y = sin ( x ) .

To solve this numerically, we first need to reduce the second-order equation to a system of first-order equations,

d y d x = z , d z d x = y .

with z(0)=1 and y(π)=0 .

Example code to solve this is given by

% Function to solve d^2ydx^2+y = 0.
function SimpleBVP()
 
solinit = bvpinit([0,pi],[0,0]);
 
sol = bvp4c(@deriv,@bcs,solinit);
 
plot(sol.x,sol.y(1,:),'b-x');
 
function dYdx = deriv(x,Y)
 
dYdx(1) = Y(2); 
dYdx(2) = -Y(1);
 
% boundary conditions y'(a)=1, y(b)=0.
function res = bcs(ya,yb)
res = [ ya(2) - 1 
        yb(1)];

Download code as SimpleBVP.m file

To run this code, download it to the current working directory and use the following command

SolveBVP

The main elements of this code are

  • solinit = bvpinit([0,pi],[0,0]);
    which defines the domain for solution through [0,π] and the initial guess for the solution at the points specified in the domain, [0,0] .
  • sol = bvp4c(@deriv,@bcs,solinit);
    which is the call to the bvp4c solver, the various parameters are:
    • @deriv, a handle to a function that returns the value of the derivative d y d x MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLnhiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=xfr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaWaaSaaaeaacaWGKbGaamyEaaqaaiaadsgacaWG4baaaaaa@39CA@ for a given x MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLnhiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=xfr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamiEaaaa@36EA@ and y MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLnhiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=xfr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamyEaaaa@36EB@ ;
    • @bcs is a handle to a function that defines the boundary conditions; and
    • solinit is the structure defining the solution domain and initial guess.
  • plot(sol.x,sol.y(1,:),'b-x');
    which plots the solution.
    function dydx = deriv(x,y)
    which is a function that returns the value of d y d x MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLnhiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=xfr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaWaaSaaaeaacaWGKbGaaCyEaaqaaiaadsgacaWG4baaaaaa@39CE@ for a given x MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLnhiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=xfr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamiEaaaa@36EA@ and y MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLnhiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=xfr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaaCyEaaaa@36EF@ at that point (which here is the vector [ z , y ] MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLnhiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=xfr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaai4waiabgkHiTiaadQhacaGGSaGaamyEaiaac2faaaa@3B47@ ).
  • function res = bcs(ya,yb)
    which is a function which defines the boundary conditions applied at a = 0 and b = 1 .

Running the code (using SolveBVP) yields the following plot

This is plotted against the exact solution, y = sin ( x ) , in the next figure.

The red line represents the actual solution and the blue crosses show the numerical solution from bvp4c.

top
Exercises 3 | Go upMATLAB | Exercises 4