RungeKutta/src/RK4PreciseStep.cc

0001 #include "RK4PreciseStep.h"
0002 #include "RK4OneStep.h"
0003 //#include "Utilities/UI/interface/SimpleConfigurable.h"
0004 #include <iostream>
0005
0006 CartesianState RK4PreciseStep::operator()(const CartesianState& start,
0007                                           const RKCartesianDerivative& deriv,
0008                                           double step,
0009                                           double eps) const {
0010   const double Safety = 0.9;
0011   double remainigStep = step;
0012   double stepSize = step;
0013   CartesianState currentStart = start;
0014   int nsteps = 0;
0015   std::pair<CartesianState, double> tryStep;
0016
0017   do {
0018     tryStep = stepWithAccuracy(currentStart, deriv, stepSize);
0019     nsteps++;
0020     if (tryStep.second < eps) {
0021       if (remainigStep - stepSize < eps / 2) {
0022         if (verbose())
0023           std::cout << "Accuracy reached, and full step taken in " << nsteps << " steps" << std::endl;
0024         return tryStep.first;  // we are there
0025       } else {
0026         remainigStep -= stepSize;
0027         // increase step size
0028         double factor = std::min(Safety * pow(fabs(eps / tryStep.second), 0.2), 4.);
0029         stepSize = std::min(stepSize * factor, remainigStep);
0030         currentStart = tryStep.first;
0031         if (verbose())
0032           std::cout << "Accuracy reached, but " << remainigStep << " remain after " << nsteps
0033                     << " steps. Step size increased by " << factor << " to " << stepSize << std::endl;
0034       }
0035     } else {
0036       // decrease step size
0037       double factor = std::max(Safety * pow(fabs(eps / tryStep.second), 0.25), 0.1);
0038       stepSize *= factor;
0039       if (verbose())
0040         std::cout << "Accuracy not yet reached: delta = " << tryStep.second << ", step reduced by " << factor << " to "
0041                   << stepSize << ", (R,z)= " << currentStart.position().perp() << ", " << currentStart.position().z()
0042                   << std::endl;
0043     }
0044   } while (remainigStep > eps / 2);
0045
0046   return tryStep.first;
0047 }
0048
0049 std::pair<CartesianState, double> RK4PreciseStep::stepWithAccuracy(const CartesianState& start,
0050                                                                    const RKCartesianDerivative& deriv,
0051                                                                    double step) const {
0052   RK4OneStep solver;
0053   CartesianState one(solver(start, deriv, step));
0054   CartesianState firstHalf(solver(start, deriv, step / 2));
0055   CartesianState secondHalf(solver(firstHalf, deriv, step / 2));
0056   double diff = distance(one, secondHalf);
0057   return std::pair<CartesianState, double>(secondHalf, diff);
0058 }
0059
0060 double RK4PreciseStep::distance(const CartesianState& a, const CartesianState& b) const {
0061   return (a.position() - b.position()).mag() + (a.momentum() - b.momentum()).mag() / b.momentum().mag();
0062 }
0063
0064 bool RK4PreciseStep::verbose() const { return true; }