Solving ordinary differential equations in C++

• Download source code - 13 KB

Introduction

This article introduces the C++ framework odeint for solving ordinary differential equations (ODEs), which is based on template meta-programming. I think this framework has some nice advantages over existing code on ODEs, and it uses templates in a very elegant way. Furthermore, odeint has a STL-like syntax, and is independent of a specific container. The whole project lives in the boost sandbox; the code presented here is more or less a snapshot of the current development.

Background

I will only very briefly describe ordinary differential equations. If there is some interest in a more detailed explanation of ODEs, I can extend this part in future versions of the article.

Ordinary differential equations

An ordinary differential equation describes the evolution of some quantity x, and has the form:

d x(t) / d t = f( x(t) )

The function f defines the ODE, and x and f can be vectors. Associated with every ODE is an initial value problem (IVP) that is the ODE, and an initial value x(t0)=x0. The solution of the initial value problem is the temporal evolution of x(t), with the additional condition that x(t0)=x0, and it can be shown that every IVP has a unique solution. Of course, ordinary differential equations are not restricted to temporal problems, hence the variable t can be replaced by another quantity, like a spatial coordinate. ODEs and their relative PDEs (partial differential equation) are very important in nearly all scientific disciplines. All major equations in physics fall in this class, like Newton's law for classical physics, the Maxwell's equations for electromagnetism, the Schrödinger equation and its relativistic generalizations for the quantum world, or Einstein's equation for general relativity. Very often, the spatial axes of a PDE are discretized, and an ODE on a lattice results.

In the following figure, an example of an ODE from chaos theory is shown: the famous Lorenz attractor. It was the first example of a deterministic chaotic system, and triggered a huge amount of scientific work. The Lorenz attractor is described by a set of coupled ordinary differential equations:

d x1 / dt = sigma * ( x2 - x1 )
d x2 / dt = R * x1 - x2 - x1 * x3
d x3 / dt = x1 * x2 - b * x3

with:

sigma = 10
R = 28
b = 8 / 3

This equation has no analytical solution, such that it can only be solved numerically.

Numerical solutions of ODEs

To solve ODEs numerically, various methods exist; all of them discretize the time. The easiest method is surely the explicit Euler scheme, which writes the derivative as the difference quotient:

d x(t) / d t = x(t+dt) - x(t) / dt

Then, some basic algebraic manipulations yield:

x(t+dt) = x(t) + dt*f(x(t))

such that x(t+dt) can be obtained from the knowledge of x at the time t. Note again, that x does not have to be a scalar, but can be a vector. The rest is very simple. We choose an initial condition x(t=0)=x0, and iteratively apply the Euler step to x(t) to obtain the evolution of x(t). Of course, more advanced solvers exist, and the most commonly used solver is probably the Runge-Kutta method of fourth order.

In the following source code, the function f(x(t)) defining the ODE is synonymously called system or dynamical system.

The main ideas

Now, we come to the design aspects of the library. We clearly need two main building blocks:

1. A stepper function, which calculates the state x(t+d t) from the knowledge of x(t).
2. An integrate function, which performs the iteration, hence which applies the stepper function many times.

Stepper methods

We try to avoid abstract classes, since this costs some extra performance. Instead, we use concepts which describe how a type or class has to behave, and whose methods have to be available. For the stepper function, we need the following ingredients:

template< class Container , class Time = double , 
          class Resizer = resizer< Container > >
class ode_step
{
 // provide basic typedefs
 public:

    typedef Container container_type;
    typedef Resizer resizer_type;
    typedef Time time_type;
    typedef const unsigned short order_type;
    typedef typename container_type::value_type value_type;
    typedef typename container_type::iterator iterator;

 // public interface
 public:

    order_type order() const;

    template< class DynamicalSystem >
    void next_step( DynamicalSystem &system ,
                    container_type &x ,
                    const container_type &dxdt ,
                    time_type t ,
                    time_type dt );

    template< class DynamicalSystem >
    void next_step( DynamicalSystem &system ,
                    container_type &x ,
                    time_type t ,
                    time_type dt );
};

