Mathematics Previous page   Next Page

Examples: Applying the ODE Initial Value Problem Solvers

This section contains several examples that illustrate the kinds of problems you can solve. For each example, there is a corresponding M-file, included in MATLAB. You can

This section presents the following examples:

Example: Simple Nonstiff Problem

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 on the example name, or type rigidode at the command line.

Example: Stiff Problem (van der Pol Equation)

vdpode illustrates the solution of the van der Pol problem described in Example: The van der Pol Equation, µ = 1000 (Stiff). The differential equations

involve a constant parameter mu.

As mu increases, the problem becomes more stiff, and the period of oscillation becomes larger. When mu is 1000 the equation is in relaxation oscillation and the problem is very stiff. The limit cycle has portions where the solution components change slowly and the problem is quite stiff, alternating with regions of very sharp change where it is not stiff (quasi-discontinuities).

By default, the solvers in the ODE suite that are intended for stiff problems approximate Jacobian matrices numerically. However, this example provides a nested function J(t,y) to evaluate the Jacobian matrix partial derivative of f with respect to y analytically at (t,y) for mu = MU. The use of an analytic Jacobian can improve the reliability and efficiency of integration.

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

Example: 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 u (0, x) = sine of x and boundary conditions u(t,0) = u(t, pi) = 0. An integer n is chosen, h is defined as pi divided by (N + 1), and the solution of the partial differential equation is approximated at x sub k = k times h for k = 0, 1, ..., N+1 by

Here phi sub k times x is a piecewise linear function that is 1 at x sub k and 0 at all the other x sub j. A Galerkin discretization leads to the system of ODEs

and the tridiagonal matrices M(t) and J are given by


The initial values c(0) are taken from the initial condition for the partial differential equation. The problem is solved on the time interval [0, pi].

In the fem1ode example, the properties

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

To run this example, click on 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.

Example: 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 alpha = 1/50 and

There are 2N equations in the system, but the Jacobian is banded with a constant width 5 if the equations are ordered as u sub 1, v sub 1, u sub 2, v sub 2, ...

In the call brussode(N), where N corresponds to N, 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 partial derivative of f with respect to y. 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 n as an argument to brussode. The default is n = 20.

Example: Simple 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

which returns

Detect zero crossings in the negative direction only
Detect all zero crossings
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 (Example: 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.

Example: 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 shown:


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 mu 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.

Example: Differential-Algebraic Problem

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.

In hb1ode, the problem is solved with initial conditions y sub 1 (0) = 1, y sub 2 (0) = 0, y sub 3 (0) = 0 to steady state. These differential equations satisfy a linear conservation law that is used to reformulate the problem as the DAE

Obviously these equations do not have a solution for y(0) with components that do not sum to 1. The problem has the form of M times y prime = f(t,y) with

M is obviously 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 y sub 3 (0) = 10 to the minus 3 power 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:

Summary of Code Examples

The following table lists the M-files for all the ODE initial value problem examples. Click the example name to see the code in an editor. Type the example name at the command line to run it.

Stiff DAE -- electrical circuit
Simple event location -- bouncing ball
ODE with time- and state-dependent mass matrix -- motion of a baton
Stiff large problem -- diffusion in a chemical reaction (the Brusselator)
ODE with strongly state-dependent mass matrix -- Burger's equation solved using a moving mesh technique
Stiff problem with a time-dependent mass matrix -- finite element method
Stiff problem with a constant mass matrix -- finite element method
Stiff ODE problem solved on a very long interval -- Robertson chemical reaction
Robertson problem -- stiff, linearly implicit DAE from a conservation law
Robertson problem -- stiff, fully implicit DAE
Burgers' equation solved as implicit ODE system
Advanced event location -- restricted three body problem
Nonstiff problem -- Euler equations of a rigid body without external forces
Parameterizable van der Pol equation (stiff for large )

Previous page  Changing ODE Integration Properties Questions and Answers, and Troubleshooting Next page

© 1994-2005 The MathWorks, Inc.