EME 171: Analysis, Simulation and Design of Mechatronic Systems

ODE Integration Best Practices With Octave/Matlab

This document describes some recommended best practices for integrating ordinary differential equations using Octave or Matlab. Following these guidelines will result in well organized, modular, readable code and provide advantages in computational efficiency [1].

Information

Extra information (asides) are in blue boxes.

Warnings

Warnings are in yellow boxes.

[1]Efficiency is not the primary purpose of this document. More terse optimized code can certainly be written to maximize computation speed if needed.

Evaluating the ODEs

There are several equations that are of potential interest, but the primary equations are the ordinary differential equations, the "ODEs". These equations always have to be considered; other equations described below are optional. Your primary goal is to solve the ODEs and discover the time dependent behavior of the system's states. You should already have the equations in explicit first order form before moving forward here. "Explicit" refers to the fact that there are no time derivatives on the right hand side of the equations and "first order" refers to their being \(m\) equations, one for each of the \(m\) state variables. The general form for these equations is:

\begin{equation*} \dot{\mathbf{x}}(t) = \mathbf{f}(t, \mathbf{x}(t), \mathbf{r}(t), \mathbf{p}) \end{equation*}

where

  • \(\mathbf{f}\) is the function defining the right hand side of the explicit first order ordinary differential equations. It defines the time derivatives of the states at any given time, i.e. how the states change with time.
  • \(\mathbf{x}(t)\) is the state trajectory vector which is implicitly a function of time, size \(m\times1\).
  • \(\mathbf{r}(t)\) is the input trajectory vector which is implicitly a function of time, size \(o\times1\).
  • \(\mathbf{p}\) is the constant parameter vector (not a function of time), size \(p\times1\).

Note that it is customary to drop the \((t)\) that expresses the implicit function of time for the various time dependent variables. That is done from here forward.

The state trajectory, \(\mathbf{x}\), is determined by integrating \(\mathbf{f}\) with respect to time:

\begin{equation*} \mathbf{x} = \int_{t_0}^{t_f} \mathbf{f}(t, \mathbf{x}, \mathbf{r}, \mathbf{p}) dt \end{equation*}

In general, these equations will be non-linear functions of the state variables and there is unlikely to be analytical symbolic solutions that describe the state trajectory in time. Thus, numerical integration is required to reach a solution. This is fine because linear ODEs really only represent a tiny percentage of all possible ODEs and the methods described here work with linear and nonlinear ODEs.

The following describes how to numerically integrate \(\mathbf{f}\) using Octave or Matlab. The first step is to write an Octave/Matlab function that evaluates \(\mathbf{f}\).

Defining the State Derivative Function

Schematic of a simple pendulum.

For example, the two explicit first order ordinary differential equations of a simple pendulum are:

\begin{equation*} \mathbf{f} = \begin{bmatrix} \dot{\theta} \\ \dot{\omega} \end{bmatrix} = \begin{bmatrix} \omega \\ \frac{-mgl\sin(\theta) + \tau}{ml^2} \end{bmatrix} \end{equation*}
  • The state vector is \(\mathbf{x} = [\theta \quad \omega]^T\), where \(\theta\) is the pendulum angle and \(\omega\) is the angular rate.
  • The parameter vector is \(\mathbf{p} = [m \quad l \quad g]^T\). The variables, respectively, are mass, length, and acceleration due to gravity.
  • The input vector is \(\mathbf{r} = [\tau]^T\), a torque acting between the pendulum and its static attachment (inertial space).

The derivation of these equations can be found on the relevant wikipedia page.

The first step is to translate these ordinary differential equations into a function that evaluates \(\mathbf{f}\) at any given time instant. Below a function named eval_rhs() is defined in an m-file named eval_rhs.m that does so. rhs stands for the "right hand side" of the differential equations.

Note that the inputs and outputs of this function are carefully documented in the lines just below the function signature. Typing help eval_rhs will print this documentation at the Octave/Matlab command prompt.

function xdot = eval_rhs(t, x, r, p)
% EVAL_RHS - Returns the time derivative of the states, i.e. evaluates the
% right hand side of the explicit ordinary differential equations.
%
% Syntax: xdot = eval_rhs(t, x, r, p)
%
% Inputs:
%   t - Scalar value of time, size 1x1.
%   x - State vector at time t, size mx1 where m is the number of states.
%   r - Input vector at time t, size ox1 were o is the number of inputs.
%   p - Constant parameter structure with p items where p is the number of
%       parameters.
% Outputs:
%   xdot - Time derivative of the states at time t, size mx1.

