Expv: Matrix Exponentials Times Vectors
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) / zIn 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.
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 computingy = 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 argumentopnorm.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 argumentishermitian.
Core API
ExponentialUtilities.expv — Functionexpv(t,A,b; kwargs) -> exp(tA)bCompute the matrix-exponential-vector product using Krylov.
A Krylov subspace is constructed using arnoldi and exp! is called on the Hessenberg matrix. Consult arnoldi for the values of the keyword arguments. An alternative algorithm, where an error estimate generated on-the-fly is used to terminate the Krylov iteration, can be employed by setting the kwarg mode=:error_estimate.
expv(t,Ks; cache) -> exp(tA)bCompute the expv product using a pre-constructed Krylov subspace.
ExponentialUtilities.phiv — Functionphiv(t,A,b,k;correct,kwargs) -> [phi_0(tA)b phi_1(tA)b ... phi_k(tA)b][, errest]Compute the matrix-phi-vector products using Krylov. k >= 1.
The phi functions are defined as
\[\varphi_0(z) = \exp(z),\quad \varphi_{k+1}(z) = \frac{\varphi_k(z) - 1}{z}\]
A Krylov subspace is constructed using arnoldi and phiv_dense is called on the Hessenberg matrix. If correct=true, then phi0 through phik-1 are updated using the last Arnoldi vector v_m+1 [1]. If errest=true then an additional error estimate for the second-to-last phi is also returned. For the additional keyword arguments, consult arnoldi.
phiv(t,Ks,k;correct,kwargs) -> [phi0(tA)b phi1(tA)b ... phi_k(tA)b][, errest]
Compute the matrix-phi-vector products using a pre-constructed Krylov subspace.
the φ-functions in exponential integrators. arXiv preprint arXiv:0907.4631. Formula (10).
ExponentialUtilities.expv! — Functionexpv!(w,t,Ks[;cache]) -> wNon-allocating version of expv that uses precomputed Krylov subspace Ks.
expv!(w, t, A, b, Ks, cache)Alternative interface for calculating the action of exp(t*A) on the vector b, storing the result in w. The Krylov iteration is terminated when an error estimate for the matrix exponential in the generated subspace is below the requested tolerance. Ks is a KrylovSubspace and typeof(cache)<:HermitianSubspaceCache, the exact type decides which algorithm is used to compute the subspace exponential.
ExponentialUtilities.phiv! — Functionphiv!(w,t,Ks,k[;cache,correct,errest]) -> w[,errest]Non-allocating version of 'phiv' that uses precomputed Krylov subspace Ks.
Missing docstring for exp_timestep. Check Documenter's build log for details.
ExponentialUtilities.phiv_timestep — Functionphiv_timestep(ts,A,B[;adaptive,tol,kwargs...]) -> UEvaluates the linear combination of phi-vector products using time stepping
\[u = \varphi_0(tA)b_0 + t\varphi_1(tA)b_1 + \cdots + t^p\varphi_p(tA)b_p\]
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.
The time stepping formula of Niesen & Wright is used [1]. If the time step tau is not specified, it is chosen according to (17) of Neisen & Wright. If adaptive==true, the time step and Krylov subsapce size adaptation scheme of Niesen & Wright is used, 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.
Set verbose=true to print out the internal steps (for debugging). For the other keyword arguments, consult arnoldi and phiv, which are used internally.
evaluating the φ-functions in exponential integrators. arXiv preprint arXiv:0907.4631.
Missing docstring for exp_timestep!. Check Documenter's build log for details.
ExponentialUtilities.phiv_timestep! — Functionphiv_timestep!(U,ts,A,B[;kwargs]) -> UNon-allocating version of phiv_timestep.
ExponentialUtilities.phi — Functionphi(z,k[;cache]) -> [phi_0(z),phi_1(z),...,phi_k(z)]Compute the scalar phi functions for all orders up to k.
The phi functions are defined as
\[\varphi_0(z) = \exp(z),\quad \varphi_{k+1}(z) = \frac{\varphi_k(z) - 1}{z}\]
Instead of using the recurrence relation, which is numerically unstable, a formula given by Sidje is used (Sidje, R. B. (1998). Expokit: a software package for computing matrix exponentials. ACM Transactions on Mathematical Software (TOMS), 24(1), 130-156. Theorem 1).
phi(A,k[;cache]) -> [phi_0(A),phi_1(A),...,phi_k(A)]Compute the matrix phi functions for all orders up to k. k >= 1.
The phi functions are defined as
\[\varphi_0(z) = \exp(z),\quad \varphi_{k+1}(z) = \frac{\varphi_k(z) - 1}{z}\]
Calls phiv_dense on each of the basis vectors to obtain the answer. If A is Diagonal, instead calls the scalar phi on each diagonal element and the return values are also Diagonals
Caches
Missing docstring for ExpvCache. Check Documenter's build log for details.
Missing docstring for PhivCache. Check Documenter's build log for details.