Merge pull request #378 from MSeeker1340/krylov-expv-timestep

Time-stepping Krylov: special case for `expv_timestep`

src/exponential_utils.jl

@@ -439,6 +439,47 @@ end
 
 ###########################################
 # Krylov phiv with internal time-stepping
+"""
+    exp_timestep(t,A,b[;adaptive,tol,kwargs...]) -> u
+
+Evaluates the matrix exponentiation-vector product using time stepping
+
+```math
+u = \\exp(tA)b 
+```
+
+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.
+
+Note that this function is just a special case of `phiv_timestep` with a more 
+intuitive interface (vector `b` instead of a n-by-1 matrix `B`).
+
+[^1]: Niesen, J., & Wright, W. (2009). A Krylov subspace algorithm for 
+evaluating the φ-functions in exponential integrators. arXiv preprint 
+arXiv:0907.4631.
+"""
+function expv_timestep(t, A, b; kwargs...)
+  u = Vector{eltype(A)}(size(A, 1))
+  expv_timestep!(u, t, A, b; kwargs...)
+end
+"""
+    expv_timestep!(u,t,A,b[;kwargs]) -> u
+
+Non-allocating version of `expv_timestep`.
+"""
+function expv_timestep!(u::Vector{T}, t::Real, A, b::Vector{T}; 
+  kwargs...) where {T <: Number}
+  B = reshape(b, length(b), 1)
+  phiv_timestep!(u, t, A, B; kwargs...)
+end
 """
     phiv_timestep(t,A,B[;adaptive,tol,kwargs...]) -> u
 

test/utility_tests.jl

@@ -1,4 +1,4 @@
-using OrdinaryDiffEq: phi, phi, phiv, phiv_timestep, expv, arnoldi, getH, getV
+using OrdinaryDiffEq: phi, phi, phiv, phiv_timestep, expv, expv_timestep, arnoldi, getH, getV
 
 @testset "Exponential Utilities" begin
   # Scalar phi
@@ -57,4 +57,8 @@ using OrdinaryDiffEq: phi, phi, phiv, phiv_timestep, expv, arnoldi, getH, getV
   u_exact = sum(t^i * Phi[i+1] * B[:,i+1] for i = 0:K)
   u = phiv_timestep(t, A, B; adaptive=true, tol=tol)
   @test norm(u - u_exact) / norm(u_exact) < tol
+  # p = 0 special case (expv_timestep)
+  u_exact = Phi[1] * B[:, 1]
+  u = expv_timestep(t, A, B[:, 1]; adaptive=true, tol=tol)
+  @test norm(u - u_exact) / norm(u_exact) < tol
 end