|
ampsci
c++ program for high-precision atomic structure calculations of single-valence systems
|
Solves a 2D system of ODEs using a K-step Adams-Moulton method. More...
#include <AdamsMoulton.hpp>
Public Member Functions | |
| ODESolver2D (Y dt, const DerivativeMatrix< T, Y > *D) | |
| Construct the solver. dt is the (constant) step size. More... | |
| constexpr std::size_t | K_steps () const |
| Returns the AM order (number of steps), K. | |
| Y | last_f () |
| Returns most recent f value. Can also access f array directly. | |
| Y | last_g () |
| Returns most recent g value. Can also access g array directly. | |
| T | last_t () |
| Returns most recent t value; last_f() := f(last_t()) | |
| Y | dt () |
| Returns the step size. | |
| Y | dfdt (Y ft, Y gt, T tt) const |
| Returns derivative, df/dt(t), given f(t),g(t),t. | |
| Y | dgdt (Y ft, Y gt, T tt) const |
| Returns derivative, dg/dt(t), given f(t),g(t),t. | |
| void | drive (T t_next) |
| Drives the DE system to next value, F(t), assuming system has already been solved for the K previous values {t-K*dt, ..., t-dt}. More... | |
| void | drive () |
| Overload of drive(T t) for 'default' case, where next t is defined as last_t + dt (for arithmetic/complex types), or last_t++/last_t– for integral types (grid index). | |
| void | solve_initial_K (T t0, Y f0, Y g0) |
| Automatically sets the first K values for F (and dF), given a single initial value for F, f0=f(t0), fg=g(t0). More... | |
Public Attributes | |
| std::array< Y, K > | f {} |
| Arrays to store the previous K values of f and g. More... | |
| std::array< Y, K > | g {} |
| std::array< Y, K > | df {} |
| Arrays to store the previous K values of derivatives, df and dg. | |
| std::array< Y, K > | dg {} |
| std::array< T, K > | t {} |
| Array to store the previous K values of t: f.at(i) = f(t.at(i)) | |
| Y | S_scale {1.0} |
Solves a 2D system of ODEs using a K-step Adams-Moulton method.
The system of ODEs is defined such that:
\[ \frac{dF(t)}{dt} = D(t) * F(t) + S(t) \]
Where F is a 2D set of functions:
\[ F(t) = \begin{pmatrix} f(t)\\ g(t) \end{pmatrix}, \]
D is the 2x2 "derivative matrix":
\[ D(t) = \begin{pmatrix} a(t) & b(t)\\ c(t) & d(t) \end{pmatrix}, \]
and S(t) is the (optional) 2D inhomogenous term:
\[ S(t) = \begin{pmatrix} s_f(t)\\ s_g(t) \end{pmatrix}. \]
See struct DerivativeMatrix - which is a pure virtual struct that must be implmented by the user in order to define the ODE.
The step-size, dt, must remain constant (since it must remain consistant between the K+1 and previous K points). It may be positive or negative, however (or complex). It's perfectly possible to have a non-uniform grid - this just introduces a Jacobian into the Derivative matrix; dt must still be constant.
The template parameter, T, is the type of the argument of the Derivative Matrix (i.e., type of t). This is often double or complex<double>, but may also be an index type (e.g., std::size_t) if the derivative matrix is only known numerically at certain grid points/stored in an array.
The template parameter, Y, is the type of the function value F(t), and the type of dt, and the return value of the Derivative Matrix. This is often double, but may also be another floating-point type, or std::complex.
The first K points of the function F, and derivative dF/dt, must be known. You may directly access the f,g (function) and df,dg (derivative) arrays, to set these points.
Alternatively, you may use the provided function
which automatically sets the first K values for F (and dF), given a single initial value for F, f0=f(t0), fg=g(t0), by using successive N-step AM methods, for N={1,2,...,K-1}.
For now, just a 2x2 system. In theory, simple to generalise to an N*N system, though requires a matrix inversion.
Example: Bessel's differential equation
\[ x^2 \frac{d^2y}{dx^2} + x \frac{dy}{dx} + (x^2-n^2)y = 0 \]
With y(0) = 1.0, the solutions are the Bessel functions, y(x) = J_n(x)
This can be re-written as a pair of coupled first-order ODEs:
\[ \begin{aligned} \frac{dy}{dx} &= p \\ \frac{dp}{dx} &= \left[\left(\frac{n}{x}\right)^2 - 1\right]y - \frac{1}{x}p \end{aligned} \]
Putting this into the notation/form required for the solver we have:
\[ F(x) = \begin{pmatrix} y(x)\\ p(x) \end{pmatrix} \]
with the "derivative matrix":
\[ D(x) = \begin{pmatrix} 0 & 1 \\ \left(\frac{n}{x}\right)^2 - 1 & \frac{-1}{x} \end{pmatrix} \]
i.e.,
for n=1:
Or, more generally (for example):
Minimal example: – see full examples included elsewhere
|
inline |
Construct the solver. dt is the (constant) step size.
| dt | – (constant) step size |
| D | – Pointer to derivative matrix |
The step-size, dt, may be positive (to drive forwards) or negative (to drive backwards); it may also be complex.
Note: a pointer to the DerivativeMatrix, tD, is stored. This may not be null, and must outlive the ODESolver2D.
|
inline |
Drives the DE system to next value, F(t), assuming system has already been solved for the K previous values {t-K*dt, ..., t-dt}.
Note: t must align with previous values: i.e., t_next = last_t + dt
We re-send t to the function to avoid large build-up of small errors.
There's an overload that avoids this, but care should be taken.
For example, the 10,000th point along t grid may not exactly line up with t0
|
inline |
Automatically sets the first K values for F (and dF), given a single initial value for F, f0=f(t0), fg=g(t0).
It does this by using successive N-step AM methods, for N={1,2,...,K-1}.
i.e.,
F[0] must be given
F[1] determined using F[0] and a 1-step AM method
F[2] determined using F[0],F[1] and a 2-step AM method
...
F[K-1] determined using F[0],F[1],...,F[K-2] and a (K-1)-step AM method
| std::array<Y, K> AdamsMoulton::ODESolver2D< K, T, Y >::f {} |
Arrays to store the previous K values of f and g.
Note: Stored in the same order, regardless of sign of dt (i.e. whether driving forwards or backwards). So f[0] is the 'oldest' value, f[K-1] is newest.
1.9.1