/*
* nr_rk54ck_lorenz.cpp
*
* Created on: May 12, 2011
* Author: mario
*/
#include <boost/array.hpp>
#include "rk_performance_test_case.hpp"
#include "lorenz.hpp"
const size_t dim = 3;
typedef boost::array< double , dim > state_type;
template< class System , typename T , size_t dim >
void rk54ck_step( System sys , boost::array< T , dim > &x , const double t , const double dt , boost::array< T , dim > &xerr )
{ // fast rk54ck implementation adapted from the book 'Numerical Recipes'
size_t i;
static const double a2 = 0.2 , a3 = 0.3 , a4 = 0.6 , a5 = 1.0 , a6 = 0.875 ,
b21 = 0.2 , b31 = 3.0 / 40.0 , b32 = 9.0 / 40.0 , b41 = 0.3 , b42 =
-0.9 , b43 = 1.2 , b51 = -11.0 / 54.0 , b52 = 2.5 , b53 =
-70.0 / 27.0 , b54 = 35.0 / 27.0 , b61 = 1631.0 / 55296.0 ,
b62 = 175.0 / 512.0 , b63 = 575.0 / 13824.0 , b64 = 44275.0
/ 110592.0 , b65 = 253.0 / 4096.0 , c1 = 37.0 / 378.0 , c3 =
250.0 / 621.0 , c4 = 125.0 / 594.0 , c6 = 512.0 / 1771.0 ,
dc1 = c1 - 2825.0 / 27648.0 , dc3 = c3 - 18575.0 / 48384.0 , dc4 =
c4 - 13525.0 / 55296.0 , dc5 = -277.00 / 14336.0 , dc6 = c6
- 0.25;
const size_t n = dim;
boost::array< T , dim > dydx , ak2 , ak3 , ak4 , ak5 , ak6 , ytemp;
sys( x , dydx , t );
for (i=0;i<n;i++)
ytemp[i] = x[i] + b21 * dt * dydx[i];
sys( ytemp , ak2 , t+a2*dt );
for (i=0;i<n;i++)
ytemp[i] = x[i] + dt*(b31*dydx[i]+b32*ak2[i]);
sys( ytemp , ak3 , t+a3*dt );
for (i=0;i<n;i++)
ytemp[i] = x[i] + dt*(b41*dydx[i]+b42*ak2[i]+b43*ak3[i]);
sys( ytemp , ak4 , t+a4*dt );
for (i=0;i<n;i++)
ytemp[i] = x[i] + dt*(b51*dydx[i]+b52*ak2[i]+b53*ak3[i]+b54*ak4[i]);
sys( ytemp, ak5 , t+a5*dt );
for (i=0;i<n;i++)
ytemp[i] = x[i] + dt*(b61*dydx[i]+b62*ak2[i]+b63*ak3[i]+b64*ak4[i]+b65*ak5[i]);
sys( ytemp , ak6 , t+a6*dt );
for (i=0;i<n;i++)
x[i] += dt*(c1*dydx[i]+c3*ak3[i]+c4*ak4[i]+c6*ak6[i]);
for (i=0;i<n;i++)
xerr[i] = dt*(dc1*dydx[i]+dc3*ak3[i]+dc4*ak4[i]+dc5*ak5[i]+dc6*ak6[i]);
}
class nr_wrapper
{
public:
void reset_init_cond()
{
m_x[0] = 10.0 * rand() / RAND_MAX;
m_x[1] = 10.0 * rand() / RAND_MAX;
m_x[2] = 10.0 * rand() / RAND_MAX;
m_t = 0.0;
}
inline void do_step( const double dt )
{
rk54ck_step( lorenz() , m_x , m_t , dt , m_x_err );
}
double state( const size_t i ) const
{ return m_x[i]; }
private:
state_type m_x;
state_type m_x_err;
double m_t;
};
int main()
{
nr_wrapper stepper;
run( stepper );
}