Manifesto on Numerical Integration of Equations of

Motion Using Matlab

C. Hall

April 11, 2002

This handout is intended to help you understand numerical integration and to put it
into practice using Matlab’s ode45 function.

1 Integration

We want to integrate ordinary differential equations (ODEs) of the form

ẋ = f (x, t), with initial conditions x(0) = x0 (1)
where the state x may be a scalar or an N × 1 matrix, and t is a scalar. This problem
is known as an initial value problem, the variable t is the independent variable (usually
time), and the state x is the dependent variable. The result of integration is to find the
solution x(t) that satisfies the differential equation and initial conditions given in Eq. (1).

The function f (x, t) is frequently referred to as the right-hand side of the ODE, and
most numerical algorithms involve repetitive calculation of the right-hand side for many
values of x and t. Since f (x, t) is the rate of change of the dependent variable x with
respect to the independent variable t, we can think of the right-hand side as the slope of
x(t), even if x is a multi-dimensional state vector.

2 Some Example Problems of the Form of Eq. (1)

2.1 A simple first-order ODE

One of the simplest possible examples is

ẋ = a, with initial conditions x(0) = b (2)

1
where a and b are constants, and x is a scalar. This equation is easily integrated analyti-
cally to obtain
x(t) = at + b (3)
which defines a straight line in the (t, x) plane, with slope a and intercept b. Notice how
the right-hand side is in fact the slope of x(t), as mentioned in the previous section.

2.2 A slightly less simple first-order ODE

Instead of a constant right-hand side, let the right-hand side have a t-dependence:

ẋ = at, with initial conditions x(0) = b (4)

where a and b are constants, and x is a scalar. This equation is also easily integrated
analytically to obtain
1
x(t) = at2 + b (5)
2
which defines a parabola in the (t, x) plane, with slope at any given point given by at and
intercept b. Notice how the right-hand side is again the slope of x(t), which is changing
as a function of t.

2.3 A more interesting, but still simple first-order ODE

Instead of a t-dependence, let the right-hand side depend on the dependent variable:

ẋ = ax, with initial conditions x(0) = b (6)

where a and b are constants, and x is a scalar. This equation is also easily integrated
analytically by using separation of variables (exercise!) to obtain

x(t) = beat (7)

which defines an exponential curve in the (t, x) plane, with slope at any given point given
by abeat = ax(t) and intercept b. Notice how the right-hand side is still the slope of x(t),
which is again changing as a function of t.

Remark 1 All of the preceding examples are of the class known as linear first-order
ordinary differential equations. The two examples where f (x, t) does not have any t-
dependence are also known as constant-coefficient or time-invariant differential equations.
The example where f (x, t) does have time-dependence is known as a time-varying differ-
ential equation.

2
2.4 Some nonlinear first-order ODEs

If the right-hand side has a nonlinear dependence on the dependent variable, then the
equation is a nonlinear differential equation. Here are some examples:

ẋ = ax2 (8)
ẋ = a sin x (9)
ẋ = aex (10)
ẋ = a log x (11)

All of these examples can be integrated analytically to obtain the solution x(t), because
the equations are all separable. As an exercise, you should separate the equations and
carry out the integration.

2.5 Second-order ODEs

Many problems of interest to engineers arise in the form of second-order ODEs, because
many problems are developed by using Newton’s second law,

f = ma e.g. f = mẍ (12)

where f is the force acting on an object of mass m, and a is the acceleration of the object,
which might be for example the second derivative of a position described by the cartesian
coordinate, x.

A familiar example is the mass-spring-damper system, where a mass m has position

described by a scalar variable x and is acted on by a spring force fs = −kx, by a damper
force fd = −cẋ, and by a control force fc . Application of f = ma to this problem leads
to the second-order ODE
mẍ + cẋ + kx = fc (13)
Assuming that the parameters m, c, and k are all constant, then this equation is a second-
order, constant-coefficient (time-invariant), linear ODE. Furthermore, if the control force
fc is zero, then the equation is called a homogeneous ODE.

