> Solving ODEs > Boundary value problems

Boundary value problems

Problem definition

Suppose we wish to solve the system of equations $\frac{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, $\frac{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:

Construct a function (here called deriv) which has input arguments $x$ and $y = (y_{1}, \dots, y_{n})$ and returns the value of the derivative $\frac{d y}{d x}$ , that is $f (x, y)$ .
Construct a function (here called bcs) which has input arguments $y (a)$ and $y (b)$ and returns the value of the residual for each specified boundary condition. For example to apply $y_{1} (a) = 1$ and $y_{1} (b) = 0$ use
```
function res = bcs(ya,yb)
res = [ ya(1) - 1
        yb(1)];
```
To apply conditions to the other variables, $y_{2}$ , etc. change the index to ya(2) or yb(2), for example.
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.)
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$ and $y$ returns the value of the derivative $\frac{d y}{d x}$ ;
- @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.
Plot the results, which are now stored as sol.x and sol.y.
```
plot(sol.x,sol.y(1,:),'b-x');
```

Walkthrough

Suppose we wish to solve the following boundary value problem.

Consider the equation $\frac{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, $\frac{d y}{d x} = z,$ $\frac{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()
    %% Set up, solve, & plot the BVP
    % Initial guess for solution on solution domain
    solinit = bvpinit([0,pi],[0,0]);
    % Solve the BVP
    sol = bvp4c(@deriv,@bcs,solinit);
    % Plot the solution
    plot(sol.x,sol.y(1,:),'b-x');
    %% Function to evaluate the right hand side of the system
    function dYdx = deriv(x,Y)
        dYdx(1) = Y(2);
        dYdx(2) = -Y(1);
    end
    %% Function to evaluate the boundary values, y'(a)=1, y(b)=0
    function res = bcs(ya,yb)
        res = [ ya(2) - 1
            yb(1)];
    end
end

You can run the code by saving in a file named SimpleBVP.m and executing it with SimpleBVP(). 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 $\frac{d y}{d x}$ for a given $x$ and $y$ ;
- @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 $\frac{d y}{d x}$ for a given $x$ and $y$ at that point (which here is the vector $[- z, y]$ ).
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 SimpleBVP()) yields the following plot:

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