odeint v2 - Solving ordinary differential equations in C++

/*
 * nr_rk4_phase_lattice.cpp
 *
 *  Created on: May 15, 2011
 *      Author: mario
 */

#include <boost/array.hpp>

#include "rk_performance_test_case.hpp"

#include "phase_lattice.hpp"

const size_t dim = 1024;

typedef boost::array< double , dim > state_type;


template< class System , typename T , size_t dim >
void rk4_step( const System sys , boost::array< T , dim > &x , const double t , const double dt )
{   // fast rk4 implementation adapted from the book 'Numerical Recipes'
    size_t i;
    const double hh = dt*0.5;
    const double h6 = dt/6.0;
    const double th = t+hh;
    boost::array< T , dim > dydx , dym , dyt , yt;

    sys( x , dydx , t );

    for( i=0 ; i<dim ; i++ )
        yt[i] = x[i] + hh*dydx[i];

    sys( yt , dyt , th );
    for( i=0 ; i<dim ; i++ )
        yt[i] = x[i] + hh*dyt[i];

    sys( yt , dym , th );
    for( i=0 ; i<dim ; i++ ) {
        yt[i] = x[i] + dt*dym[i];
        dym[i] += dyt[i];
    }
    sys( yt , dyt , t+dt );
    for( i=0 ; i<dim ; i++ )
        x[i] += h6*( dydx[i] + dyt[i] + 2.0*dym[i] );
}


class nr_wrapper
{
public:
    void reset_init_cond()
    {
        for( size_t i = 0 ; i<dim ; ++i )
			m_x[i] = 2.0*3.1415927*rand() / RAND_MAX;
		m_t = 0.0;
    }

    inline void do_step( const double dt )
    {
        rk4_step( phase_lattice<dim>() , m_x , m_t , dt );
    }

    double state( const size_t i ) const
    { return m_x[i]; }

private:
    state_type m_x;
    double m_t;
};



int main()
{
	srand( 12312354 );

    nr_wrapper stepper;

    run( stepper , 10000 , 1E-6 );
}