ODEs

Function Summary

• ODE Solvers

• Evaluation and Extension

• Solver Options

• Output Functions

ODE Solvers

The following table lists the initial value problem solvers, the kind of problem you can solve with each solver, and the method each solver uses.

Evaluation and Extension

You can use the following functions to evaluate and extend solutions to ODEs.

Solver Options

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

Output Functions

If an output function is specified, the solver calls the specified function after every successful integration step. You can use odeset to specify one of these sample functions as the OutputFcn property, or you can modify them to create your own functions.

Back to Top

Initial Value Problems

• First Order ODEs

• Higher Order ODEs

• Initial Values

First Order ODEs

An ordinary differential equation (ODE) contains one or more derivatives of a dependent variable y with respect to a single independent variable t, usually referred to as time. The derivative of y with respect to t is denoted as y ′, the second derivative as y ′′, and so on. Often y(t) is a vector, having elements y₁, y₂, ..., y_n.

MATLAB solvers handle the following types of first-order ODEs:

• Explicit ODEs of the form y ′ = f (t, y)

• Linearly implicit ODEs of the form M(t, y) y ′ = f (t, y), where M(t, y) is a matrix

• Fully implicit ODEs of the form f (t, y, y ′) = 0 (ode15i only)

Higher Order ODEs

MATLAB ODE solvers accept only first-order differential equations. To use the solvers with higher-order ODEs, you must rewrite each equation as an equivalent system of first-order differential equations of the form

You can write any ordinary differential equation

as a system of first-order equations by making the substitutions

The result is an equivalent system of n first-order ODEs.

Rewrite the second-order van der Pol equation

as a system of first-order ODEs.

Initial Values

Generally there are many functions y(t) that satisfy a given ODE, and additional information is necessary to specify the solution of interest. In an initial value problem, the solution of interest satisfies a specific initial condition, that is, y is equal to y₀ at a given initial time t₀. An initial value problem for an ODE is then

If the function f (t, y) is sufficiently smooth, this problem has one and only one solution. Generally there is no analytic expression for the solution, so it is necessary to approximate y(t) by numerical means.

Back to Top

Types of Solvers

• Nonstiff Problems

• Stiff Problems

• Fully Implicit ODEs

Nonstiff Problems

There are three solvers designed for nonstiff problems:

Stiff Problems

Not all difficult problems are stiff, but all stiff problems are difficult for solvers not specifically designed for them. Solvers for stiff problems can be used exactly like the other solvers. However, you can often significantly improve the efficiency of these solvers by providing them with additional information about the problem. (See Integrator Options.)

There are four solvers designed for stiff problems:

Fully Implicit ODEs

The solver ode15i solves fully implicit differential equations of the form

using the variable order BDF method. The basic syntax for ode15i is

[t,y] = ode15i(odefun,tspan,y0,yp0,options)

The input arguments are

The output arguments contain the solution approximated at discrete points:

See the ode15i reference page for more information about these arguments.

Back to Top

Solver Syntax

All of the ODE solver functions, except for ode15i, share a syntax that makes it easy to try any of the different numerical methods, if it is not apparent which is the most appropriate. To apply a different method to the same problem, simply change the ODE solver function name. The simplest syntax, common to all the solver functions, is

[t,y] = solver(odefun,tspan,y0,options)

where solver is one of the ODE solver functions listed previously.

The basic input arguments are

dydt = odefun(t,y)

The output arguments contain the solution approximated at discrete points:

See the reference page for the ODE solvers for more information about these arguments.

Back to Top

Integrator Options

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

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

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

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

2. Pass the options structure to the solver as follows:

   • For all solvers except ode15i, use the syntax

     [t,y] = solver(odefun,tspan,y0,options)
     

   • For ode15i, use the syntax

     [t,y] = ode15i(odefun,tspan,y0,yp0,options)

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

Back to Top

Examples

• van der Pol Equation (Nonstiff)

• van der Pol Equation (Stiff)

