Utility functions for exponential integrators
Julia TeX Fetching latest commit…
Cannot retrieve the latest commit at this time.
Type Name Latest commit message Commit time
Failed to load latest commit information. .github/workflows Nov 21, 2019 src Dec 9, 2019 test Dec 9, 2019 .codecov.yml Aug 7, 2018 .gitignore Jun 22, 2019 .travis.yml Oct 9, 2018 CITATION.bib Nov 21, 2019 LICENSE.md Aug 7, 2018 Project.toml Dec 11, 2019 README.md Dec 11, 2019 appveyor.yml Oct 9, 2018

# ExponentialUtilities

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.

## Matrix-phi-vector product

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` and `phiv`

```expv(t,A,b;kwargs) -> exp(tA)b
phiv(t,A,b,k;kwargs) -> [ϕ_0(tA)b ϕ_1(tA)b ... ϕ_k(tA)b][, errest]```

For `phiv`, all `ϕ_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 `errest=true` in `kwargs`. The error estimate is given for the second-to-last product, using the last product as an estimator. If `correct=true`, then `ϕ_0` through `ϕ_(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.

### `expv_timestep` and `phiv_timestep`

Unlike `expv` and `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 `tol`. The `delta` and `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.

### Support for matrix-free operators

You can use any object as the "matrix" `A` as long as it implements the following linear operator interface:

• `Base.eltype(A)`
• `Base.size(A, dim)`
• `LinearAlgebra.mul!(y, A, x)` (for computing `y = A * x` in 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 `opnorm`.
• `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 `ishermitian`.

## Matrix-phi function `phi`

`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 iteration `arnoldi`

`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).

For `arnoldi`, if `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 `iop` in `kwargs`.

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 `LinearAlgebra.opnorm`. Only `opnorm(A, Inf)` needs to be defined.

## `_exp!`

`_exp(A)`

A pure Julia implementation of a non-allocating matrix exponential using the Destructive matrix exponential using algorithm from Higham, 2008. Mostly generic, though the coefficients are geared towards 64-bit floating point calculations, and the use of BLAS requires a `StridedMatrix`.

## `exp_generic`

`exp(x, vk=Val{10}())`

A pure Julia generic implementation of the exponential function using the scaling and squaring method, working on any `x` for which the functions `LinearAlgebra.opnorm`, `+`, `*`, `^`, and `/` (including addition with UniformScaling objects) are defined. Use the argument `vk` to adjust the number of terms used in the Pade approximants at compile time.

"In-place" versions for `phi`, `arnoldi`, `expv`, `phiv`, `expv_timestep` and `phiv_timestep` are available as `phi!`, `arnoldi!`, `expv!`, `phiv!`, `expv_timestep!` and `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 `expv!` and `phiv!` are available as `ExpvCache` and `PhivCache`. Both of them can be dynamically resized, which is crucial for the adaptive `phiv_timestep!` method.

## References

 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.

 HIGHAM, N. J. (2005). "THE SCALING AND SQUARING METHOD FOR THE MATRIX EXPONENTIAL REVISITED." SIAM J. MATRIX ANAL. APPL.Vol. 26, No. 4, pp. 1179–1193

You can’t perform that action at this time.