NmStepperRKCK.cxx

Go to the documentation of this file.

// $Id: NmStepperRKCK.cxx,v 1.7 2009/02/28 21:46:15 gmieg Exp $
//
// Cash-Karp - Runge-Kutta stepper
//
// messier@huhepl.harvard.edu
#include <cassert>
#include <cmath>

#include "NumericalMethods/NmOdeStepperRKCK.h"
#include "NumericalMethods/NmOdeDerivFunc.h"
#include "MessageService/MsgService.h"
CVSID("$Id: NmStepperRKCK.cxx,v 1.7 2009/02/28 21:46:15 gmieg Exp $");

// MAX and MIN macros
#define MAXMACRO(x,y) (((x)>(y))?(x):(y))
#define MINMACRO(x,y) (((x)<(y))?(x):(y))

//......................................................................

NmOdeStepperRKCK::NmOdeStepperRKCK() {}

//......................................................................

NmOdeStepperRKCK::~NmOdeStepperRKCK() {}

//......................................................................

void NmOdeStepperRKCK::Step(double* y,
                            const double* dydx,
                            int n,
                            double* x, 
                            double hTry, 
                            double eps,
                            const double* yScale,
                            double* hDid,
                            double* hNext, 
                            NmOdeDerivFunc& Derivs) {
//======================================================================
// Purpose: Take a single quality-controled step
//
// Inputs: y - initial values of dependent variables. Modified on return
//         dydx - initial values of derivatives
//            n - size of y and dydx vectors
//            x - value of integration variable. Modified on return
//         hTry - suggested step size
//          eps - fractional accuracy
//       yScale - typical sizes of y variables
//       Derivs - functor for derivative calculation
// Outputs:
//         hDid - size of step taken
//        hNext - guess at next good step size
//======================================================================
  register int i;                // Loop variable
  double* yTemp = new double[n]; // tmp values for y variables
  double* yErr  = new double[n]; // tmp values for calculating errors
  double errMax;                 // Largest parameter error after step
  double h;                      // Step size
  // constants which control guess at next good step
  static const double gSafety  =  0.90;
  static const double gPgrow   = -0.20;
  static const double gPshrink = -0.25;
  static const double gMaxGrow =  5.0;
  static const double gErrCon  = pow(gMaxGrow/gSafety, 1.0/gPgrow);

  h = hTry;
  for (;;) {
    double hTemp; // Trial step size
    double xNew;  // Next x position
    
    // Make a step using the proposed step size
    this->RKCKStep(y, dydx, n, *x, h, yTemp, yErr, Derivs);

    // Evaluate estimated sizes of the errors
    errMax = 0.0;
    for (i=0; i<n; i++) {
      if (yScale[i]!=0.0) {
        errMax = MAXMACRO(errMax,fabs(yErr[i]/yScale[i]));
      }
    }
    errMax /= eps;
    
    // If errors are acceptable then we're done
    if (errMax < 1.0) break; 

    // If errors are not acceptable shrink the step size and try again
    hTemp = gSafety*h*pow(errMax, gPshrink);
    h = (h>=0.0 ? MAXMACRO(hTemp, 0.1*h) : MINMACRO(hTemp, 0.1*h));
    xNew = (*x)+h;
    if (xNew == *x) {
      MSG("Nm",Msg::kFatal)
        << "Step underflow. eps=" << eps << "\n";
      abort();
    }
  }

  // Check if we can safely increase the step size allowing for a maximum
  // increase of a factor of gMaxGrow
  if (errMax > gErrCon) {
    *hNext = gSafety*h*pow(errMax, gPgrow);
  }
  else {
    *hNext = gMaxGrow*h;
  }
  *x += h;
  *hDid = h;
  for (i=0; i<n; ++i) y[i] = yTemp[i];
  delete []yTemp; yTemp = 0;
  delete []yErr;  yErr = 0;
}

//......................................................................

