 |
|
|
|
Double precision first-order ODE (IVP, i.e. More...
#include <RunKutODESolver.h>
Public Member Functions | |
| RunKutODESolver (double(T::*fce)(double, double), T *pTemplate, float eps) | |
| Constructor - sets ODE right-hand side and required absolut accumulated precision. More... | |
| ~RunKutODESolver () | |
| Destructor. More... | |
| double | RKSolve (double x0, double y0, double x) |
| This method returns a solution y(x) to the so-called ODE initial value problem (IVP), using 4th order Runge-Kutta method. More... | |
| double | RKCSolve (double x0, double y0, double x, double hTry) |
| This method returns a solution y(x) to the so-called ODE initial value problem (IVP), using 4th order Runge-Kutta method with adaptive stepsize control. More... | |
| double | RKFSolve (double x0, double y0, double x, double hTry) |
| This method returns a solution y(x) to the so-called ODE initial value problem (IVP), using 5th(4th) order Runge-Kutta-Fehlberg method with adaptive stepsize control. More... | |
Protected Member Functions | |
| double | RK4Method (double x_i, double w_i, double h, double dydx_i) |
| This method based on 4th order Runge-Kutta algorithm returns an approximation w(x_i+1,h) to the exact solution y(x_i+1) of the so-called initial value problem (IVP) at position x_i+1 = x_i + h and using given "integration" stepsize h, where the IVP is defined as y' = f(x,y), y(x0) = y0 and for equidistant stepsizes h the coordinate x_i can be expressed as x_i = x0 + ih. More... | |
| double | RKFMethod (double x_i, double w_i, double h, double dydx_i, double *wDiffPointer) |
| This method based on 5th(4th) order Runge-Kutta-Fehlberg algorithm returns two approximations: w(x_i+1,h) (O(h6)) and v(x_i+1,h) (O(h5)), that approximate the exact solution y(x_i+1) of given initial value problem (IVP) at position x_i+1 = x_i + h and using an integration stepsize h. More... | |
Protected Attributes | |
| double(T::* | _fce )(double, double) |
| Relative pointer to an integrated function. More... | |
| T * | _pTemplate |
| Pointer to the template class. More... | |
| float | _eps |
| Absolute precision (to hold "global" accumulation of errors constant) More... | |
Double precision first-order ODE (IVP, i.e.
initial-value problem) solver, that utilizes 4th order Runge-Kutta algorithm with or without adaptive control of its stepsize or so-called 5th(4th) order Runge-Kutta-Fehlberg algorithm with adaptive control of its stepsize (The template class represents class whose method corresponds to a right-hand side FCE2(x,y) of the ODE, defined as dy/dx = FCE2(x,y)!). The template class doesn't have *.cpp file!
Definition at line 35 of file RunKutODESolver.h.
|
inline |
Constructor - sets ODE right-hand side and required absolut accumulated precision.
Definition at line 40 of file RunKutODESolver.h.
|
inline |
Destructor.
Definition at line 44 of file RunKutODESolver.h.
|
protected |
This method based on 4th order Runge-Kutta algorithm returns an approximation w(x_i+1,h) to the exact solution y(x_i+1) of the so-called initial value problem (IVP) at position x_i+1 = x_i + h and using given "integration" stepsize h, where the IVP is defined as y' = f(x,y), y(x0) = y0 and for equidistant stepsizes h the coordinate x_i can be expressed as x_i = x0 + ih.
General numerical algorithm for approximate values' calculation:
w_0 = y0 w_i+1 = w_i + h*Phi(x_i,w_i;h;f) x_i+1 = x_i + h, x_0 = x0
Euler method gives Phi(x,y;h;f) = f(x,y)
Calculation using Runge-Kutta algorithm:
k1 = f(x_i, w_i) k2 = f(x_i + h/2, w_i + h/2*k1) k3 = f(x_i + h/2, w_i + h/2*k2) k4 = f(x_i + h, w_i + h*k3)
Phi(x_i,w_i;h;f) = 1/6*[k1 + 2*k2 + 2*k3 + k4]
Parameters: position x_i, approximation w_i, stepsize h and evaluated righ-hand side of the ODE dydx_i (It may be called more than once, thus given as param!)
For more details see J.Stoer, R.Bulirsch: Introduction to Numerical Analysis, chap. 7, 3rd ed., and www.nr.com
Definition at line 465 of file RunKutODESolver.h.
| double sistrip::RunKutODESolver< T >::RKCSolve | ( | double | x0, |
| double | y0, | ||
| double | x, | ||
| double | hTry | ||
| ) |
This method returns a solution y(x) to the so-called ODE initial value problem (IVP), using 4th order Runge-Kutta method with adaptive stepsize control.
The stepsize control algorithm is performed as follows: first, an approximative solution w(x0+h;h) using stepsize h is calculated and then is compared to the solution w(x0+h;h/2), obtained in the same way, but using two consequent stepsizes h/2. If the difference is not similar to the required epsilon, then the ratio of a new proposed stepsize and the old stepsize must be much bigger than 2 (i.e. take for example bigger or equal than 3), else it's approximately equal to 2. Hence, in both cases the stepsize is set as 2 times stepsize H, that corresponds to required precision, but in the first case, one must perform the algorithm once again, while in the latter the approximative solution w_i+1 and a stepsize for new iteration are found. The epsilon corresponds to precision in each step!!! The IVP problem is defined as:
y' = dy/dx = f(x,y) with initial condition f(x0) = y0
Calculation (w(x,h) corresponds to the best approximation of y(x)):
w_0 = y0
for i = 0,1,2 ... (i-th step)
w_i+1 = w_i + h*Phi(x_i,w_i;h;f)
x_i+1 = x_i + h, x_0 = x0
Euler method gives Phi(x,y;h;f) = f(x,y)
Stepsize calculation (H denotes required stepsize, h current stepsize):
h/H = (16/15*(|w_i+1(x_i+h;h) - w_i+1(x_i+h;h/2)|/eps))^0.20
if h/H >> 2 -> h = 2*H -> required precision not achieved -> stepsize recalculation
if h/H ~ 2 -> h = 2*H -> required precision achieved -> use h for another iteration
For more details see J.Stoer, R.Bulirsch: Introduction to Numerical Analysis, chap. 7, 3rd ed., and www.nr.com
Definition at line 267 of file RunKutODESolver.h.
|
protected |
This method based on 5th(4th) order Runge-Kutta-Fehlberg algorithm returns two approximations: w(x_i+1,h) (O(h6)) and v(x_i+1,h) (O(h5)), that approximate the exact solution y(x_i+1) of given initial value problem (IVP) at position x_i+1 = x_i + h and using an integration stepsize h.
These approximative solutions are effectively used, when controlling required precision. The difference is proportional to h*C(x) and hence can be used to estimate an error and find a more suitable stepsize H. Due to a fact, that different schemes can be used to calculate w or v, we've chosen the one of Cash and Karp. The IVP is defined as y' = f(x,y), y(x0) = y0.
Rung-Kutta-Fehlberg algorithm for approximate values' calculation:
w_0 = y0
v_0 = y0
w_i+1 = w_i + h*PhiI(x_i,w_i;h;f)
v_i+1 = v_i + h*PhiII(x_i,w_i;h;f)
x_i+1 = x_i + h, x_0 = x0
PhiI(x_i,w_i;h;f) = Sum(k=0, 5){c_kI * f_k(x_i, w_i; h)}
PhiII(x_i,w_i;h;f) = Sum(k=0, 5){c_kII * f_k(x_i, w_i; h)}
f_k(x_i,w_i;h) = f(x_i + a_k*h, w_i + h * Sum(l=0, k-1){b_kl * f_l)
a_k, b_kl, c_k ...... Cash-Karp parameters
| k | a_k | b_k0 | b_k1 | b_k2 | b_k3 | b_k4 | c_kI | c_kII |
| 0 | 0 | 37/378 | 2825/27648 | |||||
| 1 | 1/5 | 1/5 | 0 | 0 | ||||
| 2 | 3/10 | 3/40 | 9/40 | 250/621 | 18575/48384 | |||
| 3 | 3/5 | 3/10 | -9/10 | 6/5 | 125/594 | 13525/55296 | ||
| 4 | 1 | -11/54 | 5/2 | -70/27 | 35/27 | 0 | 277/14336 | |
| 5 | 7/8 | 1631/55296 | 175/512 | 575/13824 | 44275/110592 | 253/4096 | 512/1771 | 1/4 |
Parameters: position x_i, approximation w_i, stepsize h, evaluated righ-hand side of the ODE dydx_i and estimated error between 5th and 4th order (i.e. the difference between w_i+1 and v_i+1, defined as wDiff).
For more details see J.Stoer, R.Bulirsch: Introduction to Numerical Analysis, chap. 7, 3rd ed., and www.nr.com
Definition at line 494 of file RunKutODESolver.h.
| double sistrip::RunKutODESolver< T >::RKFSolve | ( | double | x0, |
| double | y0, | ||
| double | x, | ||
| double | hTry | ||
| ) |
This method returns a solution y(x) to the so-called ODE initial value problem (IVP), using 5th(4th) order Runge-Kutta-Fehlberg method with adaptive stepsize control.
The stepsize control algorithm is performed as follows: first two approximations to the exact solution w(x_i+1,h) (O(h6)) and v(x_i+1,h) (O(h5)) are calculated and then the difference proportional to h*C(x) is expressed and used for calculation of a new stepsize, i.e. used as an estimate of the error. Due to a fact that the required precision _eps represents an accumulated error, permitted in maximum, C(x) is also dependent on h and thus expressed in more complicated way (Two coefficients: 0.25 and 0.20 are used), see Stepsize calculation. The IVP problem is defined as:
y' = dy/dx = f(x,y) with initial condition f(x0) = y0
Calculation (w(x,h), resp. v(x,h), correspond to the best approximation of y(x)):
w_0 = y0
v_0 = y0
for i = 0,1,2 ... (i-th step)
w_i+1 = w_i + h*PhiI(x_i,w_i;h;f)
v_i+1 = v_i + h*PhiII(x_i,w_i;h;f)
x_i+1 = x_i + h, x_0 = x0
Stepsize calculation (H denotes required stepsize, h current stepsize):
fraction = |eps*h/(x - x0)/(h*PhiI - h*PhiII)|
if fraction <= 1. -> reduce stepsize h = 0.9*h*fraction^0.25
if fraction > 1. -> increase stepsize h = 0.9*h*fraction^0.2
For more details see J.Stoer, R.Bulirsch: Introduction to Numerical Analysis, chap. 7, 3rd ed., and www.nr.com
Definition at line 366 of file RunKutODESolver.h.
| double sistrip::RunKutODESolver< T >::RKSolve | ( | double | x0, |
| double | y0, | ||
| double | x | ||
| ) |
This method returns a solution y(x) to the so-called ODE initial value problem (IVP), using 4th order Runge-Kutta method.
Where IVP is defined as follows:
y' = dy/dx = f(x,y) with initial condition f(x0) = y0
Calculation (w(x,h) corresponds to the best approximation of y(x)):
w_0 = y0
for i = 0,1,2 ... (i-th step)
w_i+1 = w_i + h*Phi(x_i,w_i;h;f)
x_i+1 = x_i + h, x_0 = x0
Euler method gives Phi(x,y;h;f) = f(x,y)
This algorithm doesn't implement adaptive stepsize control, each stepsize is calculated as relative precision * (x-xo). (Epsilon ~ 0.001, i.e 1000 iterations, gives approximately good results.)
For more details see J.Stoer, R.Bulirsch: Introduction to Numerical Analysis, chap. 7, 3rd ed., and www.nr.com
Definition at line 230 of file RunKutODESolver.h.
|
protected |
Absolute precision (to hold "global" accumulation of errors constant)
Definition at line 222 of file RunKutODESolver.h.
|
protected |
Relative pointer to an integrated function.
Definition at line 218 of file RunKutODESolver.h.
|
protected |
Pointer to the template class.
Definition at line 220 of file RunKutODESolver.h.
1.8.5