As with the simple first-order examples already presented, this equation can be solved
analytically. You should review the analytical solution to this equation in your vibrations
and ODE textbooks. One can also put this second-order ODE into the first-order form of
a system of first-order ODEs as follows. First define two states x1 and x2 as the position
and velocity; i.e.

x1 = x (14)
x2 = ẋ (15)

3
Differentiating these two states leads to
ẋ1 = ẋ = x2 (16)
fc k c
ẋ2 = ẍ = − x − ẋ = u − k̂x1 − ĉx2 (17)
m m m
where we have introduced three new symbols: u = fc /m is the acceleration due to the
control force, and k̂ and ĉ are the spring and damping coefficients divided by the mass.
These equations may be written in matrix form as
" # " #" # " #
ẋ1 0 1 x1 0
= + u (18)
ẋ2 −k̂ −ĉ x2 1
This form of the equation is so common throughout the dynamics and control literature
that it has a standard notation. If we denote the 2 × 1 state vector by x = [x1 x2 ]T , the
2 × 2 matrix that appears in the equation by A, and the 2 × 1 matrix that multiplies u
by B, then we can write the system as
ẋ = Ax + Bu (19)
In this form of linear, constant-coefficient ODEs, the variable matrix x is known as the
state vector, and the variable u is known as the control. The matrices A and B are the
plant and input matrices, respectively.

Since the input must either be constant, a function of time, a function of the state, or
a function of both the time and the state, then this system of equations may be written
in the even more general form of
ẋ = f (x, t) (20)
Note the similarity between this ODE and Eq. (1). The only difference is that we have
introduced the convention of using a bold-face font to denote that a particular variable
or parameter is a matrix.

Also note that the right-hand side is the time-rate of change of the state, which can
be interpreted as the slope of the state as a function of time.

In general, most ordinary differential equations can be written in the form of Eq. (20).
In the next subsection we give some higher-order examples.

2.6 Higher-order ODEs

The rotational equations of motion for a rotating rigid body may be written as the com-
bination of Euler’s equations and a set of kinematics equations. A typical choice for
kinematics variables is the quaternion. These equations are:
h i
ω̇ = −I−1 ω× Iω + g (21)
q̄˙ = Q(q̄)ω (22)

4
Note that the right-hand side of the ω̇ equation depends on ω and the torque g, which
(like the control in Eq. (19)) may depend on the angular velocity ω, on the attitude q̄,
and on the time t. The right-hand side of the q̄˙ equation depends on q̄ and ω. Thus the
equations are coupled. The combination of ω and q̄ forms the state; i.e.
 
ω1
ω2
 
 
 
" #  ω3 
ω  
x= = q1 
(23)
q̄
 
 

 q2 


 q3 

q4

Thus the equations can be combined into a single state vector differential equation (or
system of equations): " # " #
ω̇ −I−1 [ω× Iω + g]
= (24)
q̄˙ Q(q̄)ω
which may be written in the more concise form of

ẋ = f (x, t) (25)

As noted above, most dynamics and control problems can be written in this form, where
the right-hand side corresponds to the rate of change of the state and depends on both
the state and time.

3 Analytical vs.Numerical Solutions

We use the term analytical to refer to a solution that can be written in closed form in
terms of known function such as trigonometric functions, hyperbolic functions, elliptic
functions, and so forth. Most dynamics and control problems do not admit analytical
solutions, and so the solution x(t) corresponding to a given set of initial conditions must
be found numerically using a numerical integration algorithm. A numerical integration
algorithm can be tested by using it to integrate a system of equations that has a known
analytical solution and comparing the results.

4 Numerical Integration Algorithms

There are many popular numerical integration algorithms appropriate for solving the
initial value problems of the form

ẋ = f (x, t) with initial conditions x(0) = x0 (26)

5
The simplest algorithm of any practical use is known as Euler integration. A slightly
more complicated algorithm with substantially better performance is the Runge-Kutta
algorithm. Higher-order (more accurate, but more complicated) algorithms are widely
available, and usually one just selects an existing implementation of a suitable algorithm
and uses it for the problem of interest. In most cases, the algorithm needs to know the
initial and final times (t0 is usually 0, and tf ), the initial conditions (x0 ), and how to
calculate the right-hand side or slope f (x, t) for given values of the state and time.