% unpack the states into useful variable names
theta = x(1);
omega = x(2);

% unpack the inputs into useful variable names
tau = r(1);

% unpack the parameters into useful variable names
m = p.m;
l = p.l;
g = p.g;

% calculate the state derivatives
thetadot = omega;
omegadot = (-m*g*l*sin(theta) + tau)/(m*l^2);

% pack the state derivatives into an mx1 vector (same order as states)
xdot = [thetadot; omegadot];

end

Matlab/Octave Structures

Octave/Matlab have a data type called a "structure" (similar to a C structure). The simplest use of a structure is a key-value pair mapping. To create a structure with a single key value pair assign a value to the desired name of the structure followed by a . and the key variable name.

>> p.a = 2
p =

  scalar structure containing the fields:

    a =  2

More key value pairs can be added to the same structure by repeating what the same thing with different keys and values.

>> p.b = 3
p =

  scalar structure containing the fields:

    a =  2
    b =  3

>> p.c = 4.192
p =

  scalar structure containing the fields:

    a =  2
    b =  3
    c =  4.1920

Now you can access these values of the structure p by appending . and the key's variable name. The values in the structure can be used in computations just like other variables.

>> p.a
ans =  2
>> p.a^2 + p.b^2 - p.c^2
ans = -4.5729

You can change the values too:

>> p.c = 5.1
>> p
p =

  scalar structure containing the fields:

    a =  2
    b =  3
    c =  5.1

A structure is a convenient way to store the constant parameters of your system.

Integrating the Equations

Once the function is defined, you can integrate the differential equations with one of the available Octave/Matlab integrators or one of your own design:

%% Script that demonstrates basic integration of ODEs.

% create time values that you desire a solution at, size 1x500
% for example, 0 to 10 seconds with 500 equally spaced time values
ts = linspace(0, 10, 500);

% create a vector with the initial state values, size 2x1
theta0 = 5*pi/180;  % angle in rad
omega0 = 0.0;  % angular rate in rad/s
x0 = [theta0; omega0];  % 2x1 vector

% create a structure to hold all the constants, be careful with units!
p.m = 1.00;  % mass in kg
p.l = 1.00;  % length in m
p.g = 9.81;  % acc due to gravity in m/s^2

% create a vector to hold all of the inputs (constant torque in this case),
% size 1x1
r = [2.0];  % tau, torque in N-m

% check if the eval_rhs function works using an arbitrary time value, the
% initial state values, and the parameters. If it doesn't you'll get an
% error at this line or an unexpected output.
display('Checking eval_rhs:')
eval_rhs(5.0, x0, r, p)

% create an anonymous function with the required inputs for ode45(), i.e.
% (t, x). Note that r and p are set to the values above on creation of this
% function.
% The is needed for two reasons:
% 1. You can only store a m-file function in a variable by using anonymous
%    functions and ode45 requires that the function be anonymous
% 2. ode45 requires that the function only has t and x as arguments, but
%    we'd like to pass in the values of p and r from the variables declared
%    above. This allows us to do that.
% Point 2 above is better than using global variables to share variables in
% all function scopes or declaring these parameters directly in the right
% hand side function where they have limited access.
f_anon = @(t, x) eval_rhs(t, x, r, p);

% integrate the equations with one of the available integrators, in this
% case the Runga-Kutta 4,5 method (good for simple, non-stiff systems).
[ts, xs] = ode45(f_anon, ts, x0);

% the default output of the integrator are the times and states at each
% desired time value:
% - ts, size 500 x 1
% - xs, size 500 x 2

% plotting the states versus time should be your first check to see if the
% result seems reasonable
subplot(211)
plot(ts, xs(:, 1))
ylabel('Angle [rad]')

subplot(212)
plot(ts, xs(:, 2))
ylabel('Angular Rate [rad]')

xlabel('Time [s]')

Only define numbers once!

Note that the constant parameters are only defined in this file. This is on purpose. If you define numerical values redundantly in multiple files and functions you significantly increase your chances of having an erroroneous output due to forgetting to change them all when you make edits.

