RungeKutta/src/RK4PreciseSolver.icc

0001 #include "RK4PreciseSolver.h"
0002 #include "RK4OneStepTempl.h"
0003 #include <cmath>
0004
0005 inline double signum( double d) {return d>=0 ? 1.0 : -1.0;}
0006
0007 template <typename T, int N>
0008 typename RKSolver<T,N>::Vector
0009 RK4PreciseSolver<T,N>::operator()( Scalar startPar, const Vector& startState,
0010                    Scalar step, const RKDerivative<T,N>& deriv,
0011                    const RKDistance<T,N>& dist,
0012                    float eps)
0013 {
0014     const double Safety = 0.9;
0015     double remainigStep = step;
0016     double stepSize = step;   // attempt to solve in one step
0017     Scalar currentPar   = startPar;
0018     Vector currentStart = startState;
0019     int nsteps = 0;
0020     std::pair<Vector, Scalar> tryStep;
0021
0022     do {
0023     tryStep = stepWithAccuracy( currentPar, currentStart, deriv, dist, stepSize);
0024     nsteps++;
0025     if (tryStep.second <eps) {
0026         if (abs(remainigStep - stepSize) < eps/2) {
0027         if (verbose()) cout << "Accuracy reached, and full step taken in "
0028                     << nsteps << " steps" << endl;
0029         return tryStep.first; // we are there
0030         }
0031         else {
0032         remainigStep -= stepSize;
0033         currentPar += stepSize;
0034                 // increase step size
0035         double factor =  min( Safety * pow( fabs(eps/tryStep.second),0.2), 4.);
0036         // stepSize = min( stepSize*factor, remainigStep);
0037         double absSize =  min( abs(stepSize*factor), abs(remainigStep));
0038         stepSize = absSize * signum(stepSize);
0039         currentStart = tryStep.first;
0040         if (verbose()) cout << "Accuracy reached, but " << remainigStep
0041              << " remain after " << nsteps << " steps. Step size increased by "
0042              << factor << " to " << stepSize << endl;
0043         }
0044     }
0045     else {
0046         // decrease step size
0047         double factor =  max( Safety * pow( fabs(eps/tryStep.second),0.25), 0.1);
0048         stepSize *= factor;
0049         if (verbose()) cout << "Accuracy not yet reached: delta = " << tryStep.second
0050          << ", step reduced by " << factor << " to " << stepSize << endl;
0051     }
0052     } while (abs(remainigStep) > eps/2);
0053
0054     return tryStep.first;
0055 }
0056 template <typename T, int N>
0057 std::pair< typename RKSolver<T,N>::Vector, T>
0058 RK4PreciseSolver<T,N>::stepWithAccuracy( Scalar startPar, const Vector& startState,
0059                      const RKDerivative<T,N>& deriv,
0060                      const RKDistance<T,N>& dist, float step)
0061 {
0062   const Scalar huge = 1.e5;
0063   const Scalar hugediff = 100.;
0064
0065     RK4OneStepTempl<T,N> solver;
0066     Vector one(       solver(startPar, startState, deriv, step));
0067     if (abs(one[0])>huge || abs(one(1)>huge)) return std::pair<Vector, Scalar>(one,hugediff);
0068
0069     Vector firstHalf( solver(startPar, startState, deriv, step/2));
0070     Vector secondHalf(solver(startPar+step/2, firstHalf, deriv, step/2));
0071     Scalar diff = dist(one, secondHalf, startPar+step);
0072     return std::pair<Vector, Scalar>(secondHalf,diff);
0073 }