NmOdeStepperRKCK Class Reference

#include <NmOdeStepperRKCK.h>

Inheritance diagram for NmOdeStepperRKCK:
NmOdeStepper

List of all members.

Public Member Functions

 NmOdeStepperRKCK ()
 ~NmOdeStepperRKCK ()
void Step (double *y, const double *dydx, int n, double *x, double hTry, double eps, const double *yScale, double *hDid, double *hNext, NmOdeDerivFunc &Deriv)

Protected Member Functions

void RKCKStep (const double *y, const double *dydx, int n, double x, double h, double *yOut, double *yErr, NmOdeDerivFunc &Deriv)

Detailed Description

Definition at line 15 of file NmOdeStepperRKCK.h.


Constructor & Destructor Documentation

NmOdeStepperRKCK::NmOdeStepperRKCK (  ) 

Definition at line 23 of file NmStepperRKCK.cxx.

00023 {}

NmOdeStepperRKCK::~NmOdeStepperRKCK (  ) 

Definition at line 27 of file NmStepperRKCK.cxx.

00027 {}


Member Function Documentation

void NmOdeStepperRKCK::RKCKStep ( const double *  y,
const double *  dydx,
int  n,
double  x,
double  h,
double *  yOut,
double *  yErr,
NmOdeDerivFunc Deriv 
) [protected]

Definition at line 116 of file NmStepperRKCK.cxx.

Referenced by Step().

00123                                                         {
00124 //======================================================================
00125 // Purpose: Make a single Cash-Karp Runga-Kutta step
00126 //
00127 // Inputs:    y - initial values of dependent variables
00128 //         dydx - initial values of derivatives
00129 //            n - size of y vector
00130 //            x - value of integration variable
00131 //            h - step size
00132 //       Derivs - functor to use to calculate derivatives
00133 //
00134 // Outputs: yOut - y values after step
00135 //          yErr - estimate of erros on yOut
00136 //======================================================================
00137   // These numbers come from Numerical Recipes in C, Table in section 16.2
00138   // "Cash-Karp Parameters for Embedded Runge-Kutta Method"
00139   static const double 
00140     a2=1./ 5., a3=3./10., a4=3./5., a5=1., a6=7./8.;
00141   static const double
00142     b21= 1./5.,
00143     b31= 3./40., b32= 9./40.,
00144     b41= 3./10., b42= -9./10., b43= 6./5.,
00145     b51= -11./54., b52= 5./2., b53= -70./27., b54=35./27.,
00146     b61=1631./55296., b62=175./512., b63=575./13824., b64=44275./110592., 
00147     b65=253./4096.;
00148   static const double
00149     c1 = 37./378., 
00150 //unused:    c2 = 0.,
00151     c3 = 250./621., 
00152     c4 = 125./594., 
00153     c5 = 0.0,
00154     c6 = 512./1771.;
00155   static const double
00156     dc1 = c1 - 2825./27648., 
00157 //unused:    dc2 = c2 - 0., 
00158     dc3 = c3 - 18575./48384.,
00159     dc4 = c4 - 13525./55296., 
00160     dc5 = c5 - 277./14336., 
00161     dc6 = c6 - 1./4.;
00162   register int i;
00163 
00164   double* ak2   = new double[n];
00165   double* ak3   = new double[n];
00166   double* ak4   = new double[n];
00167   double* ak5   = new double[n];
00168   double* ak6   = new double[n];
00169   double* yTemp = new double[n];
00170   
00171   for (i=0; i<n; ++i) {
00172     yTemp[i] = y[i] + h*b21*dydx[i];
00173   }
00174   Derivs(x + a2*h, yTemp, ak2);
00175   for (i=0; i<n; ++i) {
00176     yTemp[i] = y[i] + h*(b31*dydx[i] + b32*ak2[i]);
00177   }
00178   Derivs(x + a3*h, yTemp, ak3);  
00179   for (i=0; i<n; ++i) {
00180     yTemp[i] = y[i] + h*(b41*dydx[i] + b42*ak2[i] + b43*ak3[i]);
00181   }
00182   Derivs(x+a4*h, yTemp, ak4);
00183   for (i=0; i<n; ++i) {
00184     yTemp[i] = y[i] + h*(b51*dydx[i] + b52*ak2[i] + b53*ak3[i] + b54*ak4[i]);
00185   }
00186   Derivs(x + a5*h, yTemp, ak5);
00187   for (i=0; i<n; ++i) {
00188     yTemp[i] = y[i] + 
00189       h*(b61*dydx[i] + b62*ak2[i] + b63*ak3[i] + b64*ak4[i] + b65*ak5[i]);
00190   }
00191   Derivs(x+a6*h, yTemp, ak6);
00192 
00193   // Accumulate all infomation for the full step
00194   for (i=0; i<n; ++i) {
00195     yOut[i] = 
00196       y[i]+h*(c1*dydx[i] + c3*ak3[i] + c4*ak4[i] + c6*ak6[i]);
00197   }
00198   
00199   // Estimate the error by comparing 4th and 5th order estimates
00200   for (i=0; i<n; ++i) {
00201     yErr[i] = h*(dc1*dydx[i] + dc3*ak3[i] + dc4*ak4[i] + dc5*ak5[i] + 
00202                  dc6*ak6[i]);
00203   }
00204   delete []ak2;   ak2   = 0;
00205   delete []ak3;   ak3   = 0;
00206   delete []ak4;   ak4   = 0;
00207   delete []ak5;   ak5   = 0;
00208   delete []ak6;   ak6   = 0;
00209   delete []yTemp; yTemp = 0;
00210 }

