Solve nonstiff differential equations — variable order method
tspan = [t0 tf], integrates the system of
differential equations from
y0. Each row in the solution
y corresponds to a value returned in column
All MATLAB® ODE solvers can solve systems of equations of
the form ,
or problems that involve a mass matrix, .
The solvers all use similar syntaxes. The
only can solve problems with a mass matrix if the mass matrix is constant.
solve problems with a mass matrix that is singular, known as differential-algebraic
equations (DAEs). Specify the mass matrix using the
finds where functions of (t,y),
called event functions, are zero. In the output,
the time of the event,
ye is the solution at the
time of the event, and
ie is the index of the triggered
For each event function, specify whether the integration is
to terminate at a zero and whether the direction of the zero crossing
matters. Do this by setting the
to a function, such as
and creating a corresponding function: [
For more information, see ODE Event Location.
a structure that you can use with
sol = ode113(___)
deval to evaluate
the solution at any point on the interval
You can use any of the input argument combinations in previous syntaxes.
Simple ODEs that have a single solution component can be specified as an anonymous function in the call to the solver. The anonymous function must accept two inputs
(t,y), even if one of the inputs is not used in the function.
Solve the ODE
Specify a time interval of
[0 5] and the initial condition
y0 = 0.
tspan = [0 5]; y0 = 0; [t,y] = ode113(@(t,y) 2*t, tspan, y0);
Plot the solution.
The van der Pol equation is a second-order ODE
where is a scalar parameter. Rewrite this equation as a system of first-order ODEs by making the substitution . The resulting system of first-order ODEs is
The function file
vdp1.m represents the van der Pol equation using . The variables and are the entries
y(2) of a two-element vector
function dydt = vdp1(t,y) %VDP1 Evaluate the van der Pol ODEs for mu = 1 % % See also ODE113, ODE23, ODE45. % Jacek Kierzenka and Lawrence F. Shampine % Copyright 1984-2014 The MathWorks, Inc. dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];
Solve the ODE using the
ode113 function on the time interval
[0 20] with initial values
[2 0]. 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 corresponds to .
[t,y] = ode113(@vdp1,[0 20],[2; 0]);
Plot the solutions for and against
plot(t,y(:,1),'-o',t,y(:,2),'-o') title('Solution of van der Pol Equation (\mu = 1) with ODE113'); xlabel('Time t'); ylabel('Solution y'); legend('y_1','y_2')
ode113 only works with functions that use two input arguments,
y. However, you can pass in extra parameters by defining them outside the function and passing them in when you specify the function handle.
Solve the ODE
Rewriting the equation as a first-order system yields
odefcn.m represents this system of equations as a function that accepts four input arguments:
function dydt = odefcn(t,y,A,B) dydt = zeros(2,1); dydt(1) = y(2); dydt(2) = (A/B)*t.*y(1);
Solve the ODE using
ode113. Specify the function handle such that it passes in the predefined values for
A = 1; B = 2; tspan = [0 5]; y0 = [0 0.01]; [t,y] = ode113(@(t,y) odefcn(t,y,A,B), tspan, y0);
Plot the results.
ode89 solvers are better at solving problems with stringent error tolerances. A common situation where
ode113 excels is in orbital dynamics problems, where the solution curve is smooth and requires high accuracy.
The two-body problem considers two interacting masses
m2 orbiting in a common plane. In this example, one of the masses is significantly larger than the other. With the heavy body at the origin, the equations of motion are
To solve the problem, first convert to a system of four first-order ODEs using the substitutions
The substitutions produce the first-order system
twobodyode codes the system of equations for the two-body problem.
function dy = twobodyode(t,y) % Two body problem with one mass much larger than the other. r = sqrt(y(1)^2 + y(3)^2); dy = [y(2); -y(1)/r^3; y(4); -y(3)/r^3];
twobodyode.m in your working directory, then solve the ODE using
ode113. Specify stringent error tolerances of
opts = odeset('Reltol',1e-13,'AbsTol',1e-14,'Stats','on'); tspan = [0 10*pi]; y0 = [2 0 0 0.5]; [t,y] = ode113(@twobodyode, tspan, y0, opts); plot(t,y) legend('x','x''','y','y''','Location','SouthEast') title('Position and Velocity Components')
924 successful steps 4 failed attempts 1853 function evaluations
figure plot(y(:,1),y(:,3),'-o',0,0,'ro') axis equal title('Orbit of Smaller Mass')
ode113 solver is able to obtain the solution faster and with fewer function evaluations.
ode113 is a variable-step, variable-order (VSVO) Adams-Bashforth-Moulton
PECE solver of orders 1 to 13. The highest order used appears to be 12, however, a
formula of order 13 is used to form the error estimate and the function does local
extrapolation to advance the integration at order 13.
ode113 may be more efficient than
stringent tolerances or if the ODE function is particularly expensive to evaluate.
ode113 is a multistep solver — it normally needs the
solutions at several preceding time points to compute the current solution , .
 Shampine, L. F. and M. K. Gordon, Computer Solution of Ordinary Differential Equations: the Initial Value Problem, W. H. Freeman, San Francisco, 1975.
 Shampine, L. F. and M. W. Reichelt, “The MATLAB ODE Suite,” SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1–22.