- Generated on Tue May 12 2026 01:52:49 for ampsci by
1.9.8
|
ampsci
High-precision calculations for one- and two-valence atomic systems
|
Solves a 2x2 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
#include <AdamsMoulton.hpp>
Public Member Functions | |
| ODESolver2D (Y dt, const DerivativeMatrix< T, Y > *D) | |
| Constructs the ODE solver with a given step size and derivative matrix. | |
| 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 ODE system one step to t_next, given the K previous values. | |
| 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) |
| Sets the first K values of F (and dF) given a single initial condition. | |
Public Attributes | |
| std::array< Y, K > | f {} |
| Arrays storing the previous K values of f and g. | |
| 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} |
|
inline |
Constructs the ODE solver with a given step size and derivative matrix.
The step-size dt may be positive (drive forwards), negative (drive backwards), or complex. A raw pointer to D is stored internally and must not be null; it must outlive the ODESolver2D instance.
| dt | Constant step size. |
| D | Pointer to the derivative matrix. Must not be null. |
|
inlineconstexpr |
Returns the AM order (number of steps), K.
|
inline |
Returns most recent f value. Can also access f array directly.
|
inline |
Returns most recent g value. Can also access g array directly.
|
inline |
Returns most recent t value; last_f() := f(last_t())
|
inline |
Returns the step size.
|
inline |
Returns derivative, df/dt(t), given f(t),g(t),t.
|
inline |
Returns derivative, dg/dt(t), given f(t),g(t),t.
|
inline |
Drives the ODE system one step to t_next, given the K previous values.
Assumes the system has already been solved for the K previous values {t-K*dt, ..., t-dt}. The value t_next should satisfy t_next = last_t + dt.
t is passed explicitly to avoid accumulation of floating-point errors: the 10,000th grid point may not equal t0 + 10000*dt exactly, particularly on non-linear grids. The no-argument overload drive() avoids this but should be used with care.
The type T of t_next must match the type expected by DerivativeMatrix; usually T=double for an analytic derivative, or T=std::size_t when the derivative is stored on a discrete grid.
| t_next | The target t value for the new step. |
|
inline |
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).
|
inline |
Sets the first K values of F (and dF) given a single initial condition.
Uses successive N-step AM methods for N = {1, 2, ..., K-1}:
| t0 | Initial value of t. |
| f0 | Initial value f(t0). |
| g0 | Initial value g(t0). |
| std::array<Y, K> AdamsMoulton::ODESolver2D< K, T, Y >::f {} |
Arrays storing the previous K values of f and g.
Stored in chronological order regardless of the sign of dt (i.e. whether driving forwards or backwards). f[0] is the oldest value, f[K-1] is newest.
| std::array<Y, K> AdamsMoulton::ODESolver2D< K, T, Y >::df {} |
Arrays to store the previous K values of derivatives, df and dg.
| std::array<T, K> AdamsMoulton::ODESolver2D< K, T, Y >::t {} |
Array to store the previous K values of t: f.at(i) = f(t.at(i))
1.9.8