solprop.backends.exponential#

The phi functions and exponential Adams weights, in NumPy (SEP 16 sections 4.2-4.3).

Exponential Adams integrates the polynomial through untransported source samples exactly against the linear flow. For a constant operator D, a step h and the Lagrange basis ell_j(theta) = sum_k c_jk theta**k on the normalised nodes,

w_j(x) = int_0^1 exp((1 - theta) x) ell_j(theta) dtheta
= sum_k c_jk k! phi_{k+1}(x),

with x = h D per element and phi_r(x) = sum_m x**m / (m + r)!. At a fraction xi of the step the same integral to xi is sum_k c_jk k! xi**(k+1) phi_{k+1}(xi x). At x = 0 these are the ordinary Adams weights (solprop.stepper.multistep.integrated_weights()).

The phi functions are evaluated by two algorithms, each where it is stable:

near zero (|x| <= 8): scale to y = x / 2**s with |y| <= 1/2, sum the Taylor series of the highest function to a fixed number of terms, recur down with phi_k(y) = 1/k! + y phi_{k+1}(y), which contracts an error by |y| and never divides by y, and undo the scaling by the doubling identity phi_k(2y) = [exp(y) phi_k(y) + sum_{j=1}^k phi_j(y) / (k - j)!] / 2**k, obtained by splitting the defining integral at one half [1] [2];
away from zero (|x| > 8): recur up from phi_0 = exp(x), phi_k(x) = (phi_{k-1}(x) - 1/(k-1)!) / x, whose division by |x| > 8 contracts the error of every earlier function.

Both are checked against augmented matrix exponentials and independent quadrature in tests/backends/test_exponential_weights.py, on both sides of the join. This module is the reference the compiled backends’ moment kernels are checked against.

References#

[1]

S. M. Cox and P. C. Matthews, “Exponential time differencing for stiff systems”, J. Comput. Phys. 176, 430 (2002). doi:10.1006/jcph.2002.6995

[2]

M. Hochbruck and A. Ostermann, “Exponential integrators”, Acta Numerica 19, 209 (2010), section 2. doi:10.1017/S0962492910000048

Module Attributes

`MAX_PHI_ORDER`	Highest phi function evaluated: the corrector's `MAX_ADAMS_ORDER + 1` nodes need `phi_1` to `phi_{MAX_ADAMS_ORDER + 1}`, and one more serves a neighbouring order.
`PHI_JOIN`	`\|x\|` above which the phi functions are recurred up from `exp(x)`.
`COUPLED_SECTIONS`	The coupled flow's section methods the exponential formulation steps (SEP 16 Q3).
`EXPONENTIAL_REPRESENTATIONS`	The linear-operator representations the exponential formulation steps (SEP 16).

Functions

`exponential_weights`(nodes, x[, upper, helpers])	Return the exponential Adams weights of distinct nodes at each argument.
`phi_functions`(x, order)	Return `phi_0(x)` to `phi_order(x)` for each element of x.
`refuse_formulation`(problem, backend)	Refuse the exponential formulation where it is not yet delivered (SEP 16).

solprop.backends.exponential.COUPLED_SECTIONS = frozenset({'constant_schur', 'frozen_residual'})#: The coupled flow’s section methods the exponential formulation steps (SEP 16 Q3).

solprop.backends.exponential.EXPONENTIAL_REPRESENTATIONS = ('constant', 'separable', 'expansion', 'table')#

The linear-operator representations the exponential formulation steps (SEP 16).

A constant operator’s reference flow over a step is its exact flow and its residual is the source alone, R_n = F. A varying one (expansion, table) takes the compensated reference of SEP 16 section 4.2 (phase 4): the chord over the step and the transported compensation (D(z_j) - B_n) u_j, with D the derivative of the representation itself.

solprop.backends.exponential.MAX_PHI_ORDER = 8#: Highest phi function evaluated: the corrector’s MAX_ADAMS_ORDER + 1 nodes need phi_1 to phi_{MAX_ADAMS_ORDER + 1}, and one more serves a neighbouring order.

solprop.backends.exponential.PHI_JOIN = 3.0#

|x| above which the phi functions are recurred up from exp(x).

Below it the upward recurrence loses digits to the subtraction of 1/(k-1)!: at order eight, 291 ulp of the natural scale at |x| = 2 and 32 at 3, against augmented exponentials. Above it the scaled series would need more than three doublings. Three keeps both sides within the witnesses’ 256 ulp with a margin and bounds the doublings, the most expensive part of the near branch (SEP 16 phase 3; eight, the SEP 15 study’s join, cost up to four).

solprop.backends.exponential.exponential_weights(nodes, x, upper=1.0, *, helpers=Helpers(name='numba'))[source]#

Return the exponential Adams weights of distinct nodes at each argument.

w_j = int_0^upper exp((upper - theta) x) ell_j(theta) dtheta, so that exp(upper x) u_n + h sum_j w_j F_j advances u' = D u + F by upper h when F is the polynomial interpolating the samples F_j at the nodes and x = h D.

Parameters:

nodesnumpy.ndarray or list of float: Distinct normalised positions (z_j - z_n) / h.
xnumpy.ndarray or complex: h D per element, dimensionless.
upperfloat, optional: The fraction of the step integrated, in [0, 1].
helperssolprop._helpers.Helpers, optional: The implementation of the Lagrange basis (SEP 29 section 6.6); the shared Numba one by default.

Returns:

numpy.ndarray: One weight per node along a new leading axis, each of x’s shape.

Raises:

StepperError: If the nodes are not distinct and finite, upper lies outside [0, 1] or an argument overflows.

Examples

At x = 0 the weights are the ordinary Adams-Moulton ones:

>>> import numpy as np
>>> from solprop.backends.exponential import exponential_weights
>>> (24 * exponential_weights([1, 0, -1, -2], 0.0)).real.round(12).tolist()
[9.0, 19.0, -5.0, 1.0]

solprop.backends.exponential.phi_functions(x, order)[source]#

Return phi_0(x) to phi_order(x) for each element of x.

phi_0(x) = exp(x) and phi_r(x) = int_0^1 exp((1 - theta) x) theta**(r-1) / (r-1)! dtheta.

Parameters:

xnumpy.ndarray or complex: The arguments, dimensionless (a step times the linear operator). Their real parts must not exceed 700, where exp overflows.
orderint: The highest function, in [1, MAX_PHI_ORDER].

Returns:

numpy.ndarray: phi_r(x) along a new leading axis of length order + 1.

Raises:

StepperError: If an argument is not finite or overflows, or the order is out of range.

Examples

>>> import numpy as np
>>> from solprop.backends.exponential import phi_functions
>>> values = phi_functions(np.array([0.0, 1.0]), 2)
>>> np.round(values.real, 12).tolist()
[[1.0, 2.718281828459], [1.0, 1.718281828459], [0.5, 0.718281828459]]

solprop.backends.exponential.refuse_formulation(problem, backend)[source]#

Refuse the exponential formulation where it is not yet delivered (SEP 16).

The exponential formulation steps constant, separable and varying diagonal operators (phases 2–4) and the constant-Schur and frozen-residual sections of a coupled flow (phase 5); a CF4 section is excluded (Q3) and refused by name before any source is evaluated.

Parameters:

problemsolprop.translate.problem.KernelProblem: The problem.
backendstr: The backend’s name, for the message.

Raises:

BackendError: Naming what is unsupported.