Numerical integration

Differential equations are amazing and useful stuff. Not a day passes by when you don’t solve at least one. And if you do not - don’t you feel like you have wasted a day of your life?

Sadly, as much as you’d wish to, you don’t always have your trusty Matlab (or Mathematica etc.) with you. Never fear, a sheet of paper (and later - a spreadsheet) will be more than enough!

The Euler method

You already know how to solve differential equations in memory:

A car is travelling with the velocity of 10 m/s and is now at the 100 m mark on the race track. What will its position be after 1.5 seconds?
The answer is of course: 115 m.

What took the derivative of the position of the car (the velocity) and added it to its current position after multiplying it by the time step (1.5 s). This is the essence of the Euler’s method of numerical integration:

Basically, you add to the current state (the “position”) the derivative of the state (the “velocity”) multiplied by the length of the time step. If the steps are small, the result should be good enough.

Let’s solve a simple differential equation numerically (you might recognize this as an inertial element):

We do it in the following steps:

1. Prepare the table. It should have a column for time (t), the input (u), the state (x), the state derivative column (dx) and the output column. Fill in the time values at desired time steps.

t	u	x	dx	y
0.0	…	…	…	…
0.1	…	…	…	…
0.2	…	…	…	…

2. Next, fill in the input values at different values of time. Here, we can use an unit step for example (so all the input values will be 1). We can also pick the initial value for the state. Let’s say it’s 0.

t	u	x	dx	y
0.0	1	0	…	…
0.1	1	…	…	…
0.2	1	…	…	…

3. Let’s start. First, we calculate the derivative dx at t=0.0. We use the equation inserting the values from the current row: . We can calculate the output at the same time: .

t	u	x	dx	y
0.0	1	0	1	0
0.1	1	…	…	…
0.2	1	…	…	…

4. For the second row, we have to update the state. We add the dx from the row above, multiplied by the step : .

t	u	x	dx	y
0.0	1	0	1	0
0.1	1	0.1	…	…
0.2	1	…	…	…

5. We repeat the step 3 for the second row now, calculating new dx and the new y.

t	u	x	dx	y
0.0	1	0	1	0
0.1	1	0.1	0.9	0.1
0.2	1	…	…	…

6. Keep repeating steps 3 and 4 for the subsequent rows:

t	u	x	dx	y
0.0	1	0	1	0
0.1	1	0.1	0.9	0.1
0.2	1	0.19	0.81	0.19
0.3	1	0.27	0.73	0.27
0.4	1	0.34	0.66	0.34
0.5	1	0.41	0.59	0.41
0.6	1	0.47	0.53	0.47
0.7	…	…	…	…

(If you make a plot of this it should look like the one in fig. 4.)

Solving various differential equations make for an excellent pastime for the boring and useless meetings and seminars. You will rarely get scorned for doing maths on such occasions; certainly not as if you were doing sudokus, crosswords or mind-absently doodling.

Solving differential equations numerically by hand will command great respect among your peers. Whenever I explained what I was doing to a curious inquiring soul, they never bothered me with questions again.
They evidently understand how important math is!

The Runge-Kutta method

The Euler method is fast and easy to use - and certainly sufficient for simple demonstrations. It is not very precise though.

Another popular method that helps with that problem is the classic Runge-Kutta fourth-order method. In RK4, the next state value is calculated as:

We simply calculate the derivatives four times. First - based on the current time and state values. The next approximation is based on a what-if scenario: what if we used thus calculated state value to compute derivatives at half the time step from now ()?

The third coefficient is calculated the same way, but taking the newly made prediction into account. At last, the what-if derivative is calculated using k3 prediction at the end of the time step ().

A simple implementation of RK4 in JS might look like this:

function solve(fun, t, u, x, dt) {
      // calculate derivatives k1
      var k1 = fun(t, u, x);

      // update state for calculating x2
      var x1 = x.slice();
      for (var i = 0; i < x.length; ++i) {
        x1[i] += k1[i] * dt/2;
      }

      // calculate derivatives k2
      var k2 = fun(t+dt/2, u, x1);

      // update state for calculating x3
      var x2 = x.slice();
      for (var i = 0; i < x.length; ++i) {
        x2[i] += k2[i] * dt/2;
      }

      // calculate derivatives k2
      var k3 = fun(t+dt/2, u, x2);

      // update state for calculating x4
      var x3 = x.slice();
      for (var i = 0; i < x.length; ++i) {
        x3[i] += k3[i] * dt;
      }

      // calculate derivatives k4
      var k4 = fun(t+dt, u, x3);

      // update state with the final value
      var nx = x.slice();
      for (var i = 0; i < x.length; ++i) {
        nx[i] += dt/6 * (k1[i] + 2*k2[i] + 2*k3[i] + k4[i]);
      }

      return nx;
};

You could use it like this:

// a function that calculates the derivatives
function odefun(t, u, x) {
    var dx = [ 0.0, 0.0 ];
    dx[0] = x[1];
    dx[1] = u - 0.1*x[0] - 0.1*x[1]; 

    return dx; 
};

