## ExponentialUtilities.jl

Utility functions for exponential integrators
Author SciML
Popularity
16 Stars
Updated Last
11 Months Ago
Started In
August 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

### Required Packages

No packages found.

### Used By Packages

View all packages