#include <NmOdeStepperRKCK.h>
Inheritance diagram for NmOdeStepperRKCK:
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) |
|
|
Definition at line 22 of file NmStepperRKCK.cxx. 00022 {}
|
|
|
Definition at line 26 of file NmStepperRKCK.cxx. 00026 {}
|
|
||||||||||||||||||||||||||||||||||||
|
Definition at line 115 of file NmStepperRKCK.cxx. Referenced by Step(). 00122 {
00123 //======================================================================
00124 // Purpose: Make a single Cash-Karp Runga-Kutta step
00125 //
00126 // Inputs: y - initial values of dependent variables
00127 // dydx - initial values of derivatives
00128 // n - size of y vector
00129 // x - value of integration variable
00130 // h - step size
00131 // Derivs - functor to use to calculate derivatives
00132 //
00133 // Outputs: yOut - y values after step
00134 // yErr - estimate of erros on yOut
00135 //======================================================================
00136 // These numbers come from Numerical Recipes in C, Table in section 16.2
00137 // "Cash-Karp Parameters for Embedded Runge-Kutta Method"
00138 static const double
00139 a2=1./ 5., a3=3./10., a4=3./5., a5=1., a6=7./8.;
00140 static const double
00141 b21= 1./5.,
00142 b31= 3./40., b32= 9./40.,
00143 b41= 3./10., b42= -9./10., b43= 6./5.,
00144 b51= -11./54., b52= 5./2., b53= -70./27., b54=35./27.,
00145 b61=1631./55296., b62=175./512., b63=575./13824., b64=44275./110592.,
00146 b65=253./4096.;
00147 static const double
00148 c1 = 37./378.,
00149 //unused: c2 = 0.,
00150 c3 = 250./621.,
00151 c4 = 125./594.,
00152 c5 = 0.0,
00153 c6 = 512./1771.;
00154 static const double
00155 dc1 = c1 - 2825./27648.,
00156 //unused: dc2 = c2 - 0.,
00157 dc3 = c3 - 18575./48384.,
00158 dc4 = c4 - 13525./55296.,
00159 dc5 = c5 - 277./14336.,
00160 dc6 = c6 - 1./4.;
00161 register int i;
00162
00163 double* ak2 = new double[n];
00164 double* ak3 = new double[n];
00165 double* ak4 = new double[n];
00166 double* ak5 = new double[n];
00167 double* ak6 = new double[n];
00168 double* yTemp = new double[n];
00169
00170 for (i=0; i<n; ++i) {
00171 yTemp[i] = y[i] + h*b21*dydx[i];
00172 }
00173 Derivs(x + a2*h, yTemp, ak2);
00174 for (i=0; i<n; ++i) {
00175 yTemp[i] = y[i] + h*(b31*dydx[i] + b32*ak2[i]);
00176 }
00177 Derivs(x + a3*h, yTemp, ak3);
00178 for (i=0; i<n; ++i) {
00179 yTemp[i] = y[i] + h*(b41*dydx[i] + b42*ak2[i] + b43*ak3[i]);
00180 }
00181 Derivs(x+a4*h, yTemp, ak4);
00182 for (i=0; i<n; ++i) {
00183 yTemp[i] = y[i] + h*(b51*dydx[i] + b52*ak2[i] + b53*ak3[i] + b54*ak4[i]);
00184 }
00185 Derivs(x + a5*h, yTemp, ak5);
00186 for (i=0; i<n; ++i) {
00187 yTemp[i] = y[i] +
00188 h*(b61*dydx[i] + b62*ak2[i] + b63*ak3[i] + b64*ak4[i] + b65*ak5[i]);
00189 }
00190 Derivs(x+a6*h, yTemp, ak6);
00191
00192 // Accumulate all infomation for the full step
00193 for (i=0; i<n; ++i) {
00194 yOut[i] =
00195 y[i]+h*(c1*dydx[i] + c3*ak3[i] + c4*ak4[i] + c6*ak6[i]);
00196 }
00197
00198 // Estimate the error by comparing 4th and 5th order estimates
00199 for (i=0; i<n; ++i) {
00200 yErr[i] = h*(dc1*dydx[i] + dc3*ak3[i] + dc4*ak4[i] + dc5*ak5[i] +
00201 dc6*ak6[i]);
00202 }
00203 delete []ak2; ak2 = 0;
00204 delete []ak3; ak3 = 0;
00205 delete []ak4; ak4 = 0;
00206 delete []ak5; ak5 = 0;
00207 delete []ak6; ak6 = 0;
00208 delete []yTemp; yTemp = 0;
00209 }
|
|
||||||||||||||||||||||||||||||||||||||||||||
|
Implements NmOdeStepper. Definition at line 30 of file NmStepperRKCK.cxx. References MAXMACRO, MINMACRO, MSG, and RKCKStep(). 00039 {
00040 //======================================================================
00041 // Purpose: Take a single quality-controled step
00042 //
00043 // Inputs: y - initial values of dependent variables. Modified on return
00044 // dydx - initial values of derivatives
00045 // n - size of y and dydx vectors
00046 // x - value of integration variable. Modified on return
00047 // hTry - suggested step size
00048 // eps - fractional accuracy
00049 // yScale - typical sizes of y variables
00050 // Derivs - functor for derivative calculation
00051 // Outputs:
00052 // hDid - size of step taken
00053 // hNext - guess at next good step size
00054 //======================================================================
00055 register int i; // Loop variable
00056 double* yTemp = new double[n]; // tmp values for y variables
00057 double* yErr = new double[n]; // tmp values for calculating errors
00058 double errMax; // Largest parameter error after step
00059 double h; // Step size
00060 // constants which control guess at next good step
00061 static const double gSafety = 0.90;
00062 static const double gPgrow = -0.20;
00063 static const double gPshrink = -0.25;
00064 static const double gMaxGrow = 5.0;
00065 static const double gErrCon = pow(gMaxGrow/gSafety, 1.0/gPgrow);
00066
00067 h = hTry;
00068 for (;;) {
00069 double hTemp; // Trial step size
00070 double xNew; // Next x position
00071
00072 // Make a step using the proposed step size
00073 this->RKCKStep(y, dydx, n, *x, h, yTemp, yErr, Derivs);
00074
00075 // Evaluate estimated sizes of the errors
00076 errMax = 0.0;
00077 for (i=0; i<n; i++) {
00078 if (yScale[i]!=0.0) {
00079 errMax = MAXMACRO(errMax,fabs(yErr[i]/yScale[i]));
00080 }
00081 }
00082 errMax /= eps;
00083
00084 // If errors are acceptable then we're done
00085 if (errMax < 1.0) break;
00086
00087 // If errors are not acceptable shrink the step size and try again
00088 hTemp = gSafety*h*pow(errMax, gPshrink);
00089 h = (h>=0.0 ? MAXMACRO(hTemp, 0.1*h) : MINMACRO(hTemp, 0.1*h));
00090 xNew = (*x)+h;
00091 if (xNew == *x) {
00092 MSG("Nm",Msg::kFatal)
00093 << "Step underflow. eps=" << eps << "\n";
00094 abort();
00095 }
00096 }
00097
00098 // Check if we can safely increase the step size allowing for a maximum
00099 // increase of a factor of gMaxGrow
00100 if (errMax > gErrCon) {
00101 *hNext = gSafety*h*pow(errMax, gPgrow);
00102 }
00103 else {
00104 *hNext = gMaxGrow*h;
00105 }
00106 *x += h;
00107 *hDid = h;
00108 for (i=0; i<n; ++i) y[i] = yTemp[i];
00109 delete []yTemp; yTemp = 0;
00110 delete []yErr; yErr = 0;
00111 }
|
1.3.9.1