void NmOdeStepperRKCK::Step ( double *  y,
const double *  dydx,
int  n,
double *  x,
double  hTry,
double  eps,
const double *  yScale,
double *  hDid,
double *  hNext,
NmOdeDerivFunc Deriv 
) [virtual]

Implements NmOdeStepper.

Definition at line 31 of file NmStepperRKCK.cxx.

References Msg::kFatal, MAXMACRO, MINMACRO, MSG, and RKCKStep().

00040                                                     {
00041 //======================================================================
00042 // Purpose: Take a single quality-controled step
00043 //
00044 // Inputs: y - initial values of dependent variables. Modified on return
00045 //         dydx - initial values of derivatives
00046 //            n - size of y and dydx vectors
00047 //            x - value of integration variable. Modified on return
00048 //         hTry - suggested step size
00049 //          eps - fractional accuracy
00050 //       yScale - typical sizes of y variables
00051 //       Derivs - functor for derivative calculation
00052 // Outputs:
00053 //         hDid - size of step taken
00054 //        hNext - guess at next good step size
00055 //======================================================================
00056   register int i;                // Loop variable
00057   double* yTemp = new double[n]; // tmp values for y variables
00058   double* yErr  = new double[n]; // tmp values for calculating errors
00059   double errMax;                 // Largest parameter error after step
00060   double h;                      // Step size
00061   // constants which control guess at next good step
00062   static const double gSafety  =  0.90;
00063   static const double gPgrow   = -0.20;
00064   static const double gPshrink = -0.25;
00065   static const double gMaxGrow =  5.0;
00066   static const double gErrCon  = pow(gMaxGrow/gSafety, 1.0/gPgrow);
00067 
00068   h = hTry;
00069   for (;;) {
00070     double hTemp; // Trial step size
00071     double xNew;  // Next x position
00072     
00073     // Make a step using the proposed step size
00074     this->RKCKStep(y, dydx, n, *x, h, yTemp, yErr, Derivs);
00075 
00076     // Evaluate estimated sizes of the errors
00077     errMax = 0.0;
00078     for (i=0; i<n; i++) {
00079       if (yScale[i]!=0.0) {
00080         errMax = MAXMACRO(errMax,fabs(yErr[i]/yScale[i]));
00081       }
00082     }
00083     errMax /= eps;
00084     
00085     // If errors are acceptable then we're done
00086     if (errMax < 1.0) break; 
00087 
00088     // If errors are not acceptable shrink the step size and try again
00089     hTemp = gSafety*h*pow(errMax, gPshrink);
00090     h = (h>=0.0 ? MAXMACRO(hTemp, 0.1*h) : MINMACRO(hTemp, 0.1*h));
00091     xNew = (*x)+h;
00092     if (xNew == *x) {
00093       MSG("Nm",Msg::kFatal)
00094         << "Step underflow. eps=" << eps << "\n";
00095       abort();
00096     }
00097   }
00098 
00099   // Check if we can safely increase the step size allowing for a maximum
00100   // increase of a factor of gMaxGrow
00101   if (errMax > gErrCon) {
00102     *hNext = gSafety*h*pow(errMax, gPgrow);
00103   }
00104   else {
00105     *hNext = gMaxGrow*h;
00106   }
00107   *x += h;
00108   *hDid = h;
00109   for (i=0; i<n; ++i) y[i] = yTemp[i];
00110   delete []yTemp; yTemp = 0;
00111   delete []yErr;  yErr = 0;
00112 }


The documentation for this class was generated from the following files:

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1