Adaptive stiff ODE integrator: a four-stage Rosenbrock method of order four with an embedded order-three estimate for step-size control. More...

#include <cmath>
#include <cstddef>
#include <functional>
#include <limits>
#include <string>
#include <vector>
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/lu.h"
#include "line/util/matrix.h"

Include dependency graph for ode.h:

Classes
struct	line::OdeOptions< T >
	Integration controls. More...
struct	line::OdeSolution< T >
	Result of an integration. More...

Namespaces
namespace	line

Functions
template<class T, class F>
Matrix< T >	line::ode_numeric_jacobian (const F &f, const T &t, const std::vector< T > &y, const std::vector< T > &fy)
	Numeric Jacobian by central differences.
template<class T, class F, class J>
OdeSolution< T >	line::ode_rosenbrock4 (const F &f, const J &jac, const T &t0, const T &t1, const std::vector< T > &y0, const OdeOptions< T > &opt)
	Integrate y' = f(t,y) from t0 to t1 with an analytic Jacobian.
template<class T, class F>
OdeSolution< T >	line::ode_rosenbrock4 (const F &f, const T &t0, const T &t1, const std::vector< T > &y0, const OdeOptions< T > &opt)
	Integrate y' = f(t,y) with a numeric Jacobian by central differences.
template<class T, class F>
std::vector< T >	line::ode_rosenbrock4_endpoint (const F &f, const T &t0, const T &t1, const std::vector< T > &y0)
	Integrate with the default options and return only the state at t1.

Detailed Description

Adaptive stiff ODE integrator: a four-stage Rosenbrock method of order four with an embedded order-three estimate for step-size control.

WHY IT IS HERE. Several LINE algorithms are defined by an initial value problem that MATLAB hands to ode15s or ode23s: the refined mean-field cache approximation (cache_miss_rmf.m), the RANDOM(m) multi-list mean field (cache_rrm_meanfield.m), the fluid solvers. Those right-hand sides are stiff – the mean-field drift of a cache mixes per-item request rates that differ by orders of magnitude, and the relaxation to the fixed point is integrated over a horizon of 1e4 – so an explicit method is not merely slower, it is unusable: its step is capped by the fastest time constant for the whole integration even after every fast mode has died. This header is the port's own integrator; nothing is taken from an external solver library.

THE METHOD. A Rosenbrock (linearly implicit Runge-Kutta) method replaces the nonlinear stage equations of an implicit method by linear ones built on the Jacobian, so each step costs one Jacobian, one LU factorization and s back-substitutions and no Newton iteration ever fails to converge. With J = df/dy and f_t = df/dt evaluated once per step at (t,y), the stages are

(I - h gamma J) k_i = h f(t + alpha_i h, y + sum_{j<i} a_ij k_j)

h J sum_{j<i} gamma_ij k_j
h^2 gamma_i f_t, y_{n+1} = y_n + sum_i b_i k_i, yhat_{n+1} = y_n + sum_i bhat_i k_i,

with alpha_i = sum_{j<i} a_ij and gamma_i = gamma + sum_{j<i} gamma_ij. The f_t term is exactly what the augmented system (y,t)' = (f,1) produces for the y-block, so the method is invariant under autonomization and the order conditions are the autonomous ones. The same matrix I - h gamma J serves every stage, which is the whole point of the constant diagonal gamma.

THE COEFFICIENTS ARE DERIVED HERE, NOT COPIED. They were obtained by imposing the order conditions directly rather than by quoting a published table: for random polynomial vector fields the one-step numerical solution was expanded as a truncated power series in h and matched, coefficient by coefficient, against the exact Taylor series of the solution, and the parameters solved so that the h^1 through h^4 coefficients agree. That is the definition of order four, with no intermediate rooted-tree bookkeeping to get wrong. Two further requirements were imposed at the same time: b^T B^-1 1 = 1 with B = A + Gamma the lower triangular matrix with constant diagonal gamma, which makes R(-inf) = 0 and the method L-stable – an A-stable but not L-stable method leaves the fastest modes ringing at |R| = 1 instead of damping them, precisely the failure a stiff integrator exists to avoid – and gamma pinned to a value for which |R(z)| <= 1 holds along the whole negative real axis.

The embedded estimate is SECOND order and uses the first three stages. With four stages and four weights the order-three conditions have a unique solution, which is the order-four weight vector itself, so no order-three estimate with these stages exists; the weights below instead reproduce the solution through h^2 exactly and use their remaining degree of freedom to make the h^3 mismatch as small as possible, which makes y - yhat a sharp estimate of the O(h^3) term. The step controller therefore uses the exponent 1/(2+1) = 1/3.

What the test suite checks about the coefficients is what can be checked exactly: alpha is the row sum of a, gamma_i is gamma plus the row sum of Gamma, the linear order conditions b^T B^(k-1) 1 = 1/k! for k = 1..4 hold (on a linear problem the method IS the implicit Runge-Kutta with matrix B, so these are necessary and sufficient there), the embedded weights are consistent and differ from b, |R| <= 1 on the negative axis with R(-inf) = 0, and the observed convergence rate on a nonlinear non-autonomous problem is four. A mistyped digit fails at least one of those.

ARITHMETIC. Gated on num_traits<T>::has_transcendental. The coefficients are irrational, the step size is chosen by a tolerance comparison and the result is an approximation controlled by rtol/atol no matter how the arithmetic is carried out, so an exact-rational instantiation would be a fiction. double and Real<D> both instantiate; at Real<D> the coefficients are parsed from their decimal strings, so the method keeps its order at any precision, and the numeric Jacobian's difference increment is scaled by the precision of T.

DETERMINISM. No global state, no static mutable data, no clock, no random numbers. The same inputs return the same trajectory bit for bit, including the sequence of accepted and rejected steps. Every tolerance and bound is supplied by the caller through OdeOptions.

Definition in file ode.h.

Classes

Namespaces

Functions

Detailed Description