You may be wondering what the @ symbol specifically means. This designates an anonymous function and is required by ode45(). The following section explains what an anonymous function is along with why and how it can be used.

Anonymous Functions

An anonymous function was used in the above script. The @ symbol indicates this type of function. An anonymous function has three important features that a normal function (written in a unique m-file) doesn't have:

  1. The function can be written in a single line (in fact, if your anonymous function is longer that a single line, 79 characters or so, you should move functionality into a normal function m-file).
  2. The function can be stored in a variable that can be passed to other functions. For example, ode45() requires that the right hand side function be passed in as a variable.
  3. Variables declared in the same scope as and before the anonymous function will be available in the anonymous function. This allows you to avoid the use of global variables or other bad practices at making the values available across a set of functions and scripts.

Anonymous functions are declared with the following syntax:

var_name = @(arg1, arg2, arg3, ...) expression involving the args;

You can use anonymous functions to declare simple functions that fit on one line:

>> my_func = @(x, y) x + y;
>> my_func(1, 2)
ans = 3

use and alternative name for an existing function:

>> my_mean = @mean;
my_mean = @mean
>> my_mean([1, 2, 3])
ans =  2

use anonymous functions to customize the input to existing functions:

>> my_func = @(x, y, z) mean([x, y, z]);
>> my_func(1, 2, 3)
ans = 2

and use anonymous functions to access values stored in variables in the script's scope:

>> b = 2;
>> c = 3;
>> my_func = @(x) mean([x, b, c]);
>> my_func(1)
ans = 2

Note that you have to declare the variables before declaring the anonymous function, the following code fails to compute:

>> clear all;
>> a = 1;
>> my_func = @(x) mean([x, b, c]);
>> my_func(a)
error: 'b' undefined near line 1 column 30
error: called from
    @<anonymous> at line 1 column 22
>> b = 2;
>> c = 3;
>> my_func(a)
error: 'b' undefined near line 1 column 30
error: called from
    @<anonymous> at line 1 column 22

Why not global variables?

It is possible to use global variables to simultaneously make the constant parameters available to both your primary script file and the file that defines your state derivative function. This works, but it is best to avoid global variables except for special needs. Each function provides a unique scope where all variables defined in the function are contained in the function. Using global variables increases the likelihood of programming errors when programs become more complex. A google search on "why global variables are bad" will provide you with background. Here is a Matlab specific note on them:

https://matlab.fandom.com/wiki/FAQ#Are_global_variables_bad.3F

Computation speed of eval_rhs

This function will be executed many times so it is important that this function only calculates the state derivatives and does nothing else. A simple ODE solver will evaluate the function \(n\) times, where \(n\) is the number of time instances you desire a solution at. But any quality ODE solver will execute this function more or less times than \(n\). The solvers are often adaptive and will adjust the time step during integration to ensure low integration error. Fewer time evaluations are needed for slowly changing trajectories and more evaluations are needed when the trajectories change rapidly. Systems that have rapidly changing state trajectories are referred to as "stiff systems" or "stiff equations". For example, a stiff system may require \(1000 \times n\) executions for an acceptable solution. Below, it is shown how to calculate all desired quantities that you may be tempted to calculate in eval_rhs so that you can keep this function minimal.

For example, the number of right hand side function evaluations can be obtained by turning on the stats option for the integrator. Below shows that the equations, as described above, only need to be evaluated about half the number of desired output times.

>> x0 = [5*pi/180; 0];
>> ts = linspace(0, 10, 500);
>> r = [5.0];
>> p = [1; 1; 9.81];
>> f_anon = @(t, x) eval_rhs(t, x, r, p);
>> opt = odeset('stats', 'on');
>> t_start = time();
>> solution = ode45(f_anon, ts, x0, opt);
>> time() - t_start
ans = 0.050903
>> solution.stats.nfevals
ans =  217

But notice that if the system is stiffened, significantly increaasing \(g\) does this, it now takes almost twice the number of evaluations than the desired output times.

>> p = [1; 1; 1000];
>> f_anon = @(t, x) eval_rhs_with_input(t, x, @eval_input, p);
>> t_start = time();
>> solution = ode45(f_anon, ts, x0, opt);
>> time() - t_start
ans = 0.35892
>> solution.stats.nfevals
ans =  1975

