ExponentialUtilities is a package of utility functions used by the exponential integrators in OrdinaryDiffEq. It is a lightweight pure Julia package with no external dependencies, so it can also be used independently.
The main functionality of ExponentialUtilities is the computation of matrix-phi-vector products. The phi functions are defined as
ϕ_0(z) = exp(z) ϕ_(k+1)(z) = (ϕ_k(z) - 1) / z
In exponential algorithms, products in the form of
ϕ_m(tA)b is frequently encountered. Instead of computing the matrix function first and then computing the matrix-vector product, the common alternative is to construct a Krylov subspace
K_m(A,b) and then approximate the matrix-phi-vector product.
expv(t,A,b;kwargs) -> exp(tA)b phiv(t,A,b,k;kwargs) -> [ϕ_0(tA)b ϕ_1(tA)b ... ϕ_k(tA)b][, errest]
ϕ_m(tA)b products up to order
k is returned as a matrix. This is because it's more economical to produce all the results at once than individually. A second output is returned if
kwargs. The error estimate is given for the second-to-last product, using the last product as an estimator. If
ϕ_(k-1) are updated using the last Arnoldi vector. The correction algorithm is described in .
You can adjust how the Krylov subspace is constructed by setting various keyword arguments. See the Arnoldi iteration section for more details.
phiv, the timestepping methods divide
t into smaller time steps and compute the product step-by-step. By doing this in smaller chunks, the methods allow for finer error control as well as adaptation. The timestepping algorithm is described in , which is based upon the numerical package Expokit .
exp_timestep(ts,A,b;kwargs) -> U
Evaluates the matrix exponentiation-vector product
u = exp(tA)b using time stepping.
phiv_timestep(ts,A,[b_0 b_1 ... b_p];kwargs) -> U
Evaluates the linear combination of phi-vector products
u = ϕ_0(tA)b_0 + tϕ_1(tA)b_1 + ... + t^pϕ_p(tA)b_p using time stepping.
In both cases,
ts is an array of time snapshots for u, with
U[:,j] ≈ u(ts[j]).
ts can also be just one value, in which case only the end result is returned and
U is a vector.
Apart from keyword arguments that affect the computation of Krylov subspaces (see the Arnoldi iteration section), you can also adjust the timestepping behavior using the arguments. By setting
adaptive=true, the time step and Krylov subsapce size adaptation scheme of Niesen & Wright is used and the relative tolerance of which can be set using the keyword parameter
gamma parameter of the adaptation scheme can also be adjusted. The
tau parameter adjusts the size of the timestep (and for
adaptive=true, the initial timestep). By default, it is calculated using a heuristic formula by Niesen & Wright.
You can use any object as the "matrix"
A as long as it implements the following linear operator interface:
LinearAlgebra.mul!(y, A, x)(for computing
y = A * xin place).
LinearAlgebra.opnorm(A, p=Inf). If this is not implemented or the default implementation can be slow, you can manually pass in the operator norm (a rough estimate is fine) using the keyword argument
LinearAlgebra.ishermitian(A). If this is not implemented or the default implementation can be slow, you can manually pass in the value using the keyword argument
phi(z,k[;cache]) -> [ϕ_0(z),ϕ_1(z),...,ϕ_k(z)]
Compute ϕ function directly.
z can both be a scalar or a
AbstractMatrix (note that unlike the previous functions, you need to use a concrete matrix). This is used by the caching versions of the ExpRK integrators to precompute the operators.
Instead of using the recurrence relation, which is numerically unstable, a formula given by Sidje is used .
arnoldi(A,b[;m,tol,opnorm,iop,cache]) -> Ks
Performs Arnoldi iterations to obtain the Krylov subspace
K_m(A,b). The result is a
KrylovSubspace that can be used by
phiv via the alternative interface
phiv(t,Ks,k;kwargs) -> [ϕ_0(tA)b ϕ_1(tA)b ... ϕ_k(tA)b][, errest]
The reason for having this alternative interface is that we may want to compute
ϕ_m(tA)b for different values of
t. In this case, we can compute
Ks just once (which is expensive) and follow up with several
phiv calls using
Ks (which is not as expensive).
A is hermitian, then the more efficient Lanczos algorithm is used instead. For cases when
A is almost hermitian or when accuracy is not important, the incomplete orthogonalization procedure (IOP) can be used by setting the IOP length
For the other keyword arguments,
m determines the dimension of the Krylov subspace and
tol is the relative tolerance used to determine the "happy-breakdown" condition. You can also set custom operator norm in
opnorm, e.g. efficient norm estimation functions instead of the default
opnorm(A, Inf) needs to be defined.
"In-place" versions for
phiv_timestep are available as
phiv_timestep!. You can refer to the docstrings for more information.
In addition, you may provide pre-allocated caches to the functions to further improve efficiency. In particular, dedicated cache types for
phiv! are available as
PhivCache. Both of them can be dynamically resized, which is crucial for the adaptive
 Niesen, J., & Wright, W. (2009). A Krylov subspace algorithm for evaluating the φ-functions in exponential integrators. arXiv preprint arXiv:0907.4631.
 Sidje, R. B. (1998). Expokit: a software package for computing matrix exponentials. ACM Transactions on Mathematical Software (TOMS), 24(1), 130-15.
 Koskela, A. (2015). Approximating the matrix exponential of an advection-diffusion operator using the incomplete orthogonalization method. In Numerical Mathematics and Advanced Applications-ENUMATH 2013 (pp. 345-353). Springer, Cham.
15 days ago