This is a simple header-only C++ library for solving ordinary differential equations using the fixed-step 4th order Runge-Kutta method. It offers a rudimentary RK4 solver that is light-weight, fast and hopefully easy to use. The objective is to solve a system of ordinary differential equations (ODEs) given as
Some event functionality is included for hybrid dynamics, but it should not be trusted as this is a fixed-step size method. Kahan summation algorithm is used in the integration to compensate for the accumulation of floating point errors.
This library and its only dependency matrix_op was written with embedded systems in mind, and it has been tested on hard real-time embedded systems. It does not use RTTI or exceptions, and you can disable all dynamic heap allocations using the DO_NOT_USE_HEAP
compiler flag, but beware of stack overflows. Use the USE_SINGLE_PRECISION
compiler flag to compile for single precision systems.
Include all or some of the headers in include/
into your project.
Alternatively, you can use FetchContent() in your CMakeLists.txt
:
FetchContent_Declare(rk4_solver URL https://github.com/cinaral/rk4_solver/releases/download/<RELEASE_TAG>/src.zip)
FetchContent_MakeAvailable(rk4_solver)
Use CTest to test the library before using. Currently, four tests are available:
- Integrating a sine function,
- Solving a first-order system,
- Solving for the time response of a motor that is driven by a sinusoidal input,
- Solving for the motion of a bouncing ball (hybrid dynamics).
Some of the tests require reference solutions which is included in this repository, see Testing for details. For the minimum examples, see examples/
.
You must first create an integrator object:
rk4_solver::Integrator<x_dim, Dynamics> integrator(dynamics, &Dynamics::ode_fun, time_step, t_init = 0);
ODE function should be a member of Dynamics
, which is a class you will provide for the dynamics of your system. The ODE function must be of type OdeFun_T
:
//* dx/dt = ode_fun(t, x)
void Dynamics::ode_fun(t, x, OUT: dt_x);
If you need an event function, you can use an event function of type EventFun_T
. It is fine to use the same object for both the integrator and the event function:
//* x_plus = event_fun(t, x)
void Events::check(t, x, OUT: x_plus);
WARNING: This is a fixed-step size method and therefore the event detection will be approximate within the time step size.
OdeFun_T
and EventFun_T
are function pointers defined in types.hpp:
using OdeFun_T = void (T::*)(const Real_T t, const Real_T (&x)[X_DIM], Real_T (&dt_x)[X_DIM]);
using EventFun_T = bool (T::*)(const Real_T t, const Real_T (&x)[X_DIM], Real_T (&x_plus)[X_DIM]);
WARNING: By default, Real_T
is double
. Use USE_SINGLE_PRECISION
compiler flag to set to float
.
Call step(...)
to forward the integration by one step:
void
step(
const Real_T t,
const Real_T (&x)[X_DIM],
Real_T (&x_next)[X_DIM]
);
Call loop(...)
to loop through the integration steps as follows if you want to discard all intermediate values except the final state:
void
loop(
Integrator integrator,
OPTIONAL: Event event,
const Real_T t_init,
const Real_T (&x0)[X_DIM],
Real_T &t,
Real_T (&x)[X_DIM]
);
Otherwise, you can use the overloaded version of loop(...)
that takes an array for time and another array for state history to save all intermediate values of the solution:
void
loop(
Integrator integrator,
OPTIONAL: Event event,
const Real_T t_init,
const Real_T (&x0)[X_DIM],
Real_T (&t)[T_DIM],
Real_T (&x)[T_DIM][X_DIM]
);
This is a minimum example of a single integration step. The dynamics of the system is defined in the Dynamics
class. The ODE function is defined as a member of Dynamics
and is passed to the integrator as a function pointer. The integrator is then used to integrate the system for a single step.
#include "rk4_solver/integrator.hpp"
//...
Dynamics::ode_fun(const Real_T, const Real_T (&x)[x_dim], Real_T (&dt_x)[x_dim])
{
dt_x[0] = x[1];
dt_x[1] = u;
}
//...
Dynamics dynamics;
int main()
{
rk4_solver::Integrator<x_dim, Dynamics> integrator(dynamics, &Dynamics::ode_fun, time_step);
Real_T t;
Real_T x[x_dim];
integrator.step(t_init, x_init, t, x); //* integration step
//...
}
See step-example.cpp for details.
Here is an example of an integration loop where we save all intermediate values.
#include "rk4_solver/cum_loop.hpp"
//...
int main()
{
rk4_solver::Integrator<x_dim, Dynamics> integrator(dynamics, &Dynamics::ode_fun, time_step);
Real_T t_arr[t_dim];
Real_T x_arr[t_dim][x_dim];
//* save the intermediate values
rk4_solver::loop(integrator, t_init, x_init, t_arr, x_arr);
//...
}
See loop-example.cpp for details.
In this example we discard the intermediate values, and stop at the first ground contact. Ground contact is checked using the current vertical displacement, and the displacement is reset to zero at impact.
#include "rk4_solver.hpp"
//...
struct Dynamics {
//...
bool
event_fun(const Real_T, const Real_T (&x)[x_dim], Real_T (&x_plus)[x_dim])
{
bool did_occur = false;
if (x[0] <= 0) {
x_plus[0] = 0;
x_plus[1] = -e_restitution * x[1];
did_occur = true;
}
return did_occur;
}
};
Dynamics dynamics;
int
main()
{
rk4_solver::Integrator<x_dim, Dynamics> integrator(dynamics, &Dynamics::ode_fun, time_step);
rk4_solver::Event<x_dim, Dynamics> event(dynamics, &Dynamics::event_fun);
Real_T t;
Real_T x[x_dim];
//* discard the intermediate values and stop at first event
rk4_solver::loop<t_dim>(integrator, event, t_init, x_init, t, x, true);
//...
}
See event-example.cpp for details.
There are two benchmark tests:
- A step integration loop without final time, and intermediate values are discarded.
- A cumulative integration loop with final time, and intermediate values are saved.
The benchmark test is a 3rd order linear system compiled using g++ with -O3
optimization level. A desktop Intel i7-9700K at 3.60 GHz processor with 32 GB of memory was used to obtain the following sample benchmarks:
Flags | Loop (million steps per second) | Cumulative Loop (million steps per second) |
---|---|---|
None (Default) | 18.5 | 28.6 |
USE_SINGLE_PRECISION |
18.5 | 30.5 |
DO_NOT_USE_HEAP |
18.6 | 28.2 |
DO_NOT_USE_HEAP and USE_SINGLE_PRECISION |
18.6 | 30.7 |
Using the USE_SINGLE_PRECISION
or the DO_NOT_USE_HEAP
flag does not affect the performance meaningfully. Use them if your target platform can benefit from them, remember to test and benchmark.
WARNING: Your stack can easily overflow for large problems with the DO_NOT_USE_HEAP
flag and should be avoided if it is not necessary.
Reference solutions are required for some tests, which are included in test/dat/
. matrix_rw library is used to read and write from file, the *.dat
files are comma and newline delimited. If you have access to MATLAB, the formatting is compatible with writematrix
and readmatrix
.
You may need to generate new reference solutions to update the existing tests or add new tests. run_test.m can be used to (re)generate the reference solutions if you have MATLAB. The MATLAB tests in test/matlab/
are optional, but they can be used to visually verify the results.