Re: [boost] [odeint] gsoc 2011: finalization of odeint

Subject: Re: [boost] [odeint] gsoc 2011: finalization of odeint
From: Matthias Schabel (boost_at_[hidden])
Date: 2011-03-22 15:47:27

> Just a little remark,
> 
> the current (development) version with all the features Mario described
> is in
> 
> https://svn.boost.org/svn/boost/sandbox/odeint/branches/karsten
> 
> Some Runge-Kutta steppers are missing but it is easy to add them.

I would love to see this in Boost. Of course, odeint needs to play well with Boost.Units. Here is how I would expect the code to look (stripped down version of numeric/odeint/examples/lorenz_integrator.cpp) :

#include <iostream>
#include <iterator>

#include <boost/tr1/array.hpp>

#include <boost/units/systems/si.hpp>

#include <boost/numeric/odeint.hpp>
#include <boost/numeric/odeint/container_traits_tr1_array.hpp>

#define tab "\t"

using namespace std;
using namespace boost;
using namespace boost::units;
using namespace boost::numeric::odeint;

typedef quantity<si::time>				time_type;
typedef array<quantity<si::length>, 3>	position_type;
typedef array<quantity<si::velocity>, 3>	velocity_type;

const double eps_abs = 1E-6;
const double eps_rel = 1E-7;

const size_t time_points = 101;

typedef velocity_type (*function_type)(const position_type&,const time_type&);

velocity_type lorenz(const position_type &x,const time_type& t )
{
        const quantity<si::frequency>	sigma(10.0*si::hertz);
        const quantity<si::length>	        rho(28.0*si::meter),
                                                                beta((8.0/3.0)*si::meter);
        
    const array<quantity<si::velocity>,3>	dxdt = { 	(x[1] - x[0])*sigma,
                                                                                         (x[0]*(rho - x[2])/(1.0*si::meter) - x[1])/si::second,
                                                                                         (x[0]*x[1] - beta * x[2])/si::meter/si::second };
                                                                                                         
        return dxdt;
}

int main( int argc , char **argv )
{
    position_type x1;

    x1[0] = 1.0*si::meters;
    x1[1] = 0.0*si::meters;
    x1[2] = 20.0*si::meters;

    vector<position_type>	x_vec;
    vector<time_type>		t_vec;

    stepper_half_step< stepper_euler< function_type > > euler(lorenz);

    size_t steps = integrate( euler, 
                                                x1, 
                                                      0.0*si::seconds, 
                                                10.0*si::seconds, 
                                                1E-4*si::seconds, 
                                                      back_inserter(t_vec),
                                                      back_inserter(x_vec));

    clog << "Euler Half Step: #steps " << steps << endl;

    cout.precision(16);
    cout.setf(ios::fixed,ios::floatfield);

    cout << t_vec.size() << endl;
    
    cout << (x_vec.back())[0] << tab << (x_vec.back())[1] << tab << (x_vec.back())[2] << tab;

    return 0;
}

In particular, the steppers should be able to take a function pointer or other function object and derive the container types directly from it rather than requiring that these be provided as template parameters. Then the user need only provide the function to integrate (optionally providing a Jacobian for stiff solvers).

Matthias