blob: d1f9d556deb05f8d8ac75b940605b3a49ca22d40 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
|
//============================================================================
//
// Ordinary differential equation solver using the Runge-Kutta method.
// $Id$
// Copyright (C) 2004 Georg Drenkhahn
//
// This file is free software; you can redistribute it and/or modify it under
// the terms of the GNU General Public License version 2 as published by the
// Free Software Foundation.
//
//============================================================================
#ifndef RKODESOLVER_H
#define RKODESOLVER_H
// STL headers
#include <valarray>
/** @brief Solver class to integrate a first-order ordinary differential
* equation (ODE) by means of a 6. order Runge-Kutta method.
*
* The ODE system must be given as the derivative
* dy/dx = f(x,y)
* with x in R and y in R^n.
*
* Within this class the function f() is a pure virtual function, which must be
* reimplemented in a derived class.
*
* No other special data type for vectors or matrices are needed besides the STL
* class std::valarray. */
template<typename T>
class RkOdeSolver
{
public:
/** @brief Constructor
* @param x Initial integration parameter
* @param y Initial function values of function to integrate
* @param dx Initial guess for step size. Will be automatically adjusted to
* guarantee required precision.
* @param eps Relative precision
*
* Initialises the solver with start conditions. */
RkOdeSolver(const T& x=0.0,
const std::valarray<T>& y=std::valarray<T>(0),
const T& dx=0,
const T& eps=1e-6);
/** @brief Destructor */
virtual ~RkOdeSolver(void);
/** @brief Integrates the ordinary differential equation from the current x
* value to x+@a dx.
* @param dx x-interval size to integrate over starting from x. dx may be
* negative.
*
* The integration is performed by calling rkStepCheck() repeatedly until the
* desired x value is reached. */
void integrate(const T& dx);
// Accessors
// get/set x value
/** @brief Get current x value.
* @return Reference of x value. */
const T& X(void) const;
/** @brief Set current x value.
* @param a The value to be set. */
void X(const T& a);
// get/set y value
/** @brief Get current y value.
* @return Reference of y vector. */
const std::valarray<T>& Y(void) const;
/** @brief Set current y values.
* @param a The vector to be set. */
void Y(const std::valarray<T>& a);
/** @brief Get current dy/dx value.
* @return Reference of dy/dx vector. */
const std::valarray<T>& dYdX(void) const;
// get/set dx value
/** @brief Get current estimated step size dX.
* @return Reference of dX value. */
const T& dX(void) const;
/** @brief Set estimated step size dX.
* @param a The value to be set. */
void dX(const T& a);
// get/set eps value
/** @brief Get current presision.
* @return Reference of precision value. */
const T& Eps(void) const;
/** @brief Set estimated presision.
* @param a The value to be set. */
void Eps(const T& a);
protected:
// purely virtual function which is integrated
/** @brief ODE function
* @param x Integration value
* @param y Function value
* @return Derivation
*
* This purely virtual function returns the value of dy/dx for the given
* parameter values of x and y. */
virtual std::valarray<T>
f(const T& x, const std::valarray<T>& y) const = 0;
private:
/** @brief Perform one integration step with a tolerable relative error given
* by ::mErr.
* @param dx Maximal step size, may be positive or negative depending on
* integration direction.
* @return Flag indicating if made absolute integration step was equal |@a dx
* | (true) less than |@a dx | (false).
*
* A new estimate for the maximum next step size is saved to ::mStep. The
* new values for x, y and f are saved in ::mX, ::mY and ::mDydx. */
bool rkStepCheck(const T& dx);
/** @brief Perform one Runge-Kutta integration step forward with step size
* ::mStep
* @param dx Step size relative to current x value.
* @param yerr Reference to vector in which the estimated error made in y is
* returned.
* @return The y value after the step at x+@a dx.
*
* Stored current x,y values are not adjusted. */
std::valarray<T> rkStep(const T& dx, std::valarray<T>& yerr) const;
/** current x value */
T mX;
/** current y value */
std::valarray<T> mY;
/** current value of dy/dx */
std::valarray<T> mDydx;
/** allowed relative error */
T mEps;
/** estimated step size for next Runge-Kutta step */
T mStep;
};
// inline accessors
template<typename T>
inline const T&
RkOdeSolver<T>::X(void) const
{
return mX;
}
template<typename T>
inline void
RkOdeSolver<T>::X(const T &a)
{
mX = a;
}
template<typename T>
inline const std::valarray<T>&
RkOdeSolver<T>::Y(void) const
{
return mY;
}
template<typename T>
inline const std::valarray<T>&
RkOdeSolver<T>::dYdX(void) const
{
return mDydx;
}
template<typename T>
inline const T&
RkOdeSolver<T>::dX(void) const
{
return mStep;
}
template<typename T>
inline const T&
RkOdeSolver<T>::Eps(void) const
{
return mEps;
}
#endif
|