This results in the stiff system integration taking about 7 times that of the less stiff system. If the eval_rhs takes a long time to execute by itself this can easily cause longer integration times.

Time Varying Inputs

In the above example, a constant input for the torque was used. This is sometimes desired but in general is quite limiting. What if you want the input to be a function of time, the state, or the parameters (which are all valid choices)?

\begin{equation*} \mathbf{r} = \mathbf{w}(t, \mathbf{x}, \mathbf{p}) \end{equation*}

Similarly to the function that evaluates the differential equations, create an Octave/Matlab function that returns the input vector given the current time, state, and constant parameter values. Save this as eval_input.m.

function r = eval_input(t, x, p)
% EVAL_INPUT - Returns the input vector at any given time.
%
% Syntax: r = eval_input(t, x, p)
%
% Inputs:
%   t - A scalar value of time, size 1x1.
%   x - State vector at time t, size mx1 were m is the number of states.
%   p - Constant parameter structure with p items where p is the number of
%       parameters.
% Outputs:
%   r - Input vector at time t, size ox1 where o is the number of inputs.

% NOTE : x is not needed in this case, so no need to unpack it

% unpack the structure
m = p.m;
l = p.l;
g = p.g;

% define a sinusoidal time varying input
tau = m*g*l/2*sin(pi/20*t);

% pack the inputs into an 1x1 vector
r = [tau];

end

For this function to be useful a slight adjustment to eval_rhs.m needs to be made so that it accepts the input function instead of the values directly. Save this as eval_rhs_with_input.m.

function xdot = eval_rhs_with_input(t, x, w, p)
% EVAL_RHS_WITH_INPUT - Returns the time derivative of the states, i.e.
% evaluates the right hand side of the explicit ordinary differential
% equations.
%
% Syntax: xdot = eval_rhs_with_input(t, x, w, p)
%
% Inputs:
%   t - Scalar value of time, size 1x1.
%   x - State vector at time t, size mx1 where m is the number of states.
%   w - Anonymous function, w(t, x, p), that returns the input vector at
%       time t, size ox1 were o is the number of inputs.
%   p - Constant parameter structure with p items where p is the number of
%       parameters.
% Outputs:
%   xdot - Time derivative of the states at time t, size mx1.

% unpack the states into useful variable names
theta = x(1);
omega = x(2);

% evaluate the input function
r = w(t, x, p);

% unpack the inputs into useful variable names
tau = r(1);

% unpack the parameters into useful variable names
m = p.m;
l = p.l;
g = p.g;

% calculate the state derivatives
thetadot = omega;
omegadot = (-m*g*l*sin(theta) + tau)/(m*l^2);

% pack the state derivatives into an mx1 vector
xdot = [thetadot; omegadot];

end

Now you can pass in the input function as an anonymous function in similar fashion as shown earlier for eval_rhs(). Save as integrate_with_input_function.m.

% this setup of initial conditions, time values, and parameters is the same
% as before
x0 = [5*pi/180;  % [rad]
      0];        % [rad/s]
ts = linspace(0, 10, 500);  % [s]
p.m = 1;  % kg
p.l = 1;  % m
p.g = 9.81;  % m/s^2

% check if the input function works for the initial condition
display('Checking if eval_input returns an expected result:')
eval_input(5.0, x0, p)

% check if the eval_rhs function works with the input function passed in as
% an anoymous function
display('Checking if eval_rhs_with_input returns an expected result:')
eval_rhs_with_input(5.0, x0, @eval_input, p)

% create an anonymous function with the required inputs for ode45(), i.e. (t, x).
f = @(t, x) eval_rhs_with_input(t, x, @eval_input, p);

% Now the equations can be integrated.
[ts, xs] = ode45(f, ts, x0);

This design sets you up to easily swap out input functions. You can create an input function for each desired input type. For example, here is a step function, eval_step_input.m.

function r = eval_step_input(t, x, p)
% EVAL_STEP_INPUT - Returns the input vector at any given time.
%
% Syntax: r = eval_step_input(t, x, p)
%
% Inputs:
%   t - Scalar value of time, size 1x1.
%   x - State vector at time t, size mx1 where m is the number of states.
%   p - Constant parameter structure with p items where p is the number of
%       parameters.
% Outputs:
%   r - Input vector at time t, size ox1 where o is the number of inputs.

