1 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4 Author: Thomas Kreitler
5 Date started: 04/9/2010
7 ------------- Copyright (C) -------------
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
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
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.
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.
27 --------------------------------------------------------------------------------
30 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
32 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
34 #ifndef FGRUNGEKUTTA_H
35 #define FGRUNGEKUTTA_H
37 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
39 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
41 // #include "FGJSBBase.h" // later
43 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
45 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
47 #define ID_RUNGEKUTTA "$Id: FGRungeKutta.h,v 1.1 2010/06/02 04:05:13 jberndt Exp $"
51 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
53 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
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.
62 For more powerfull routines see GNU Scientific Library (GSL)
63 or GNU Plotutils 'ode'.
67 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
68 DECLARATION: FGRungeKuttaProblem
69 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
72 Abstract base for the function to solve.
74 class FGRungeKuttaProblem {
76 virtual double pFunc(double x, double y) = 0;
79 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
80 DECLARATION: FGRungeKutta
81 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
90 enum eStates { eNoError=0, eMathError=1, eFaultyInit=2, eEvolve=4, eUnknown=8} ;
92 int init(double x_start, double x_end, int intervals = 4);
94 double evolve(double y_0, FGRungeKuttaProblem *pf);
96 double getXEnd() { return x_end; }
97 double getError() { return err; }
99 int getStatus() { return status; }
100 int getIterations() { return iterations; }
101 void clearStatus() { status = eNoError; }
102 void setTrace(bool t) { trace_values = t; }
106 FGRungeKutta(): status(eNoError), trace_values(false), iterations(0) {};
107 virtual ~FGRungeKutta();
109 FGRungeKuttaProblem *pfo;
112 double h05; // h*0.5, halfwidth
117 virtual double approximate(double x, double y) = 0;
119 bool sane_val(double x);
121 static const double RealLimit;
134 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
136 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
141 class FGRK4 : public FGRungeKutta {
144 double approximate(double x, double y);
148 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
149 DECLARATION: FGRKFehlberg
150 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
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.
164 class FGRKFehlberg : public FGRungeKutta {
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; }
176 double approximate(double x, double y);
181 static const double A2[], A3[], A4[], A5[], A6[];
182 static const double B[], Bs[], C[];
187 } // namespace JSBSim