4.1 Euler integration

Euler integration is based directly on the notion that the right-hand side defines the slope
of the function x(t). Since we know the initial conditions, x0 , the initial slope is simply
f (x0 , t0 ), and one might reasonably suppose that the slope is nearly constant for a short
“enough” timespan following the initial time. Thus, one chooses a small “enough” time
step, ∆t, and estimates the state at time t0 + ∆t as

x1 = x(t0 + ∆t) ≈ x0 + f (x0 , t0 )∆t (27)

Continuing in this manner, one obtains a sequence of approximations at increasing times

as

x1 = x0 + f (x0 , t0 )∆t (28)

x2 = x1 + f (x1 , t1 )∆t (29)
x3 = x2 + f (x2 , t2 )∆t (30)
· (31)
· (32)
· (33)
xn = xn−1 + f (xn−1 , tn−1 )∆t (34)

where tj = tj−1 +∆t. This method is easy to understand, and easy to implement; however
it is not accurate enough for integrating over any reasonably long time intervals. It is
accurate enough to use in short time integrations though.

4.2 Runge-Kutta integration

5 Numerical Integration Using Matlab

Matlab has several numerical integration algorithms implemented as .m files. The most
commonly used of these algorithms is ode45, and some examples using this function are
provided below. In most implementations, the use of the integration function requires

6
the development of two Matlab .m files: a driver file and a right-hand side file. For the
examples here, we use the filenames driver.m and rhs.m. The driver file is a script file,
meaning that it contains a sequence of Matlab statements that are executed in order. The
rhs file is a function file, meaning that it defines a function with arguments (generally t
and x) that returns the value of the right-hand side (or slope) at the given values of time
and state.

The driver file is responsible for defining parameters and initial conditions, calling
the integration function, and then completing any post-processing of the results, such as
making plots of the state variables as functions of time.

First we present a simple implementation of the Runge-Kutta algorithm.

5.1 An implementation of Runge-Kutta

Here is a sample integration algorithm that works similarly to the ode45 algorithm, but
is simpler. It implements the standard Runge-Kutta algorithm. It is a Matlab function,
and has 3 arguments: frhs is a string containing the name of the .m file that computes
the right-hand side of the differential equations; tspan is a 1 × N matrix of the times at
which the algorithm will sample the right-hand side (t0 , t0 + ∆t, ...tf ); and x0 is the initial
state vector.

% oderk.m
% C. Hall
% [t, x] = oderk(frhs,tspan,x0)
% implements runge-kutta 4th order algorithm
% frhs = ’filename’ where filename.m is m-file
% tspan = 1xN matrix of t values
% x0 = nx1 matrix of initial conditions
% returns t=tspan, x=nxN matrix
% x(:,i) = x(t(i))
function [tspan,x] = oderk(frhs,tspan,x0)
N=length(tspan);
n=length(x0);
x0=reshape(x0,n,1);
x=[x0 zeros(n,N-1)];
w=x0;
for i=1:N-1
h=tspan(i+1)-tspan(i);
t=tspan(i);
K1=h*feval(frhs,t,w);
K2=h*feval(frhs,t+h/2,w+K1/2);
K3=h*feval(frhs,t+h/2,w+K2/2);

7
 K4=h*feval(frhs,t+h,w+K3);
w=w+(K1+2*K2+2*K3+K4)/6;
x(:,i+1)=w;
end
x=x’;

As should be evident from the function, it divides each time interval into 2 segments
and computes the right-hand side (slope) at the initial, final, and midpoint times of each
interval. You should make an effort to know how these algorithms work, but you can use
them as black boxes just as you might use an eigenvalue solver. That is, you can think of
oderk as a function that takes its three arguments and returns the state vector calculated
at each of the times in the tspan argument.

Here is an example illustrating the integration of a single first-order equation, ẋ = ax,

x(0) = x0 , from t = 0 to t = 10 with a ∆t of 0.1 s. The driver file is presented first.