• van der Pol Equation (Parametrizing the ODE)

• van der Pol Equation (Evaluating the Solution)

• Euler Equations (Nonstiff)

• Fully Implicit ODE

• Finite Element Discretization

• Large Stiff Sparse Problem

• Event Location

• Advanced Event Location

• Differential-Algebraic Equations

• Nonnegative Solutions

• Additional Examples

van der Pol Equation (Nonstiff)

This example illustrates the steps for solving an initial value ODE problem:

1. Rewrite the problem as a system of first-order ODEs. Rewrite the van der Pol equation (second-order)

   

   where is a scalar parameter, by making the substitution . The resulting system of first-order ODEs is

   

2. Code the system of first-order ODEs. Once you represent the equation as a system of first-order ODEs, you can code it as a function that an ODE solver can use. The function must be of the form

   dydt = odefun(t,y)
   

   Although t and y must be the function's two arguments, the function does not need to use them. The output dydt, a column vector, is the derivative of y.

   The code below represents the van der Pol system in the function, vdp1. The vdp1 function assumes that . The variables and are the entries y(1) and y(2) of a two-element vector.

   function dydt = vdp1(t,y)
   dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];
   

   Note that, although vdp1 must accept the arguments t and y, it does not use t in its computations.

3. Apply a solver to the problem.

   Decide which solver you want to use to solve the problem. Then call the solver and pass it the function you created to describe the first-order system of ODEs, the time interval on which you want to solve the problem, and an initial condition vector.

   For the van der Pol system, you can use ode45 on time interval [0 20] with initial values y(1) = 2 and y(2) = 0.

   [t,y] = ode45(@vdp1,[0 20],[2; 0]);
   

   This example uses @ to pass vdp1 as a function handle to ode45. The resulting output is a column vector of time points t and a solution array y. Each row in y corresponds to a time returned in the corresponding row of t. The first column of y corresponds to , and the second column to .

   Note   For information on function handles, see the function_handle (@), func2str, and str2func reference pages, and the Function Handles section in MATLAB Programming Fundamentals.

4. View the solver output. You can simply use the plot command to view the solver output.

   plot(t,y(:,1),'-',t,y(:,2),'--')
   title('Solution of van der Pol Equation, \mu = 1');
   xlabel('time t');
   ylabel('solution y');
   legend('y_1','y_2')
   

   

As an alternative, you can use a solver output function to process the output. The solver calls the function specified in the integration property OutputFcn after each successful time step. Use odeset to set OutputFcn to the desired function. See Solver Output Properties, in the reference page for odeset, for more information about OutputFcn.

van der Pol Equation (Stiff)

This example presents a stiff problem. For a stiff problem, solutions can change on a time scale that is very short compared to the interval of integration, but the solution of interest changes on a much longer time scale. Methods not designed for stiff problems are ineffective on intervals where the solution changes slowly because they use time steps small enough to resolve the fastest possible change.

When is increased to 1000, the solution to the van der Pol equation changes dramatically and exhibits oscillation on a much longer time scale. Approximating the solution of the initial value problem becomes a more difficult task. Because this particular problem is stiff, a solver intended for nonstiff problems, such as ode45, is too inefficient to be practical. A solver such as ode15s is intended for such stiff problems.

The vdp1000 function evaluates the van der Pol system from the previous example, but with = 1000.

function dydt = vdp1000(t,y)
dydt = [y(2); 1000*(1-y(1)^2)*y(2)-y(1)];

Note   This example hardcodes in the ODE function. The vdpode example solves the same problem, but passes a user-specified as a parameter to the ODE function.

Now use the ode15s function to solve the problem with the initial condition vector of [2; 0], but a time interval of [0 3000]. For scaling reasons, plot just the first component of y(t).

[t,y] = ode15s(@vdp1000,[0 3000],[2; 0]);
plot(t,y(:,1),'-');
title('Solution of van der Pol Equation, \mu = 1000');
xlabel('time t');
ylabel('solution y_1');

