src/FDM/JSBSim/math/FGRungeKutta.h

   1 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   2
   3  Header:       FGRungeKutta.h
   4  Author:       Thomas Kreitler
   5  Date started: 04/9/2010
   6
   7  ------------- Copyright (C)  -------------
   8
   9  This program is free software; you can redistribute it and/or modify it under
  10  the terms of the GNU Lesser General Public License as published by the Free Software
  11  Foundation; either version 2 of the License, or (at your option) any later
  12  version.
  13
  14  This program is distributed in the hope that it will be useful, but WITHOUT
  15  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  16  FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
  17  details.
  18
  19  You should have received a copy of the GNU Lesser General Public License along with
  20  this program; if not, write to the Free Software Foundation, Inc., 59 Temple
  21  Place - Suite 330, Boston, MA  02111-1307, USA.
  22
  23  Further information about the GNU Lesser General Public License can also be found on
  24  the world wide web at http://www.gnu.org.
  25
  26 HISTORY
  27 --------------------------------------------------------------------------------
  28
  29
  30 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  31 SENTRY
  32 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  33
  34 #ifndef FGRUNGEKUTTA_H
  35 #define FGRUNGEKUTTA_H
  36
  37 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  38   INCLUDES
  39   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  40
  41 // #include "FGJSBBase.h" // later
  42
  43 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  44   DEFINITIONS
  45   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  46
  47 #define ID_RUNGEKUTTA "$Id: FGRungeKutta.h,v 1.1 2010/06/02 04:05:13 jberndt Exp $"
  48
  49 namespace JSBSim {
  50
  51 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  52 CLASS DOCUMENTATION
  53 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  54
  55
  56 /**
  57    Minimalistic implementation of some Runge-Kutta methods. Runge-Kutta methods
  58    are a standard for solving ordinary differential equation (ODE) initial
  59    value problems. The code follows closely  the description given on
  60    Wikipedia, see http://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods.
  61
  62    For more powerfull routines see GNU Scientific Library (GSL)
  63    or GNU Plotutils 'ode'.
  64 */
  65
  66
  67 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  68 DECLARATION: FGRungeKuttaProblem
  69 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  70
  71 /**
  72    Abstract base for the function to solve.
  73 */
  74 class FGRungeKuttaProblem {
  75   public:
  76     virtual double pFunc(double x, double y) = 0;
  77 };
  78
  79 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  80 DECLARATION: FGRungeKutta
  81 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  82 /**
  83    Abstract base.
  84 */
  85
  86 class FGRungeKutta {
  87
  88   public:
  89
  90     enum eStates { eNoError=0, eMathError=1, eFaultyInit=2, eEvolve=4, eUnknown=8} ;
  91
  92     int init(double x_start, double x_end, int intervals = 4);
  93
  94     double evolve(double y_0, FGRungeKuttaProblem *pf);
  95
  96     double getXEnd()      { return x_end; }
  97     double getError()     { return err; }
  98
  99     int  getStatus()      { return status; }
 100     int  getIterations()  { return iterations; }
 101     void clearStatus()    { status = eNoError; }
 102     void setTrace(bool t) { trace_values = t; }
 103
 104   protected:
 105     // avoid accidents
 106     FGRungeKutta():  status(eNoError), trace_values(false), iterations(0) {};
 107     virtual ~FGRungeKutta();
 108
 109     FGRungeKuttaProblem *pfo;
 110
 111     double h;
 112     double h05;  // h*0.5, halfwidth
 113     double err;
 114
 115   private:
 116
 117     virtual double approximate(double x, double y) = 0;
 118
 119     bool sane_val(double x);
 120
 121     static const double RealLimit;
 122
 123     double x0, x1;
 124     double safer_x1;
 125     double x_end;
 126
 127     int status;
 128     bool trace_values;
 129     int iterations;
 130
 131 };
 132
 133
 134 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 135 DECLARATION: FGRK4
 136 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 137 /**
 138   Classical RK4.
 139 */
 140
 141 class FGRK4 : public FGRungeKutta {
 142     virtual ~FGRK4();
 143   private:
 144     double approximate(double x, double y);
 145 };
 146
 147
 148 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 149 DECLARATION: FGRKFehlberg
 150 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 151
 152 /**
 153    Runge-Kutta-Fehlberg method.
 154    This is a semi adaptive implementation of rkf - the interval only
 155    shrinks. As a result interval calculations remain trivial, but
 156    sometimes too many calculations are performed.
 157    Rationale: this code is not meant to be a universal pain-reliever
 158    for ode's. Rather it provides some safety if the number of
 159    intervals is set too low, or the problem function behaves a bit
 160    nasty in rare conditions.
 161 */
 162
 163
 164 class FGRKFehlberg : public FGRungeKutta {
 165
 166   public:
 167     FGRKFehlberg() : shrink_avail(4), epsilon(1e-12) { };
 168     virtual ~FGRKFehlberg();
 169     double getEpsilon()          { return epsilon; }
 170     int    getShrinkAvail()      { return shrink_avail; }
 171     void   setEpsilon(double e)  { epsilon = e; }
 172     void   setShrinkAvail(int s) { shrink_avail = s; }
 173
 174   private:
 175
 176     double approximate(double x, double y);
 177
 178     int    shrink_avail;
 179     double epsilon;
 180
 181     static const double A2[], A3[], A4[], A5[], A6[];
 182     static const double B[],  Bs[], C[];
 183
 184 };
 185
 186
 187 } // namespace JSBSim
 188
 189 #endif