ode_rkck.cc

Go to the documentation of this file.

 
 
#include "ode_rkck.h"
 
namespace ode {
 
OdeRKCK::OdeRKCK (unsigned long neq) :
    OdeEmbedded (neq, false, 4),
    OdeRK (neq, 6),
    OdeERK (neq) {
 
    method_ = "RKCK";
    //coefficents of tableau
    c2 =  1.0/5; a21 =        1.0/5;
    c3 = 3.0/10; a31 =       3.0/40; a32 =    9.0/40;
    c4 =  3.0/5; a41 =       3.0/10; a42 =   -9.0/10; a43 =        6.0/5;
    c5 =    1.0; a51 =     -11.0/54; a52 =     5.0/2; a53 =     -70.0/27; a54 = 35.0/27;
    c6 =  7.0/8; a61 = 1631.0/55296; a62 = 175.0/512; a63 =  575.0/13824; a64 = 44275.0/110592; a65 = 253.0/4096;
                  b1 =     37.0/378;                  b3 =     250.0/621; b4  =      125.0/594;                   b6 = 512.0/1771;
                  d1 = 2825.0/27648;                  d3 = 18575.0/48384; d4 =   13525.0/55296; d5 = 277.0/14336; d6 = 1.0/4;
                  e1 =      19.0/54;                  e3 =      -10.0/27; e4 =         55.0/54;
                  f1 =       -3.0/2; f2 =      5.0/2;
                  g1 =          1.0;
}
 
void OdeRKCK::step_ (double dt) {
 
    unsigned long i;
 
    //------------------------------------------------------------------
    ode_fun_(sol_, k_[0]);
 
    //------------------------------------------------------------------
    for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*a21*k_[0][i];
    ode_fun_(soltemp_, k_[1]);
 
    //------------------------------------------------------------------
    for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a31*k_[0][i]
                                                     + a32*k_[1][i]);
    ode_fun_(soltemp_, k_[2]);
 
    //------------------------------------------------------------------
    for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a41*k_[0][i]
                                                     + a42*k_[1][i]
                                                     + a43*k_[2][i]);
    ode_fun_(soltemp_, k_[3]);
 
    //------------------------------------------------------------------
    for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a51*k_[0][i]
                                                     + a52*k_[1][i]
                                                     + a53*k_[2][i]
                                                     + a54*k_[3][i]);
    ode_fun_(soltemp_, k_[4]);
 
    //------------------------------------------------------------------
    for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a61*k_[0][i]
                                                     + a62*k_[1][i]
                                                     + a63*k_[2][i]
                                                     + a64*k_[3][i]
                                                     + a65*k_[4][i]);
    ode_fun_(soltemp_, k_[5]);
 
    //------------------------------------------------------------------
    for (i=0; i<neq_; i++) {
        solemb_[i] = sol_[i] + dt*(d1*k_[0][i]
                                 + d3*k_[2][i]
                                 + d4*k_[3][i]
                                 + d5*k_[4][i]
                                 + d6*k_[5][i]);
        sol_[i] = sol_[i] + dt*(b1*k_[0][i]
                              + b3*k_[2][i]
                              + b4*k_[3][i]
                              + b6*k_[5][i]);
    }
}
 
} // namespace ode