Branch data Line data Source code
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 : 0 : 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 : 0 : int getStatus() { return status; }
100 : 0 : int getIterations() { return iterations; }
101 : 0 : void clearStatus() { status = eNoError; }
102 : : void setTrace(bool t) { trace_values = t; }
103 : :
104 : : protected:
105 : : // avoid accidents
106 : 0 : 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 : 0 : 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
|
