This is just a small update on my experiments with the Arduino. I implemented the runge-kutta-method for solving a multibody system a few weeks ago. So this is a working implementation of the standard 4th-order runge-kutta ODE (ordinary differential equations) solver for the arduino platform, something I haven’t seen elsewhere.

(Note: If you are just looking for the arduino sketch, you’ll find a link at the end of this blog post. If you are wondering about related questions or have trouble adopting it to your needs, please leave a comment. Also see this post on how numerical integration of differential equations works.

Update in August 2016: See also my new post on achievable simulation rates with an Arduino Uno/Nano and Due)

My main goal was to get a better grip on simulation speeds. *Why is my simulation so slow?* is a question I’m really thinking about a lot. Since I wasn’t able to find a good example of an arduino sketch for ODE solving online, this may be interesting anyway.

The sketch can be used to simulate systems of ordinary differential equations (ODEs) in state form:

I currently include a simple 2-mass multibody system (→ 4 state variables) but this can be changed easily. My Arduino Uno achieves something around 1000 simulation steps per second which is less than I expected. This seems to result mainly from the many floating point operations the sketch requires (each +, – or * with floating point numbers uses approx. 100 CPU cycles).

For comparison: The same code compiled in C++ on my laptop achieves around 2 million simulation steps per second. Also C++ allows me to use double precision for the variables, resulting in smaller rounding arrors (arguably a significant difference for some simulations).

Nonetheless, the arduino sketch is still fast enough if you want to do some simple physics modelling. For example something like:

- take analog input from a human,
- have the arduino simulate a physical system and
- output the behavior to either a screen or an actuator.

(Update: I now did exactly that – see the update here)

Here’s the link to the arduino sketch (updated version from August 2016).

*Update on April 27th 2015: Fixed some minor grammar issues and also changed the caption.*

Update on Feb. 26th 2016: The link to the sketch at least didn’t work for me. –> now updated.

Update on May. 7th 2016: Added some clarifications in the introduction, making it easier to find the code. I’m currently looking at ODE solvers in C++, mainly the odeint and GSL packages, so maybe there will be an update some time in the future.

Update on Aug. 24th 2016: Updated the arduino sketch and removed the bonus link (blog no longer active).

Many thanks!

How will be the void ODE(double t, double Y[], double Y_p[] for a Lotka-Volterra equations seen, please?

dx/dt = a*x – b*x*y

dy/dt = d*x*y – g*y

Andrew

Hey Andrew,

try this:

col 33: change the number of state variables to two:

const int n=2;

col 37: modify initial vector to a suitable one (with only two values instead of four, add constants for a, b, c, g instead of c1, c2, … below.

And the ODE function – try this:

void ODE(double t, double Y[], double Y_p[]) {

// x = Y[0], y = Y[1] in your equation

Y_p[0]= a*Y[0] – b*Y[0]*Y[1];

Y_p[1]=d*Y[0]*Y[1] – g*Y[1];

}

You might as well want to add intermediate variables just for clarity:

void ODE(double t, double Y[], double Y_p[]) {

double x = Y[0];

double y = Y[1]; // C++ is case sensitive, so this should work.

Y_p[0]= a*x – b*x*y; // Y_p[0] = dx / dt

Y_p[1]=d*x*y – g*y; // Y_p[1] = dy / dt

}

Please note – I haven’t tested this and please compare your results to some reference. If this doesn’t help, please write a short note.

Cheers, Daniel

Many thanks for your post!

I wonder if it is possible to solve Riccati equation (second order ode) in this way. Have you ever tried doing that? Do you have any suggestion?

Dave

Hi Dave,

sorry for the late reply. Riccati ODEs appear to be a family of ODEs, not a single one – are you sure you are looking for a numerical solver instead of a general analytic approach? (just taking guesses here). If you do, this is possible:

Transform the second-order ODE into a first-order one with two state variables. The same problem arises with dynamics where newtons law leads to second order ODEs with F/m = acceleration. This is solved by using both position and velocity as state variables. If x is your state vector, you use

x = (pos, vel)

dx/dt = (vel, accel) = (x[2], F/m)

(with F usually being a function of your current state).

I use exactly this approach in the example provided.

There’s some general info on how numerical integration works in my post here, but maybe you have someone available to discuss options with (this would be so much easier if we were both in the same room!).

Cheers, Daniel

Hey,

Amazing idea.

I was wondering if it can be used for prediction. Per example, a set of temperatures is given and I want to predict in a minute or hour. Would this library and if so then how??

Thank you for reading.