I was recently working with SplitODEProblem using exponential integrators and I noticed a huge error when doing the interpolation of the solution, which I think is wrong.
I began to notice this since, by default, the solution gives the values at each time step and always ends at tspan[2], regardless if this value is smaller than the last step for a given dt. This last interpolated step seemed to be way off from the solution, so I started experimenting with it and I have the following MWE.
using DifferentialEquations, GLMakie
K = [1.0 -1.0 0.0;
-1.0 2.0 -1.0;
0.0 -1.0 2.0]
M = [1.0 0.0 0.0;
0.0 1.0 0.0;
0.0 0.0 1.0]
F = [0.1, 0.0, 0.0]
F₀ = [20.0, 20.0, 0.0]
M⁻¹F, M⁻¹F₀, M⁻¹K = M \ F, M \ F₀, M \ K
ξ = 0.01
function f2!(du, u, p, t)
N = length(u) ÷ 2
M⁻¹F, M⁻¹F₀, Ω = p
du[(N+1):end] .= M⁻¹F * sin(Ω * t) .+ M⁻¹F₀
end
A = MatrixOperator([zeros(size(M)) I; -M⁻¹K -M⁻¹K*ξ])
Ω = 0.445
p = (M⁻¹F=M⁻¹F, M⁻¹F₀=M⁻¹F₀, Ω=Ω)
tspan = (0.0, 2π / Ω * 10)
u0 = zeros(6)
prob = SplitODEProblem(A, f2!, u0, tspan, p)
sol = solve(prob, LawsonEuler(), dt=1e-1, progress=true)
t = (2π/Ω*9):(2π/Ω/100):(2π/Ω*10)
lines(t, sol(t)[1, :])
This example creates an oscillating dynamical system and I solve it with an exponential integrator. However, if I do the interpolation as in the last lines, I get the following plot
whereas examining the complete solution from the integrator gives
which seems to be smooth and well resolved. Therefore, is this a bug of the interpolation functions for these kinds of scheme? Also, is it intended that the last time step from the time vector is always tspan[2] even though there could be a time instant that goes after this point?