% NOTE : x and p are not used in this case

% calculate the step input
if t > 1.0
    tau = 0.5;  % [N]
else
    tau = 0.0; % [N]
end

% pack the input values into a ox1 vector
r = [tau];

end

Now integrating with the new input only requires changing the name of the anonymous function in the main script, named here as integrate_with_step_function.m.

% intial condition, time values, and parameter setup is the same
x0 = [5*pi/180;  % [rad]
      0];        % [rad/s]
ts = linspace(0, 10, 500);  % [s]
p.m = 1;  % kg
p.l = 1;  % m
p.g = 9.81;  % m/s^2

% check if the input function works for the initial condition
display('Checking step input function:');
eval_step_input(5.0, x0, p)

% check if the eval_rhs funciton works with the input function passed in as
% an anoymous function
display('Checking rhs input function:');
eval_rhs_with_input(5.0, x0, @eval_step_input, p)

% changing the input only requires changing the function name here
f = @(t, x) eval_rhs_with_input(t, x, @eval_step_input, p);

% Now the equations can be integrated.
[ts, xs] = ode45(f, ts, x0);

Outputs Other Than The States

The first type of outputs you may be interested in are functions of the states, time, inputs, and constant parameters. It is useful to create a function that can calculate these. It is typically best to do this after integration for both an organizational standpoint and computational efficiency purposes (e.g. you an leverage vectorization and broadcasting, as shown below).

\begin{equation*} \mathbf{y} = \mathbf{g}(t, \mathbf{x}, \mathbf{r}, \mathbf{p}) \end{equation*}

Example outputs for the pendulum might be the Cartesian coordinates of the pendulum bob and the energy, kinetic and potential. The equations below describe these computations:

\begin{align*} x_p = l \sin(\theta) \\ y_p = l - l \cos(\theta) \\ E_k = ml^2\omega/2 \\ E_p = mgy_p \end{align*}

Create a new function file, eval_output.m, that encodes these mathematical operations.

function y = eval_output(t, x, r, p)
% EVAL_OUTPUT - Returns the output vector at the specified time.
%
% Syntax: y = eval_output(t, x, r, p)
%
% Inputs:
%   t - Scalar value of time, size 1x1.
%   x - State vector at time t, size mx1 where m is the number of states.
%   r - Input vector at time t, size ox1 where o is the number of
%       inputs.
%   p - Constant parameter structure with p items where p is the number of
%       parameters.
% Outputs:
%   y - Output vector at time t, size qx1 where q is the number of outputs.

% unpacke the states
theta = x(1);
omega = x(2);

% unpack the parameters
m = p.m;
l = p.l;
g = p.g;

% calculate the Cartesian position of the pendulum bob
x_pos = l*sin(theta);
y_pos = l - l*cos(theta);  % position relative to stable equilibrium

% calculate the kinetic and potential energies
kinetic_energy = m*l^2*omega^2/2;
potential_energy = m*g*y_pos;

% pack the outputs into a qx1 vector
y = [x_pos; y_pos; kinetic_energy; potential_energy];

end

Now this function can be used after integrating the ODEs to compute any desired outputs. The following file, integrate_with_output.m, shows how this is done.

% same integration code as above
x0 = [5*pi/180; 0];
ts = linspace(0, 10, 500);
p.m = 1;
p.l = 1;
p.g = 9.81;
f_anon = @(t, x) eval_rhs_with_input(t, x, @eval_input, p);
[ts, xs] = ode45(f_anon, ts, x0);

% to compute the outputs at each time, you must iterate through time
ys = zeros(length(ts), 4);  % create a matrix to store the values, nxq
for i=1:length(ts)
    % the r input isn't used, so nan can be set as a placeholder
    ys(i, :) = eval_output(ts(i), xs(i, :), nan, p);
end

Vectorizing functions

It is also worth noting that Octave/Matlab code can generally be written to avoid loops, like in the above example. Slight adjustments to the output function will allow batch calculations of the outputs, as shown below in eval_output_vectorized.m:

