Large Interpolation errors for exponential integrators in SplitODEProblems

[RayleighLord]

This is a duplicate from here

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?

[oscardssmith]

Wait, looking at this, it looks like you are never assigning into the first N components of du. Isn't that a bug in your code? I've confirmed that assigning those to 0 fixes this.

[RayleighLord]

Oh, I thought that du was automatically created to be a vector of zeros. And since this ODE comes from turning a second order system to a first order, that means that the first N components are zero except for the linear part which is handled with the A matrix. Is this not the case then? Is there no elegant way of keeping those du to zero without enforcing that in every timestep?

Also that raises me one doubt. How is it possible that the ODE is solved correctly without setting it to zero but then the interpolation fails because of that?

[ChrisRackauckas]

We normally try to zero all unassigned memory instead of using similar. We may have missed one here, or in ExponentialUtilities.jl. So while I would recommend not relying on that behavior (you can get minor performance optimizations, but they are pretty minor... and normally that means you could have generated your function in a better way), we should still track down where this is coming from and zero it.