var x = [ 0.0, 0.0 ]; // initial conditions
for (var t = 0.0; t < 10.0; t += 0.1) {
    u = 1.0; // the input
    x = solve(odefun, t, u, x, 0.1);
    console.log(t, x);
};

RK4 spreadsheet

While Euler’s method is perfectly portable and easy on the hand while drawing it out on the paper, RK4 is somewhat more complicated. I made a little spreadsheet that can be used to solve simple LTI systems in LibreOffice Calc or Excel.

The spreadsheet (and other materials mentioned here) is available HERE

LTI systems (linear time-invariant systems) are a class of systems taking the following form:

A, B, C, D are matrices here, which describe the linear relationships between the input signals u, states x and output signals y (x, u and y are vectors). The above equations are presented in condensed notations, but are really quite simple. To be more verbose:

The A matrix (state matrix) describes how the states influence the state derivatives. The B matrix (input matrix) relates the influence of input signals on the states. The output signal y is generated based on C matrix (iutput matrix) multiplied by the states and the D matrix (feedthrough matrix) directly passing through the input signals to the output.

The spreadsheet layout is presented in fig. 1. The top part contains the data you can specify as the user (marked in yellow). This data consists of the A, B, C and D matrices elements, the initial conditions (upper left corner), time step (bottom right corner) and the input signal generator.

Fig. 1. Solving 5x5 MIMO LTI state-space system in a spreadsheet.

The dark yellow cells take values of 0 or 1 and are used to enable and disable plotting of the corresponding input and output channels.

You can solve systems of up to 5 states and use 5 input and 5 output channels as well. The inputs are generated as a sum of 5 sinusoidal functions, which you can parametrize in the panel at the botom of the user input area. The parameters for the signle channel are: the amplitude A, the frequency ω [rad/s], the phase φ [rad] and the starting time t0 [s]:

The input panel makes it possible to generate a range of popular signals:

Unit step

Unit step is easily made by setting the amplitude A to 1, the frequency ω to 0, and the phase φ to π/2. This effectively makes the generated input to be a cosine calculated at 0 - or simply constant 1.
You can also use the fourth row (t0) to specify at which time the signal should be enabled.

Sine wave

Simply insert the desired values of the amplitude, frequency and phase.

Square wave

We can use Fourier expansion to approximate a square wave of period L as:

An example using 5 input channels to generate a square wave approximation (see fig. 2) is located at examples/square_wave.ods.

Fig. 2. Square wave approximated with Fourier expansion.

Examples

Below are a few examples of the spreadsheet usage.

Integrator

The simplest example is an integrator element, described with the following equation:

Its state-space realization is:

where u is the input, x is the integrator state and y is the output. Ti is the integrator time constant - the time at which the integrator output reaches the level of 1 when subjected to an unit step input.

The state-space matrices, based on the equation above are:

The spreadsheet for this example is available in the examples/integration.ods file. The output of the integrator with the unit step input is shown in fig. 3.

Fig. 3. Integrating an unit step (Ti = 50 s).

Try changing the input to a sine function. The output should then take form of the sine integrated - or a cosine. Hopefully the numerical errors won’t accumulate too much!

Inertia

Another classical example is an inertial element. It represents an object that takes a while to reach its set-point value (e.g. a temperature in a heater, rotations of the engine etc.).

It is described with the following equations:

You only need one state to deal with the inertial element. The matrices look like this:

Fig. 4 shows the inertial element response. The example is located in examples/inertia.ods.

Fig. 4. Inertial element response to an unit step (k = 1, T = 10 s).

Mass-spring-damper system

Another classic. A harmonic damped oscillator is one of the cornerstones of physics. It is defined by the following equations:

We have to use two state variables (x1 representing the position and x2 representing the velocity) to reduce the order of that differential equation and arrive at the following representation:

The state representation matrices are:

The system response in shown in fig. 5. The example is to be found in examples/mass.ods file.

Fig. 5. Mass-spring-damper system response to an unit step (parameters: m = 1, k = 1, b = 2.5). The blue line is the mass position and the red one is the velocity.

PID

The ideal PID controller response is described by the following differential equation:

It is not unfortunately possible to implement an ideal derivative, and thus the following approximation is used (here in the form of a transfer function):

Such a PID controller could be realized in the state-space in the following form:

so that the state-space matrices are:

The x1 variable handles the integrator state, while the x2 contains the derivative’s inertial element state. The proportional gain is implemented with the use of the feedthrough matrix D. A sample PID response simulation is presented in fig. 6. The PID controller example can be found in the examples/pid.ods spreadsheet.
Too bad it is somewhat tricky to use spreadsheets for online control!

Fig. 6. Simulating PID controller response (Kp = 1, Ti = 100, Td = 1, T = 1).

It should be easy enough to use the spreadsheet to simulate the behaviour of other models. One idea could be to use it for the simulation of other models I have presented here yet: the cart system, the tanks or the oscillating pendulum (some of these models would have to be linearized first though).

​