The container type defines how the stated is stored. It has to fulfill some basic requirements, like the begin() and end() iterator and the value type. The resizer is an abstraction of the resize() functionality of the sequence concept. Its usage will be shown below. The next_step functions calculate x(t+dt) from the knowledge of x(t), and they transform x(t) in-place, meaning that x(t) is replaces by x(t+dt) within the next_step function. The argument system defines the ODE. This can be a simple function with three arguments: system( state , dxdt , t ), or a class with the appropriate brackets operator()(state,dxdt,t). Two versions of next_step exist: one which calculates x(t+dt) with the explicit knowledge of f(x(t)), and one which calculates this value internally. In principle, the second version only calls the first version. Both functions exist, because in many situations, the user needs the knowledge of f(x(t)) and will call the system method with x(t). To avoid a double call of system, the next_step version with the explicit declaration of f(x(t)) is introduced.

Note that the above class definition is not found in the code; it is only used here, to show the stepper functionality. All specific classes possessing these typedefs and methods are said to be stepper classes; we also say that all classes with these definitions fulfill the stepper concept.

Now, we look at a specific method, the Euler solver:

template< class Container , class Time = double , 
          class Resizer = resizer< Container > >
class ode_step_euler
{
 // provide basic typedefs
 public:

    typedef Container container_type;
    typedef Resizer resizer_type;
    typedef Time time_type;
    typedef const unsigned short order_type;
    typedef typename container_type::value_type value_type;
    typedef typename container_type::iterator iterator;

 // private members
 private:

    container_type m_dxdt;
    resizer_type m_resizer;

 // public interface
 public:

    order_type order() const { return 1; }

    template< class DynamicalSystem >
    void next_step( DynamicalSystem &system ,
                    container_type &x ,
                    const container_type &dxdt ,
                    time_type t ,
                    time_type dt )
    {
        // computes x = x + dt * dxdt
        detail::it_algebra::increment( x.begin() , x.end() , dxdt.begin() , dt );
    }

    template< class DynamicalSystem >
    void next_step( DynamicalSystem &system ,
                    container_type &x ,
                    time_type t ,
                    time_type dt )
    {
        m_resizer.adjust_size( x , m_dxdt );
        system( x , m_dxdt , t );
        next_step( system , x , m_dxdt , t , dt );
    }
};

The definition is straightforward, only three notes:

• All iterator operations have been put in their own functions. This has the advantage that the code is small and clear, since no iterators need to be defined in the stepper functions. Furthermore, it could be possible to specialize and optimize iterator functions for specific iterators.
• The resizer checks that the internally stored derivative m_dxdt has the same size as the state. Only, using the resize() function of the sequence concept is not sufficient here. For example, the std::tr1::array has no resize() and could not be used with this code. But a specialized version of the resizer for arrays exists, which will avoid calling resize(). Furthermore, the resizer concepts could allow the usage of matrix types or types from Blitz++ or MTL.
• The resizer functionality is necessary if the system size changes during the evolution; think, for example, of an ODE defined on a network with a non-constant number of nodes.

This is not the end of the story of the stepper functions. odeint also provides an advanced stepper concept which calculates the numerical errors made during each step. With this error calculation, it is possible to adapt the step-size dt during each step in the iteration. But this concept will not be described here.

The following code demonstrates the usage of the stepper function:

#include <iostream>

#include <tr1/array>

#include <odeint.hpp>


#define tab "\t"

using namespace std;
using namespace odeint;

typedef std::tr1::array< double , 3 > state_type;

const double sigma = 10.0;
const double R = 28.0;    
const double b = 8.0 / 3.0;

void lorenz( state_type &x , state_type &dxdt , double t )
{                                                         
    dxdt[0] = sigma * ( x[1] - x[0] );                    
    dxdt[1] = R * x[0] - x[1] - x[0] * x[2];              
    dxdt[2] = x[0]*x[1] - b * x[2];                       
}                                                         

int main( int argc , char **argv )
{
    const double dt = 0.01;

    state_type x = {{ 1.0 , 0.0 , 0.0 }};
    ode_step_euler< state_type > stepper;

    double t = 0.0;
    for( size_t oi=0 ; oi<10000 ; ++oi,t+=dt )
    {
        stepper.next_step( lorenz , x , t , dt );
        cout << x[0] << tab << x[1] << tab << x[2] << endl;
    }

    return 0;
}