void NmOdeStepperRKCK::RKCKStep(const double* y, 
                                const double* dydx, 
                                int n, 
                                double x, 
                                double h, 
                                double* yOut, 
                                double* yErr, 
                                NmOdeDerivFunc& Derivs) {
//======================================================================
// Purpose: Make a single Cash-Karp Runga-Kutta step
//
// Inputs:    y - initial values of dependent variables
//         dydx - initial values of derivatives
//            n - size of y vector
//            x - value of integration variable
//            h - step size
//       Derivs - functor to use to calculate derivatives
//
// Outputs: yOut - y values after step
//          yErr - estimate of erros on yOut
//======================================================================
  // These numbers come from Numerical Recipes in C, Table in section 16.2
  // "Cash-Karp Parameters for Embedded Runge-Kutta Method"
  static const double 
    a2=1./ 5., a3=3./10., a4=3./5., a5=1., a6=7./8.;
  static const double
    b21= 1./5.,
    b31= 3./40., b32= 9./40.,
    b41= 3./10., b42= -9./10., b43= 6./5.,
    b51= -11./54., b52= 5./2., b53= -70./27., b54=35./27.,
    b61=1631./55296., b62=175./512., b63=575./13824., b64=44275./110592., 
    b65=253./4096.;
  static const double
    c1 = 37./378., 
//unused:    c2 = 0.,
    c3 = 250./621., 
    c4 = 125./594., 
    c5 = 0.0,
    c6 = 512./1771.;
  static const double
    dc1 = c1 - 2825./27648., 
//unused:    dc2 = c2 - 0., 
    dc3 = c3 - 18575./48384.,
    dc4 = c4 - 13525./55296., 
    dc5 = c5 - 277./14336., 
    dc6 = c6 - 1./4.;
  register int i;

  double* ak2   = new double[n];
  double* ak3   = new double[n];
  double* ak4   = new double[n];
  double* ak5   = new double[n];
  double* ak6   = new double[n];
  double* yTemp = new double[n];
  
  for (i=0; i<n; ++i) {
    yTemp[i] = y[i] + h*b21*dydx[i];
  }
  Derivs(x + a2*h, yTemp, ak2);
  for (i=0; i<n; ++i) {
    yTemp[i] = y[i] + h*(b31*dydx[i] + b32*ak2[i]);
  }
  Derivs(x + a3*h, yTemp, ak3);  
  for (i=0; i<n; ++i) {
    yTemp[i] = y[i] + h*(b41*dydx[i] + b42*ak2[i] + b43*ak3[i]);
  }
  Derivs(x+a4*h, yTemp, ak4);
  for (i=0; i<n; ++i) {
    yTemp[i] = y[i] + h*(b51*dydx[i] + b52*ak2[i] + b53*ak3[i] + b54*ak4[i]);
  }
  Derivs(x + a5*h, yTemp, ak5);
  for (i=0; i<n; ++i) {
    yTemp[i] = y[i] + 
      h*(b61*dydx[i] + b62*ak2[i] + b63*ak3[i] + b64*ak4[i] + b65*ak5[i]);
  }
  Derivs(x+a6*h, yTemp, ak6);

  // Accumulate all infomation for the full step
  for (i=0; i<n; ++i) {
    yOut[i] = 
      y[i]+h*(c1*dydx[i] + c3*ak3[i] + c4*ak4[i] + c6*ak6[i]);
  }
  
  // Estimate the error by comparing 4th and 5th order estimates
  for (i=0; i<n; ++i) {
    yErr[i] = h*(dc1*dydx[i] + dc3*ak3[i] + dc4*ak4[i] + dc5*ak5[i] + 
                 dc6*ak6[i]);
  }
  delete []ak2;   ak2   = 0;
  delete []ak3;   ak3   = 0;
  delete []ak4;   ak4   = 0;
  delete []ak5;   ak5   = 0;
  delete []ak6;   ak6   = 0;
  delete []yTemp; yTemp = 0;
}

---