Open jagot opened 5 years ago
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.
I think any matrix exponential method which can be used in exponential integrators is worth having here. The integrators themselves go to OrdinaryDiffEq.jl
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.
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.
Is this package the place where we'd like to put other algorithms as well? Such as
I imagine being able to do the following.