% driver.m
% this file sets the parameters and initial conditions,
% calls oderk, and makes a plot of x(t)
global a % make the parameter a global
a = 0.1; % set the parameter a
x0 = 1.0; % initial condition of state
t = 0:0.1:10; % [0 0.1 0.2 ... 10]
fname = ’rhs’; % the name of the rhs file (.m is assumed)
[t,x] = oderk(fname,t,x0); % turn the work over to oderk
plot(t,x) % make the plot
xlabel(’t’); % add some labels
ylabel(’x’);
title(’Sample Runge-Kutta Solution of First-Order ODE’);

The right-hand side file is quite simple.

% rhs.m
% xdot = rhs(t,x)
function xdot = rhs(t,x)
global a
xdot = a*x;
% that’s all!

The same integrator (oderk) can be used to integrate any system of ODEs. For
example, suppose we have the system of three first-order equations:

ẋ = σ(y − x) (35)

8
 ẏ = −xz + rx − y (36)
ż = xy − bz (37)

These equations are called the Lorenz equations and are considered to be the first widely
known system of equations admitting chaotic solutions1 . The behavior of solutions to
these equations is strongly dependent on the values of the three parameters b, σ, and
r. An interesting chaotic attractor (also known as a strange attractor ) occurs for the
parameter values b = 8/3, σ = 10, r = 28.

Here are driver and rhs files to implement these equations and compute the strange
attractor. Try it out.

% lorenzdriver.m
% this file sets the parameters and initial conditions,
% calls oderk, and makes plots
global b sig r % make the parameters global
b = 8/3; % set the parameters
sig = 10;
r = 28;
x0 = [0; 1; 0]; % initial conditions of state vector
t = 0:0.01:50; % [0 0.01 0.02 ... 50]
fname = ’lorenzrhs’; % the name of the rhs file (.m is assumed)
[t,x] = oderk(fname,t,x0); % turn the work over to oderk
plot3(x(:,1),x(:,2),x(:,3));% make a 3-d plot (cool!)
xlabel(’x’); % add some labels
ylabel(’y’);
zlabel(’z’);
title(’Sample Runge-Kutta Solution of Lorenz System’);

The right-hand side file is quite simple.

% lorenzrhs.m
% xdot = rhs(t,x)
function xdot = rhs(t,x)
global b sig r
xx = x(1); % the x component
y = x(2); % and y ...
z = x(3);
xdot = [ sig*(y-xx); -xx*z+r*xx-y; xx*y-b*z];
% that’s all!
1
For more information on chaotic systems, see Gleick’s book Chaos, Moon’s Chaotic and Fractal
Dynamics, or Tabor’s Chaos and Integrability in Nonlinear Dynamics.

9
5.2 Using ode45

The calling format for using ode45 is identical to that for calling oderk. However, there
are numerous additional arguments that one can use in calling ode45 that I have not
included in my implementation of oderk. You can find out what the options are by
typing help ode45 at the Matlab prompt. Probably the most useful of these additional
arguments are the options and parameters arguments. For example, instead of making
the parameters in the Lorenz equations example global variables, we can pass them to
the right-hand side file as parameters. And we can also set the desired precision we want
from the integration. The following driver and rhs files implement these two options.

% lorenzdriverp.m
% this file sets the parameters and initial conditions,
% calls ode45, and makes plots

b = 8/3; % set the parameters

sig = 10;
r = 28;
x0 = [0; 1; 0]; % initial conditions of state vector
t = 0:0.01:50; % [0 0.01 0.02 ... 50]
fname = ’lorenzrhs’; % the name of the rhs file (.m is assumed)
options = odeset(’abstol’,1e-9,’reltol’,1e-9); % set tolerances
[t,x] = ode45(fname,t,x0,options,b,sig,r);
plot3(x(:,1),x(:,2),x(:,3)); % make the plot
xlabel(’x’); % add some labels
ylabel(’y’);
zlabel(’z’);
title(’Sample ode45 Solution of Lorenz System’);

% lorenzrhsp.m
% xdot = rhs(t,x,b,sig,r)
function xdot = rhs(t,x,b,sig,r)
xx = x(1); % the x component
y = x(2); % and y ...
z = x(3);
xdot = [ sig*(y-xx); -xx*z+r*xx-y; xx*y-b*z];
% that’s all!

10