Changeset 6936 for lang/cplusplus/boost-supplement

Timestamp:

02/19/08 20:41:04 (5 years ago)

Author:

mrkn

Message:

lang/cplusplus/boost-supplement/branches/numeric-ode/: libs/numeric/ode/test_runge_kutta_rk4.cpp (nr_rkck, nr_rkqs, nr_odeint): These functions are implemented based onf Numerical Recipes in C, 2nd ed. and nr_rkqc has been needed.

Location:

lang/cplusplus/boost-supplement/branches/numeric-ode

Files:

• ChangeLog (modified) (1 diff)
• libs/numeric/ode/test_runge_kutta_rk4.cpp (modified) (3 diffs)

lang/cplusplus/boost-supplement/branches/numeric-ode/ChangeLog

@@ +1 @@
+2008-02-19  Kenta Murata  <mrkn@gmail.com>
+
+        * libs/numeric/ode/test_runge_kutta_rk4.cpp
+        (nr_rkck, nr_rkqs, nr_odeint): These functions are implemented
+        based onf Numerical Recipes in C, 2nd ed. and nr_rkqc has been
+        needed.
+
 2008-02-18  Kenta Murata  <mrkn@gmail.com>
 

lang/cplusplus/boost-supplement/branches/numeric-ode/libs/numeric/ode/test_runge_kutta_rk4.cpp

@@ -8 +8 @@
 #include <boost/test/floating_point_comparison.hpp>
 #include <boost/tuple/tuple.hpp>
+
+#include <stdexcept>
 
 typedef boost::numeric::ublas::vector<double> dvector;
@@ -67 +69 @@
     nr_rk4(v, dv, nvar, x, h, vout, derivs);
     if ((double)(x+h) == x)
-      throw "Step size too small in routine rkdumb";
+      throw std::runtime_error("Step size too small in routine rkdumb");
     x += h;
     xx[k+1] = x;
@@ -82 +84 @@
 
 void