Here, lorenz could also be defined as a class:

class lorenz_class
{
 public:
    void operator()( state_type &x , state_type &dxdt , double t )
    {
        dxdt[0] = sigma * ( x[1] - x[0] );
        dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
        dxdt[2] = x[0]*x[1] - b * x[2];
    }
}

which is then called via:

lorenz_class lorenzo;
stepper.next_step( lorenzo , x , t , dt );

With classes more complicated, dynamical systems (ODEs) can be created.

Integration function

The next main idea is to provide one or more functions to automatically iterate the state of the ODE, with the possibility to observe the state in each step. This is done with the integrate function:

template< class Stepper , class DynamicalSystem , class Observer >
void integrate( Stepper &stepper ,
                DynamicalSystem &system ,
                typename Stepper::time_type start_time ,
                typename Stepper::time_type dt ,
                typename Stepper::container_type &state ,
                typename Stepper::time_type end_time ,
                Observer &observer )
{
    if( start_time > end_time )
        throw std::invalid_argument( "integrate() : start_time > end_time" );

    observer( start_time , state , system );
    while( start_time < end_time )
    {
        stepper.next_step( system , state , start_time , dt );
        start_time += dt;
        observer( start_time , state , system );
    }
}

which solves the ODE with constant steps dt in the time interval start_time, end_time. The initial state has also to be provided. Furthermore, the concept of an observer is introduced, although this is not exactly an Observer pattern. Again, the observer can have any type, it only has to have the bracket operator, operator()( time , state , system ). For example, lambda expressions from boost can be used like:

integrate( stepper , system , 0.0 , dt , x , 10.0 , 
           cout << _1 << "\t" << _2[0] << "\n" );

or:

vector<double> x;
back_insert_iterator< vector<double> > iter( x );
integrate( stepper , system , 0.0 , dt , x , 10.0 , var(*iter++) = _2[0] );

Furthermore, an integrate version with an adaptive step-size exists.

Summary and outlook

I have shown you some design aspects, and the usage of odeint - a C++ framework for solving ordinary differential equations. It is very easy to use and very flexible. The performance is also very good. I have made one test run, where I compared odeint with the gsl. As a test system, the Lorenz system was used. It turned out that odeint is about two times faster than the gsl routines.

The source code provided with this article is a snapshot of the current development. It is not complete, and will be extended in the near future. For example, now odeint contains only the Euler method and Runge-Kutta 4, but higher order methods will be introduced very soon. We have a small Todo list, which has to be completed before the final release:

• More methods: RK8, RK78, Midpoint, RK45
• An observer container
• Adaptors for Blitz++, MTL, and boost::ublas
• Symplectic integrators for separable Hamiltonian system
• Documentation

If you like to contribute to this library, or have some interesting points, criticism, or suggestions, you are cordially invited.

History

• 9.11.2009 - Initial document.

You must Sign In to use this message board.

FAQ    Noise Tolerance Very HighHighMediumLowVery Low   Layout Expand Posts & RepliesThread ViewNo JavascriptNo JS + Preview   Per page 102550

Msgs 1 to 10 of 10 (Total in Forum: 10) (Refresh) FirstPrevNext