van der Pol Equation (Parametrizing the ODE)

The preceding sections showed how to solve the van der Pol equation for two different values of the parameter µ. In those examples, the values µ = 1 and µ=1000 are hard-coded in the ODE functions. If you are solving an ODE for several different parameter values, it might be more convenient to include the parameter in the ODE function and assign a value to the parameter each time you run the ODE solver. This section explains how to do this for the van der Pol equation.

One way to provide parameter values to the ODE function is to write an M-file that

• Accepts the parameters as inputs.

• Contains ODE function as a nested function, internally using the input parameters.

• Calls the ODE solver.

The following code illustrates this:

function [t,y] = solve_vdp(mu)
tspan = [0 max(20, 3*mu)];
y0 = [2; 0];

% Call the ODE solver ode15s.
[t,y] = ode15s(@vdp,tspan,y0);

    % Define the ODE function as nested function, 
    % using the parameter mu.
    function dydt = vdp(t,y)
    dydt = [y(2); mu*(1-y(1)^2)*y(2)-y(1)];
    end
end

Because the ODE function vdp is a nested function, the value of the parameter mu is available to it.

To run the M-file for mu = 1, enter

[t,y] = solve_vdp(1);

To run the code for µ = 1000, enter

[t,y] = solve_vdp(1000);

See the vdpode code for a complete example based on these functions.

van der Pol Equation (Evaluating the Solution)

The numerical methods implemented in the ODE solvers produce a continuous solution over the interval of integration . You can evaluate the approximate solution, , at any point in using the function deval and the structure sol returned by the solver. For example, if you solve the problem described in van der Pol Equation (Nonstiff) by calling ode45 with a single output argument sol,

sol = ode45(@vdp1,[0 20],[2; 0]);

ode45 returns the solution as a structure. You can then evaluate the approximate solution at points in the vector xint = 1:5 as follows:

xint = 1:5;
Sxint = deval(sol,xint)

Sxint =

    1.5081    0.3235   -1.8686   -1.7407   -0.8344
   -0.7803   -1.8320   -1.0220    0.6260    1.3095

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

Euler Equations (Nonstiff)

rigidode illustrates the solution of a standard test problem proposed by Krogh for solvers intended for nonstiff problems [8].

The ODEs are the Euler equations of a rigid body without external forces.

For your convenience, the entire problem is defined and solved in a single M-file. The differential equations are coded as a subfunction f. Because the example calls the ode45 solver without output arguments, the solver uses the default output function odeplot to plot the solution components.

To run this example, click the example name, or type rigidode at the command line.

function rigidode 
%RIGIDODE  Euler equations: rigid body without external forces
tspan = [0 12];
y0 = [0; 1; 1];

% Solve the problem using ode45
ode45(@f,tspan,y0);
% ------------------------------------------------------------
function dydt = f(t,y)
dydt = [ y(2)*y(3) 
        -y(1)*y(3) 
        -0.51*y(1)*y(2) ];

Fully Implicit ODE

The following example shows how to use the function ode15i to solve the implicit ODE problem defined by Weissinger's equation

with the initial value . The exact solution of the ODE is

The example uses the function weissinger, which is provided with MATLAB, to compute the left-hand side of the equation.

Before calling ode15i, the example uses a helper function decic to compute a consistent initial value for . In the following call, the given initial value is held fixed and a guess of 0 is specified for . See the reference page for decic for more information.

t0 = 1;
y0 = sqrt(3/2);
yp0 = 0;
[y0,yp0] = decic(@weissinger,t0,y0,1,yp0,0);

You can now call ode15i to solve the ODE and then plot the numerical solution against the analytical solution with the following commands.

[t,y] = ode15i(@weissinger,[1 10],y0,yp0);
ytrue = sqrt(t.^2 + 0.5);
plot(t,y,t,ytrue,'o');

Finite Element Discretization