+nr_rkck(double y[], double dydx[], int n, double x, double h,
+        double yout[], double yerr[],
+        void (*derivs)(double, double[], double[]))
+// Given value for n variables y[1..n] and their derivatives
+// dydx[1..n] known at x, use the fifth-order Cash-Karp Runge-Kutta
+// method to advance the solution over an interval h and return the
+// incremented variables as yout[1..n].  Also return an estimate of
+// the local truncation error in yout using the embedded fourth-order
+// method.  The user supplies the routine derivs(x, y, dydx), which
+// returns derivatives dydx at x.
+{
+  int 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,
+    dc5 = -277.0/14336.0;
+  static cdonst double dc1=c1-2825.0/27648.0, dc3=c3-18575.0/48384.0,
+    dc4=c4-13525.0/55296.0, dc6=c6-0.25;
+  double *ak2, *ak3, *ak4, *ak5, *ak6, *ytemp;
+
+  ak2 = new double[n+1];
+  ak3 = new double[n+1];
+  ak4 = new double[n+1];
+  ak5 = new double[n+1];
+  ak6 = new double[n+1];
+  ytemp = new double[n+1];
+
+  for (i = 1; i <= n; i++)      // First step.
+    ytemp[i]=y[i]+b21*h*dydx[i];
+
+  (*derivs)(x+a2*h, ytemp, ak2); // Second step.
+  for (i = 1; i <= n; i++)
+    ytemp[i]=y[i]+h*(b31*dydx[i]+b32*ak2[i]);
+
+  (*derivs)(x+a3*h, ytemp, ak3); // Third step.
+  for (i = 1; i <= n; i++)
+    ytemp[i]=y[i]+h*(b41*dydx[i]+b42*ak2[i]+b43*ak3[i]);
+
+  (*derivs)(x+a4*h, ytemp, ak4); // Fourth step.
+  for (i = 1; 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); // Fifth step.
+  for (i = 1; 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); // Sixth step.
+  for (i = 1; i <= n; i++) // Accumulate increments with proper weights.
+    yout[i]=y[i]+h*(c1*dydx[i]+c3*ak3[i]+c4*ak4[i]+c6*ak6[i]);
+  for (i = 1; i <= n; i++)
+    // Estimate error as difference between fourth and fifth order methods.
+    yerr[i]=h*(dc1*dydx[i]+dc3*ak3[i]+dc4*ak4[i]+dc5*ak5[i]+dc6*ak6[i]);
+
+  delete[] ytemp;
+  delete[] ak6;
+  delete[] ak5;
+  delete[] ak4;
+  delete[] ak3;
+  delete[] ak2;
+}
+
+void
 nr_rkqc(double y[], double dydx[], int n, double* x, double htry,
         double eps, double yscal[], double* hdid, double* hnext,
-        void (*derivs)(double, double*, double*))
-{
+        void (*derivs)(double, double[], double[]))
+// Fifth-order Runge-Kutta step with monitoring of local truncation
+// error to ensure accuracy and adjust stepsize.  Input are the
+// dependent variable vector y[1..n] and its derivative dydx[1..n] at
+// the starting value of the independent variable x.  Also input are
+// the stepsize to be attempted htry, the required accuracy eps, and
+// the vector yscal[1..n] against which the error is scaled.  On
+// output, y and x are replaced by their new values, hdid is the
+// stepsize that was actually accomplished, and hnext is the estimated
+// next stepsize.  derivs is the user-supplied routine that computes
+// the right-hand side derivatives.
+{
+  static const double SAFETY = 0.9;
   static const double PGROW  = -0.20;
   static const double PSHRNK = -0.25;
-  static const double FCOR   = 0.06666666; // 1.0 / 15.0
-  static const double SAFETY = 0.9;
-  static const double ERRCON = 6.0e-4; // == (4.0/SAFETY)**(1/PGROW)
+
+  // The value ERRCON equals (5/SAFETY) raised to the power (1/PGROW).
+  static const double ERRCON = 1.89e-4;
 
   int i;
-  double xsav, hh, h, temp, errmax;
-  double *dysav, *ysav, *ytemp;
-
-  dysav = new double[n+1];
-  ysav  = new double[n+1];
+  double errmax, h, htemp, xnew, *yerr, *ytemp;
+
+  yerr = new double[n+1];
   ytemp = new double[n+1];
-
-  xsav = *x;
-  for (i = 1; i <= n; ++i) {
-    ysav[i] = y[i];
-    dysav[i] = dydx[i];
-  }
-  h = htry;
+  h = htry;                // Set stepsize to the initial trial value.
   for (;;) {
-    // two half-steps
-    hh = 0.5*h;
-    rk4(ysav, dysav, n, xsav, hh, ytemp, derivs);
-    *x = xsav + hh;
-    (*derivs)(*x, ytemp, dydx);
-    rk4(ytemp, dydx, n, *x, hh, y, derivs);
-
-    // the full-step
-    *x = xsav + h;
-    if (*x == xsav)
-      throw "Step size too small in routine rkqc";
-    rk4(ysav, dysav, n, xsav, h, ytemp, derivs);
-
-    // calculate precision
-    errmax = 0.0;
-    for (i = 1; i <= n; ++i) {
-      ytemp[i] = y[i] - ytemp[i];
-      temp = fabs(ytemp[i] / yscal[i]);
-      if (errmax < temp) errmax = temp;
+    nr_rkck(y, dydx, n, *x, h, ytemp, yerr, derivs); // Take a step.
+    errmax = 0.0;               // Evaluate accuracy.
+    for (i = 1; i <= n; i++)
+      errmax = FMAX(errmax, fabs(yerr[i]/yscal[i]));
+    errmax /= eps;            // Scale relative to required tolerance.
+    if (errmax <= 1.0) break; // Step succeeded. Compute size of next step.
+    htemp = SAFETY*h*pow(errmax, PSHRNK);
+    // Truncation error too large, reduce stepsize.
+    h = (h >= 0.0 ? FMAX(htemp, 0.1*h) : FMIN(htemp, 0.1*h));
+    // No more than a factor of 10.
+    xnew = (*x)+h;
+    if (xnew == *x)
+      throw std::runtime_error("stepsize underflow in nr_rkqs");
+  }
+
+  if (errmax > ERRCON) *hnext = SAFETY*h*pow(errmax, PGROW);
+  else *hnext = 5.0*h;         // No more than a factor of 5 increase.
+  *x += (*hdid=h);
+  for (i = 1; i <= n; i++) y[i] = ytemp[i];
+
+  delete[] ytemp;
+  delete[] yerr;
+}
+
+const int MAXSTP = 10000;
+const double TINY = 1.e-30;
+
+// User storage for intermediate results.  Preset kmax and dxsav in
+// the calling program.  If kmax != 0 results are stored at
+// approximate intervals dxsav in the arrays xp[1..kount],
+// yp[1..nvar][1..kount], where kount is output by odeint.  Defining
+// declarations for these variables, with memory allocations
+// xp[1..kmax] and yp[1..nvar][1..kmax] for the arrays, should be in
+// the calling program.
+struct nr_odeint_data 
+{
+  int kmax, kount;
+  double *xp, **yp, dxsav;
+};
+
+void
+nr_odeint(double ystart[], int nvar, double x1, double x2, double eps,
+          double h1, double hmin, int* nok, int* nbad, nr_odeint_data* data,
+          void (*derivs)(double, double[], double[]),
+          void (*rkqs)(double[], double[], int, double*, double,
+                       double, double[], double*, double*,
+                       void (*)(double, double[], double[])))
+// Runge-Kutta driver with adaptive stepsize control.  Integrate
+// starting values ystart[1..nvar] from x1 to x2 with accuracy eps,
+// storing intermediate results into the variable data.  h1 should be
+// set as a guessed first stepsize, hmin as the minimum allowed
+// stepsize (can be zero).  On output nok and nbad are the number of
+// good and bad (but retried and fixed) steps taken, and ystart is
+// replaced by values at the end of the integration interval.  derivs
+// is the user-supplied routine for calculating the right-hand side
+// derivative, while rkqs is the name of the stepper routine to be
+// used.
+{
+  int nstp, i;
+  double xsav, x, hnext, hdid, h;
+  double *yscal, *y, *dydx;
+
+  yscal = new double[nvar+1];
+  y = new double[nvar+1];
+  dydx = new double[nvar+1];
+
+  x = x1;
+  h = SIGN(h1, x2-x1);
+  *nok = (*nbad) = data->kount = 0;
+  for (i = 1; i <= nvar; ++i) y[i] = ystart[i];
+  if (data->kmax > 0)
+    xsav = x - data->dxsav*2.0; // Assures storage of first step.
+  for (nstp = 1; nstp <= MAXSTP; nstp++) { // Take at most MAXSTP steps.
+    (*derivs)(x, y, dydx);
+    for (i = 1; i <= nvar; i++)
+      // Scaling used to monitor accuracy.  This general-purpose
+      // choice can be modified if need be.
+      yscal[i] = fabs(y[i])+fabs(dydx[i]*h)+TINY;
+    if (data->kmax > 0 && data->kount < data->kmax-1 &&
+        fabs(x-xsav) > fabs(data->dxsav)) {
+      data->xp[++data->kount] = x;          // Store intermediate results.
+      for (i = 1; i <= nvar; i++)
+        data->yp[i][data->kount] = y[i];
+      xsav = x;
     }
-    errmax /= eps;              // scale for allowing error
-    if (errmax <= 1.0) {
-      *hdid = h;
-      *hnext = (errmax > ERRCON ?
-                SAFETY*h*exp(PGROW*log(errmax)) : 4.0*h);
-      break;
+    if ((x + h - x2)*(x+h-x1) > 0.0) h = x2 - x; // If stepsize can
+                                                 // overshoot,
+                                                 // decrese.
+    (*rkqs)(y, dydx, nvar, &x, h, eps, yscal, &hdid, &hnext, derivs);
+    if (hdid == h) ++(*nok); else ++(*nbad);
+    if ((x-x2)*(x2-x1) >= 0.0) { // Are we done?
+      for (i = 1; i <= nvar; i++) ystart[i] = y[i];
+      if (data->kmax) {
+        data->xp[++data->kount] = x; // Save final step.
+        for (i = 1; i <= nvar; i++)
+          data->yp[i][data->kount] = y[i];
+      }
+      delete[] dydx;
+      delete[] y;
+      delete[] yscal;
+      return;                   // Normal exit.
     }
-    h = SAFETY*h*exp(PSHRNK*log(errmax));
-  }
-  for (i = 1; i <= n; ++i)
-    y[i] += ytemp[i]*FCOR;
-
-  delete[] ytemp;
-  delete[] ysav;
-  delete[] dysav;
-}
-
+    if (fabs(hnext) <= hmin)
+      throw std::runtime_error("Step size too small in nr_odeint");
+    h = hnext;
+  }
+  throw std::runtime_error("Too many steps in routine nr_odeint");
 }