Other algorithms

[jagot]

Is this package the place where we'd like to put other algorithms as well? Such as

• Chebyshev (requires prior knowledge of spectral range of the operator)
• Padé (explicit & implicit Euler, Crank–Nicolson).
• Spectral (when the operator is diagonal in the basis chosen, or even if you have precomputed the diagonalization)

I imagine being able to do the following.

CrankNicolson{T} where T = Padé{1,1,T}

A = SymTridiagonal(...) # or whatever
# For a time-independent, this would precompute the factorizations necessary, that is I - A/2
# This requires fixed time-step, though
Ae = CrankNicolson(dt*A)

expv!(w, Ae, b) # Amounts to w = (I-A/2)\(I+A/2)*b

# Alternative interface, for time-dependent operators/adaptive time-steps, which cannot reuse previous factorizations
expv!(w, dt, A, P)

[jagot]

To illustrate where I'd like to go, I want to be able to write splitting algorithms in the following way (approximately):

macro swap!(x,y)
     quote
         local tmp = $(esc(x))
         $(esc(x)) = $(esc(y))
         $(esc(y)) = tmp
     end
end

function strang_splitting(fun::Function,
                           steps::Integer,
                           dt::AbstractFloat,
                           μ::Number,
                           v::V,
                           As::Vector) where V
     A₁ = As[1]
     tmp = similar(v)
     a,b = 1,2
     ws = (v,tmp)
     τ = μ*dt
     
     for i in 1:steps
         t = dt*(i-0.5) # 2nd order splitting, evaluate at midpoint of time-step
         for A in reverse(As[2:end])
             expv!(ws[b], τ/2, A(t), ws[a])
             @swap! a b
         end

         expv!(ws[b], τ, A₁(t), ws[a])
         @swap! a b

         for A in As[2:end]
             expv!(ws[b], τ/2, A(t), ws[a])
             @swap! a b
         end
     end
     copyto!(v, tmp)
 end

That is, the splitting algorithm should not need to bother how the exponentiation is performed, which caches are necessary, etc. For time-independent A, and if the time-step is the same as previously used, the exponentiator should be able to figure out it can keep previously computed factorizations.

[ChrisRackauckas]

I think any matrix exponential method which can be used in exponential integrators is worth having here. The integrators themselves go to OrdinaryDiffEq.jl

[MSeeker1340]

Expokit.jl (or the original Expokit) can be a good reference for Pade and Chebyshev.

@jagot Do you have any special idea for spectral algorithms? Since we already have Diagonal handling capability in phi (https://github.com/JuliaDiffEq/ExponentialUtilities.jl/blob/master/src/phi.jl#L139) I think this should satisfy all that we need.

[jagot]

Spectral I just added for completeness. It is good that Diagonal is already covered, but I figured it would be useful to have an expv! implementation for the case when you've explicitly pre-diagonalized your operator: exp(A) = S*exp(D)*S', where S is the matrix of eigenvectors, and D the diagonal matrix of eigenvalues.