How to compile it? caviiar 13:38 30 Nov '09   Hi. First of all, thanks for such an elegant and apparently simple method to solve ODEs numerically. Only one question: Could you add a little bit about how to compile it? I am not very proficient in C++, and I am getting many errors trying to compile it with MingW. A simple paragraph with 1. do this, 2. do that, 3... etc to compile one of your examples could help very much. In essence, many of us just need a model to imitate. Thanks! Sign In·View Thread·PermaLink Re: How to compile it? headmyshoulder 1:25 1 Dec '09   Ok, thanks for this hint. I will do this very soon. The whole library is header-only and you only have to include the odeint.hpp. Nevertheless, a Makefile is provided with the sources, doing the job of compiling the examples, at least at most unixoid os. Alternativly you can compile everything by yourself and you only have to provide the current working directory as a header search path. With g++ or gcc the command is the following: g++ -I. expample.cc For lorenz_integrate_constant_step.cpp you also need the boost library. Here, the command is g++ -I. -I$(BOOST_ROOT) expample.cc Sign In·View Thread·PermaLink Re: How to compile it? caviiar 11:04 1 Dec '09   Thanks so much for your response! I had tried the first option before. This is what I get when trying to compile from a folder in which I unzipped all files you provided: C:\odeint>g++ -I. lorenz_simple.cpp lorenz_simple.cpp:20:21: tr1/array: No such file or directory In file included from ./odeint/euler.hpp:17, from ./odeint.hpp:12, from lorenz_simple.cpp:22: ./odeint/resizer.hpp:47: error: `std::tr1' has not been declared ./odeint/resizer.hpp:47: error: `array' was not declared in this scope ./odeint/resizer.hpp:47: error: wrong number of template arguments (2, should be 1) ./odeint/resizer.hpp:18: error: provided for `template<class Container> class od eint::resizer' ./odeint/resizer.hpp:47: error: expected unqualified-id before '>' token lorenz_simple.cpp:30: error: `std::tr1' has not been declared lorenz_simple.cpp:30: error: expected initializer before '<' token lorenz_simple.cpp:36: error: variable or field `lorenz' declared void lorenz_simple.cpp:36: error: `state_type' was not declared in this scope lorenz_simple.cpp:36: error: `x' was not declared in this scope lorenz_simple.cpp:36: error: `state_type' was not declared in this scope lorenz_simple.cpp:36: error: `dxdt' was not declared in this scope lorenz_simple.cpp:36: error: expected primary-expression before "double" lorenz_simple.cpp:37: error: initializer expression list treated as compound exp ression lorenz_simple.cpp:37: error: expected `,' or `;' before '{' token lorenz_simple.cpp:46: error: `state_type' has not been declared lorenz_simple.cpp:46: error: `state_type' has not been declared lorenz_simple.cpp:47: error: ISO C++ forbids declaration of `x' with no type lorenz_simple.cpp:47: error: ISO C++ forbids declaration of `dxdt' with no type lorenz_simple.cpp: In member function `void lorenz_class::operator()(int&, int&, double)': lorenz_simple.cpp:48: error: invalid types `int[int]' for array subscript lorenz_simple.cpp:48: error: invalid types `int[int]' for array subscript lorenz_simple.cpp:48: error: invalid types `int[int]' for array subscript lorenz_simple.cpp:49: error: invalid types `int[int]' for array subscript lorenz_simple.cpp:49: error: invalid types `int[int]' for array subscript lorenz_simple.cpp:49: error: invalid types `int[int]' for array subscript lorenz_simple.cpp:49: error: invalid types `int[int]' for array subscript lorenz_simple.cpp:49: error: invalid types `int[int]' for array subscript lorenz_simple.cpp:50: error: invalid types `int[int]' for array subscript lorenz_simple.cpp:50: error: invalid types `int[int]' for array subscript lorenz_simple.cpp:50: error: invalid types `int[int]' for array subscript lorenz_simple.cpp:50: error: invalid types `int[int]' for array subscript lorenz_simple.cpp: In function `int main(int, char**)': lorenz_simple.cpp:58: error: `state_type' was not declared in this scope lorenz_simple.cpp:58: error: expected `;' before "x" lorenz_simple.cpp:59: error: `state_type' cannot appear in a constant-expression lorenz_simple.cpp:59: error: template argument 1 is invalid lorenz_simple.cpp:59: error: template argument 3 is invalid lorenz_simple.cpp:59: error: invalid type in declaration before ';' token lorenz_simple.cpp:64: error: request for member `next_step' in `stepper', which is of non-class type `int' lorenz_simple.cpp:64: error: `x' was not declared in this scope Any ideas? Thanks in advance! BTW, running Win7, MingW 5.1.6 Sign In·View Thread·PermaLink Re: How to compile it? headmyshoulder 11:57 1 Dec '09   Ok, your platform does not support tr1 headers [1], or at least does not find them. I have no ad-hoc solution for this problem right now, but will change this soon. In the meantime you can do is the following: 1. remove the resizer< array > specialization in odeint/resizer.hpp 2. change typedef array< double , n > state_type to typedef vector< double > state_type 3. allocate enough memory for the vectors. One examples shows, that arrays are faster due to their fixed size. Sign In·View Thread·PermaLink Re: How to compile it? caviiar 11:13 3 Dec '09   Thanks for that hint. Now there seems to be a smaller problem, which I am not sure how to solve: C:\odeint>g++ -I. lorenz_simple.cpp lorenz_simple.cpp: In function `int main(int, char**)': lorenz_simple.cpp:59: error: `x' must be initialized by constructor, not by `{...}' Any additional hint? I hope I am no offending you with my lack of expertise in C++. The bright side is that finding and solving these difficulties will make your library more accessible to people like me: competent mathematicians but mediocre programmers. Thanks! Sign In·View Thread·PermaLink Re: How to compile it? headmyshoulder 11:21 3 Dec '09   This is still an relict from tr1/array. If you have changed state_type to vector (what I suppose, since no other errors occur) you only have to replace the line state_type x = {{ 1.0 , 0.0 , 0.0 }}; by state_type x(3); x[0] = 1.0; x[1] = 0.0; x[2] = 0.0 This initializes a vector of length 3 with the values (1,0,0). The double braces can only be used to initialize tr1/array's. Hope this helps. Sign In·View Thread·PermaLink Re: How to compile it? caviiar 14:35 3 Dec '09   It worked!!! THANKS for your help. I will solve my problem, document a step-by-step solution in my web page, and post the link here relatively soon. Hopefully that will help others too. I had to quit on Matlab because the time it took to produce a bifurcation manifold was unacceptable. I was looking simultaneously at the classic DLSODE and a couple other options. I could produce results first with with this library. GREAT JOB! It is soooo simple to manipulate the dynamical system I am studying... Please keep improving this library and expanding it. I have plenty of research models for testing and will be glad to contribute (my focus is on analysis, but numerics is unavoidable). All best, J. Sign In·View Thread·PermaLink I suggest you remind yourself about numerical analysis. MicroImaging 16:45 12 Nov '09   When I took Numerical Analysis which included error analysis of various techniques for numerical solution of ODES, From Atkinson "An Introduction to Numerical Analysis". The Runge-Kutta Fehlberg 4 order combined with the extra evaluation of the 5th order gave an automatic error analysis and thus allowed for automatic adjustment of the step size used in the solution of the ODE's. This method could be used for many many ODE's without fear of loss of error, since the epsilon for the computer 1.0+Epsilon > 1.0 where Epsilon is as small as possible. The RKF45 method was also included in the "Computer Methods for Mathematical Computations" by Forsythe Malcom and Moler. My point is that you should never be using a straight Runge-Kutta method, but should always be using the Fehlberg variation to improve error analysis and allow your program to adjust the step size, according to the error results by comparing the 4th and 5th order methods. Of course if the objective function is prohibitively expensive, then this method may be expensive, in which case other methods such as collocation, or spectral methods may be better. For a general ODE solver the RKF45 method with step adjustment simply has stood the test of time. Sign In·View Thread·PermaLink Re: I suggest you remind yourself about numerical analysis. headmyshoulder 21:17 12 Nov '09   Error analysis and step-size control is already included in the current version of odeint. Up to now, the RK45 method is missing, but a method which estimates the error by calculating two steps with dt/2 is present. RK45 and RK78 will be included very soon. Nevertheless, in many cases one does not need adaptive step-size control, think for example on creating artificial time series from an ODE. Sign In·View Thread·PermaLink Re: I suggest you remind yourself about numerical analysis. Mario Mulansky 2:18 13 Nov '09   I can assure you that we are aware of the standard numerical methods to solve ode's. The article is focusing more on the structural part of the code and the basic interfaces to some (abstract) step and error adaption routines. We are currently implementing more sophisticated integration methods, like the rk45 with automatic error evaluation (I've just commited the code for that). RK89 will follow as well as Bulirsch-Stoer methods, I think. Sign In·View Thread·PermaLink 5.00/5 (1 vote)

Last Visit: 19:17 9 Feb '10     Last Update: 19:17 9 Feb '10 1