fem1ode illustrates the solution of ODEs that result from a finite element discretization of a partial differential equation. The value of N in the call fem1ode(N) controls the discretization, and the resulting system consists of N equations. By default, N is 19.

This example involves a mass matrix. The system of ODEs comes from a method of lines solution of the partial differential equation

with initial condition and boundary conditions . An integer is chosen, is defined as , and the solution of the partial differential equation is approximated at for k = 0, 1, ..., N+1 by

Here is a piecewise linear function that is 1 at and 0 at all the other . A Galerkin discretization leads to the system of ODEs

and the tridiagonal matrices and are given by

and

The initial values are taken from the initial condition for the partial differential equation. The problem is solved on the time interval .

In the fem1ode example, the properties

options = odeset('Mass',@mass,'MStateDep','none','Jacobian',J) 

indicate that the problem is of the form . The nested function mass(t) evaluates the time-dependent mass matrix and J is the constant Jacobian.

To run this example, click the example name, or type fem1ode at the command line. From the command line, you can specify a value of N as an argument to fem1ode. The default is N = 19.

function fem1ode(N)
%FEM1ODE  Stiff problem with a time-dependent mass matrix 

if nargin < 1
  N = 19;
end
h = pi/(N+1);
y0 = sin(h*(1:N)');
tspan = [0; pi];

% The Jacobian is constant.
e = repmat(1/h,N,1);    %  e=[(1/h) ... (1/h)];
d = repmat(-2/h,N,1);   %  d=[(-2/h) ... (-2/h)]; 
% J is shared with the derivative function.
J = spdiags([e d e], -1:1, N, N);

d = repmat(h/6,N,1);  
% M is shared with the mass matrix function.
M = spdiags([d 4*d d], -1:1, N, N);

options = odeset('Mass',@mass,'MStateDep','none', ...
                 'Jacobian',J);

[t,y] = ode15s(@f,tspan,y0,options);

figure;
surf((1:N)/(N+1),t,y);
set(gca,'ZLim',[0 1]);
view(142.5,30);
title(['Finite element problem with time-dependent mass ' ...
       'matrix, solved by ODE15S']);
xlabel('space ( x/\pi )');
ylabel('time');
zlabel('solution');
%--------------------------------------------------------------
function yp = f(t,y)
% Derivative function.
   yp = J*y;    % Constant Jacobian provided by outer function
end             % End nested function f
%--------------------------------------------------------------
function Mt = mass(t)
% Mass matrix function.
   Mt = exp(-t)*M;    % M is provided by outer function
end                   % End nested function mass
%--------------------------------------------------------------
end

Large Stiff Sparse Problem

brussode illustrates the solution of a potentially large stiff sparse problem. The problem is the classic "Brusselator" system [3] that models diffusion in a chemical reaction

and is solved on the time interval [0,10] with = 1/50 and

There are equations in the system, but the Jacobian is banded with a constant width 5 if the equations are ordered as

In the call brussode(N), where N corresponds to , the parameter N ≥ 2 specifies the number of grid points. The resulting system consists of 2N equations. By default, N is 20. The problem becomes increasingly stiff and the Jacobian increasingly sparse as N increases.

The nested function f(t,y) returns the derivatives vector for the Brusselator problem. The subfunction jpattern(N) returns a sparse matrix of 1s and 0s showing the locations of nonzeros in the Jacobian . The example assigns this matrix to the property JPattern, and the solver uses the sparsity pattern to generate the Jacobian numerically as a sparse matrix. Providing a sparsity pattern can significantly reduce the number of function evaluations required to generate the Jacobian and can accelerate integration.

For the Brusselator problem, if the sparsity pattern is not supplied, 2N evaluations of the function are needed to compute the 2N-by-2N Jacobian matrix. If the sparsity pattern is supplied, only four evaluations are needed, regardless of the value of N.

To run this example, click on the example name, or type brussode at the command line. From the command line, you can specify a value of as an argument to brussode. The default is  = 20.

function brussode(N)
%BRUSSODE  Stiff problem modeling a chemical reaction 

if nargin < 1
  N = 20;
end

tspan = [0; 10];
y0 = [1+sin((2*pi/(N+1))*(1:N));
repmat(3,1,N)];

options = odeset('Vectorized','on','JPattern',jpattern(N));

[t,y] = ode15s(@f,tspan,y0,options);

u = y(:,1:2:end);
x = (1:N)/(N+1);
surf(x,t,u);
view(-40,30);
xlabel('space');
ylabel('time');
zlabel('solution u');
title(['The Brusselator for N = ' num2str(N)]);
% --------------------------------------------------------------
function dydt = f(t,y)
c = 0.02 * (N+1)^2;
dydt = zeros(2*N,size(y,2));      % preallocate dy/dt
% Evaluate the two components of the function at one edge of 
% the grid (with edge conditions).
i = 1;
dydt(i,:) = 1 + y(i+1,:).*y(i,:).^2 - 4*y(i,:) + ...
            c*(1-2*y(i,:)+y(i+2,:));
dydt(i+1,:) = 3*y(i,:) - y(i+1,:).*y(i,:).^2 + ...
              c*(3-2*y(i+1,:)+y(i+3,:));
% Evaluate the two components of the function at all interior 
% grid points.
i = 3:2:2*N-3;
dydt(i,:) = 1 + y(i+1,:).*y(i,:).^2 - 4*y(i,:) + ...
            c*(y(i-2,:)-2*y(i,:)+y(i+2,:));
dydt(i+1,:) = 3*y(i,:) - y(i+1,:).*y(i,:).^2 + ...
              c*(y(i-1,:)-2*y(i+1,:)+y(i+3,:));
% Evaluate the two components of the function at the other edge 
% of the grid (with edge conditions).
i = 2*N-1;
dydt(i,:) = 1 + y(i+1,:).*y(i,:).^2 - 4*y(i,:) + ...
            c*(y(i-2,:)-2*y(i,:)+1);
dydt(i+1,:) = 3*y(i,:) - y(i+1,:).*y(i,:).^2 + ...
              c*(y(i-1,:)-2*y(i+1,:)+3);
end   % End nested function f
end   % End function brussode
% --------------------------------------------------------------
function S = jpattern(N)
B = ones(2*N,5);
B(2:2:2*N,2) = zeros(N,1);
B(1:2:2*N-1,4) = zeros(N,1);
S = spdiags(B,-2:2,2*N,2*N);
end;

Event Location

ballode models the motion of a bouncing ball. This example illustrates the event location capabilities of the ODE solvers.

The equations for the bouncing ball are:

In this example, the event function is coded in a subfunction events

[value,isterminal,direction] = events(t,y)

which returns

• A value of the event function

• The information whether or not the integration should stop when value = 0 (isterminal = 1 or 0, respectively)

• The desired directionality of the zero crossings:

-1	Detect zero crossings in the negative direction only
0	Detect all zero crossings
1	Detect zero crossings in the positive direction only

The length of value, isterminal, and direction is the same as the number of event functions. The ith element of each vector, corresponds to the ith event function. For an example of more advanced event location, see orbitode (Advanced Event Location).

In ballode, setting the Events property to @events causes the solver to stop the integration (isterminal = 1) when the ball hits the ground (the height y(1) is 0) during its fall (direction = -1). The example then restarts the integration with initial conditions corresponding to a ball that bounced.

To run this example, click on the example name, or type ballode at the command line.

function ballode
%BALLODE  Run a demo of a bouncing ball.

tstart = 0;
tfinal = 30;
y0 = [0; 20];
refine = 4;
options = odeset('Events',@events,'OutputFcn', @odeplot,...
                 'OutputSel',1,'Refine',refine);

set(gca,'xlim',[0 30],'ylim',[0 25]);
box on
hold on;

tout = tstart;
yout = y0.';
teout = [];
yeout = [];
ieout = [];
for i = 1:10
  % Solve until the first terminal event.
  [t,y,te,ye,ie] = ode23(@f,[tstart tfinal],y0,options);
  if ~ishold
    hold on
  end
  % Accumulate output.  
  nt = length(t);
  tout = [tout; t(2:nt)];
  yout = [yout; y(2:nt,:)];
  teout = [teout; te];    % Events at tstart are never reported.
  yeout = [yeout; ye];
  ieout = [ieout; ie];

  ud = get(gcf,'UserData');
  if ud.stop
    break;
  end
  
  % Set the new initial conditions, with .9 attenuation.
  y0(1) = 0;
  y0(2) = -.9*y(nt,2);

  % A good guess of a valid first time step is the length of 
  % the last valid time step, so use it for faster computation.
  options = odeset(options,'InitialStep',t(nt)-t(nt-refine),...
                           'MaxStep',t(nt)-t(1));
  tstart = t(nt);
end

plot(teout,yeout(:,1),'ro')
xlabel('time');
ylabel('height');
title('Ball trajectory and the events');
hold off
odeplot([],[],'done');
% --------------------------------------------------------------
function dydt = f(t,y)
dydt = [y(2); -9.8];
% --------------------------------------------------------------
function [value,isterminal,direction] = events(t,y)
% Locate the time when height passes through zero in a 
% decreasing direction and stop integration.
value = y(1);     % Detect height = 0
isterminal = 1;   % Stop the integration
direction = -1;   % Negative direction only

Advanced Event Location

orbitode illustrates the solution of a standard test problem for those solvers that are intended for nonstiff problems. It traces the path of a spaceship traveling around the moon and returning to the earth (Shampine and Gordon [8], p. 246).

The orbitode problem is a system of the following four equations:

where

The first two solution components are coordinates of the body of infinitesimal mass, so plotting one against the other gives the orbit of the body. The initial conditions have been chosen to make the orbit periodic. The value of corresponds to a spaceship traveling around the moon and the earth. Moderately stringent tolerances are necessary to reproduce the qualitative behavior of the orbit. Suitable values are 1e-5 for RelTol and 1e-4 for AbsTol.

The nested events function includes event functions that locate the point of maximum distance from the starting point and the time the spaceship returns to the starting point. Note that the events are located accurately, even though the step sizes used by the integrator are not determined by the location of the events. In this example, the ability to specify the direction of the zero crossing is critical. Both the point of return to the initial point and the point of maximum distance have the same event function value, and the direction of the crossing is used to distinguish them.

To run this example, click on the example name, or type orbitode at the command line. The example uses the output function odephas2 to produce the two-dimensional phase plane plot and let you to see the progress of the integration.

function orbitode
%ORBITODE  Restricted three-body problem

mu = 1 / 82.45;
mustar = 1 - mu;
y0 = [1.2; 0; 0; -1.04935750983031990726];
tspan = [0 7];

options = odeset('RelTol',1e-5,'AbsTol',1e-4,...
                 'OutputFcn',@odephas2,'Events',@events);

[t,y,te,ye,ie] = ode45(@f,tspan,y0,options);

plot(y(:,1),y(:,2),ye(:,1),ye(:,2),'o');
title ('Restricted three body problem')
ylabel ('y(t)')
xlabel ('x(t)')
% --------------------------------------------------------------
function dydt = f(t,y)
r13 = ((y(1) + mu)^2 + y(2)^2) ^ 1.5;
r23 = ((y(1) - mustar)^2 + y(2)^2) ^ 1.5;
dydt = [ y(3)
         y(4)
         2*y(4) + y(1) - mustar*((y(1)+mu)/r13) - ...
                         mu*((y(1)-mustar)/r23)
        -2*y(3) + y(2) - mustar*(y(2)/r13) - mu*(y(2)/r23) ];
end   % End nested function f
% --------------------------------------------------------------
function [value,isterminal,direction] = events(t,y)
% Locate the time when the object returns closest to the
% initial point y0 and starts to move away; stop integration.
% Also locate the time when the object is farthest from the 
% initial point y0 and starts to move closer.
% 
% The current distance of the body is
% 
%   DSQ = (y(1)-y0(1))^2 + (y(2)-y0(2))^2 
%       = <y(1:2)-y0(1:2),y(1:2)-y0(1:2)>
% 
% A local minimum of DSQ occurs when d/dt DSQ crosses zero 
% heading in the positive direction. Compute d(DSQ)/dt as
% 
%  d(DSQ)/dt = 2*(y(1:2)-y0(1:2))'*dy(1:2)/dt = ... 
%                 2*(y(1:2)-y0(1:2))'*y(3:4)
% 
dDSQdt = 2 * ((y(1:2)-y0(1:2))' * y(3:4));
value = [dDSQdt; dDSQdt];
isterminal = [1; 0];            % Stop at local minimum
direction = [1; -1];            % [local minimum, local maximum]
end   % End nested function events
end 

Differential-Algebraic Equations

hb1dae reformulates the hb1ode example as a differential-algebraic equation (DAE) problem. The Robertson problem coded in hb1ode is a classic test problem for codes that solve stiff ODEs.

Note   The Robertson problem appears as an example in the prolog to LSODI [4].

In hb1ode, the problem is solved with initial conditions , , to steady state. These differential equations satisfy a linear conservation law that is used to reformulate the problem as the DAE

These equations do not have a solution for with components that do not sum to 1. The problem has the form of with

is singular, but hb1dae does not inform the solver of this. The solver must recognize that the problem is a DAE, not an ODE. Similarly, although consistent initial conditions are obvious, the example uses an inconsistent value to illustrate computation of consistent initial conditions.

To run this example, click on the example name, or type hb1dae at the command line. Note that hb1dae

• Imposes a much smaller absolute error tolerance on than on the other components. This is because is much smaller than the other components and its major change takes place in a relatively short time.

• Specifies additional points at which the solution is computed to more clearly show the behavior of .

• Multiplies by 10⁴ to make visible when plotting it with the rest of the solution.

• Uses a logarithmic scale to plot the solution on the long time interval.

  function hb1dae
  %HB1DAE  Stiff differential-algebraic equation (DAE)
  
  % A constant, singular mass matrix
  M = [1 0 0
       0 1 0 
       0 0 0];
  
  % Use inconsistent initial condition to test initialization.
  y0 = [1; 0; 1e-3];
  tspan = [0 4*logspace(-6,6)];
  
  % Use the LSODI example tolerances.'MassSingular' is left
  % at its default 'maybe' to test the automatic detection
  % of a DAE.
  options = odeset('Mass',M,'RelTol',1e-4,...
                   'AbsTol',[1e-6 1e-10 1e-6],...
                   'Vectorized','on');
  
  [t,y] = ode15s(@f,tspan,y0,options);
  
  y(:,2) = 1e4*y(:,2);
  
  semilogx(t,y);
  ylabel('1e4 * y(:,2)');
  title(['Robertson DAE problem with a Conservation Law, '...
         'solved by ODE15S']);
  xlabel('This is equivalent to the stiff ODEs in HB1ODE.');
  % --------------------------------------------------------
  function out = f(t,y)
  out = [ -0.04*y(1,:) + 1e4*y(2,:).*y(3,:)
           0.04*y(1,:) - 1e4*y(2,:).*y(3,:) - 3e7*y(2,:).^2
           y(1,:) + y(2,:) + y(3,:) - 1 ];
  

  

Nonnegative Solutions

If certain components of the solution must be nonnegative, use odeset to set the NonNegative property for the indices of these components.

Note   This option is not available for ode23s, ode15i, or for implicit solvers (ode15s, ode23t, ode23tb) applied to problems where there is a mass matrix.

Imposing nonnegativity is not always a trivial task. We suggest that you use this option only when necessary, for example in instances in which the application of a solution or integration will fail otherwise.

Consider the following initial value problem solved on the interval [0, 40]:

y' = - |y|, y(0) = 1

The solution of this problem decays to zero. If a solver produces a negative approximate solution, it begins to track the solution of the ODE through this value, the solution goes off to minus infinity, and the computation fails. Using the NonNegative property prevents this from happening.

In this example, the first call to ode45 uses the defaults for the solver parameters:

ode = @(t,y) -abs(y);
[t0,y0] = ode45(ode,[0, 40], 1);

The second uses options to impose nonnegativity conditions:

options = odeset('NonNegative',1);
[t1,y1] = ode45(ode,[0, 40], 1, options);

This plot compares the numerical solution to the exact solution.

Here is a more complete view of the code used to obtain this plot:

ode = @(t,y) -abs(y);
options = odeset('Refine',1);
[t0,y0] = ode45(ode,[0, 40], 1,options);
options = odeset(options,'NonNegative',1);
[t1,y1] = ode45(ode,[0, 40], 1, options);
t = linspace(0,40,1000);
y = exp(-t);
plot(t,y,'b-',t0,y0,'ro',t1,y1,'b*');
legend('Exact solution','No constraints','Nonnegativity', ...
       'Location','SouthWest')

The kneeode Demo.   The MATLAB kneeode demo solves the "knee problem" by imposing a nonnegativity constraint on the numerical solution. The initial value problem is

ε*y' = (1-x)*y - y^2,    y(0) = 1 

For 0 < ε < 1, the solution of this problem approaches null isoclines y = 1 - x and y = 0 for x < 1 and x > 1, respectively. The numerical solution, when computed with default tolerances, follows the y = 1 - x isocline for the whole interval of integration. Imposing nonnegativity constraints results in the correct solution.

Here is the code that makes up the kneeode demo:

function kneeode
%KNEEODE  The "knee problem" with Nonnegativity constraints.

% Problem parameter
epsilon = 1e-6;

y0 = 1;
xspan = [0, 2];
 
% Solve without imposing constraints 
options = [];
[x1,y1] = ode15s(@odefcn,xspan,y0,options);
 
% Impose nonnegativity constraint
options = odeset('NonNegative',1);
[x2,y2] = ode15s(@odefcn,xspan,y0,options);
 
figure
plot(x1,y1,'b.-',x2,y2,'g-')
axis([0,2,-1,1]);
title('The "knee problem"');
legend('No constraints','nonnegativity')
xlabel('x');
ylabel('solution y')

   function yp = odefcn(x,y)
      yp = ((1 - x)*y - y^2)/epsilon;
   end  
end  % kneeode

The derivative function is defined within nested function odefcn. The value of epsilon used in odefcn is obtained from the outer function:

function yp = odefcn(x,y)
yp = ((1 - x)*y - y^2)/epsilon;
end

The demo solves the problem using the ode15s function, first with the default options, and then by imposing a nonnegativity constraint. To run the demo, type kneeode at the MATLAB command prompt.

Here is the output plot. The plot confirms correct solution behavior after imposing constraints.

Additional Examples

The following additional examples are available. Type

edit examplename

to view the code and

examplename

to run the example.

Back to Top

Troubleshooting

• General Questions

• Memory and Computational Efficiency

• Time Steps for Integration

• Error Tolerance and Other Options

• Other Types of Equations

• Other Common Problems

General Questions

Memory and Computational Efficiency

Time Steps for Integration

Error Tolerance and Other Options

Other Types of Equations

Other Common Problems

t0 ≤ t ≤ tf

Back to Top

Differential Equations

DDEs

Recommended Products

• MATLAB
• Signal Processing Toolbox
• Optimization Toolbox
• Statistics Toolbox
• Control System Toolbox

Get MATLAB trial software

Includes the most popular MATLAB recorded presentations with Q&A sessions led by MATLAB experts.

Get the Interactive Kit