function ys = eval_output_vectorized(ts, xs, rs, p)
% EVAL_OUTPUT_VECTORIZED - Returns the output vector at the specified times.
%
% Syntax: ys = eval_output_vectorized(ts, xs, rs, p)
%
% Inputs:
%   ts - Scalar values of time, size nx1.
%   xs - State vector at each time, size nxm where m is the number of states
%        or mx1 if n=1.
%   rs - Input vector at each time, size nxo where o is the number of inputs
%        or ox1 if n=1.
%   p - Constant parameter structure with p items where p is the number of
%       parameters.
% Outputs:
%   ys - Output vector at each time, size nxq where q is the number of
%        outputs or qx1 if n=1.

% for n=1, tranpose
if length(ts) == 1
    xs = xs';  % 1xm
    rs = rs';  % 1xo
end

% unpack the state trajectories
theta = xs(:, 1);  % size n
omega = xs(:, 2);  % size n

% unpack the parameters
m = p.m;
l = p.l;
g = p.g;

% NOTE : rs is not used in this case

% make sure to use elementwise operators, e.g. .* instead of * for
% vectorized calculations
x_pos = l.*sin(theta);
y_pos = l - l.*cos(theta);

kinetic_energy = m.*l.^2.*omega.^2./2;
potential_energy = m.*g.*y_pos;

% pack the results so the function returns a size nxq result (use commas)
ys = [x_pos, y_pos, kinetic_energy, potential_energy];

% ensures this function returns a column vector if n=1
if length(ts) == 1
    ys = ys';
end

end

Now, instead of the for loop, you can type:

ys = eval_output_vectorized(ts, xs, nan, p);

These batch, or "vectorized", calculations can be significantly faster than the loops, if that is desirable.

Outputs Involving State Derivatives

Additional outputs you may desire can also depend on the value of the time derivative of the states, i.e. \(\mathbf{\dot{x}}\), and the output function then takes this form:

\begin{equation*} \mathbf{z} = \mathbf{h}(t, \dot{\mathbf{x}}, \mathbf{x}, \mathbf{r}, \mathbf{p}) \end{equation*}

For example, the following function, eval_output_with_state_derivatives.m, calculates the radial and tangential acceleration of the pendulum bob. The tangential acceleration depends on \(\dot{omega}\).

function z = eval_output_with_state_derivatives(t, xdot, x, r, p)
% EVAL_OUTPUT_WITH_STATE_DERIVATIVES - Returns the output vector at the
% specified time.
%
% Syntax: z = eval_output_with_state_derivatives(t, xdot, x, r, p)
%
% Inputs:
%   t - Scalar value of time, size 1x1.
%   xdot - State derivative vector as time t, size mx1 where m is the number
%          of states.
%   x - State vector at time t, size mx1 where m is the number of states.
%   r - Input vector at time t, size ox1 were o is the number of inputs.
%   p - Constant parameter structure with p items where p is the number of
%       parameters.
% Outputs:
%   z - Output vector at time t, size qx1.

% unpack all of the vectors
thetadot = xdot(1);
omegadot = xdot(2);

theta = x(1);
omega = x(2);

m = p.m;
l = p.l;
g = p.g;

% calculate the radial and tangential accelerations
radial_acc = omega^2 * l;
tangential_acc = omegadot * l;

% pack the result into a qx1 vector
z = [radial_acc; tangential_acc];

end

The state derivatives are calculated internally when ode45() is called and are not stored during integration. These can be recalculated after integration for use in you primary script, e.g. as in integrate_with_derivative_output.m.

% integrate the ODEs as defined above
x0 = [5*pi/180; 0];
ts = linspace(0, 10, 500);
p.m = 1;
p.l = 1;
p.g = 9.81;
f_anon = @(t, x) eval_rhs_with_input(t, x, @eval_input, p);
[ts, xs] = ode45(f_anon, ts, x0);

% to compute the outputs at each time, loop through time evaluating f(t, x,
% r, p) and then h(t, x', x, r, p)
zs = zeros(length(ts), 2);  % place to store outputs
for i=1:length(ts)
    % calculate x' at the given time
    xdot = eval_rhs_with_input(ts(i), xs(i, :), @eval_input, p);
    % calculate the outputs that depend on x' and store them
    zs(i, :) = eval_output_with_state_derivatives(ts(i), xdot, xs(i, :), @eval_input, p);
end