`sol(t, Val{2})` on FBDF includes spurious discontinuities for (differential) variable in DAE system

[topolarity]

using ModelingToolkit

@variables t
D = Differential(t)

@mtkmodel Foo begin
    @parameters begin
        ω₁
        ω₂
        l
    end
    @variables begin
        x(t) # differential variable
        y(t) # algebraic variable
    end
    @equations begin
        D(D(x)) ~ -ω₁^2 * x     # solution: x(t) = cos(ω₁t)
        y^l     ~ sin(ω₂ * t)   # solution: y(t) = sin(ω₂t) when l == 1.0
    end
end

@mtkbuild foo = Foo();

prob = ODEProblem(foo,
                  [foo.x => 1.0, D(foo.x) => 0.0, foo.y => 0.0],
                  (0.0, 2pi),
                  [foo.ω₁ => 1.0, foo.ω₂ => 10.0, foo.l => 1.0]);
               
using OrdinaryDiffEq: solve, FBDF
using Plots

sol = solve(prob, FBDF())

t = 0.01:1e-5:2pi
plot(t, sol(t, Val{0}, idxs=foo.x).u; label="primal")
plot!(t, sol(t, Val{1}, idxs=foo.x).u; label="first-derivative")
plot!(t, sol(t, Val{2}, idxs=foo.x).u; label="second-derivative")

yields the following plot:

The discontinuities appear to align with the zero-crossings of the algebraic variable:

This gets worse if foo.ω₂ is increased to 100.0:

[staticfloat]

We looked into this a bit more and found that there are actually issues in both the first derivative and the primal, just much smaller (but in the case of the first derivative, still much larger than abstol).

Code:

using OrdinaryDiffEq: solve, FBDF
using WGLMakie

deftol = 1e-8
sol = solve(prob, FBDF(); abstol=deftol, reltol=deftol)

t = 0.01:1e-5:2pi

# Signal plots
begin
    fig = Figure()
    ax = Axis(fig[1,1]; title="Signals")
    lines!(ax, t, sol(t, Val{0}, idxs=foo.x).u; label="primal")
    lines!(ax, t, sol(t, Val{1}, idxs=foo.x).u; label="first derivative")
    lines!(ax, t, sol(t, Val{2}, idxs=foo.x).u; label="second derivative")
    axislegend(ax)
    fig
end

# Error plots
begin
    fig = Figure()
    ax = Axis(fig[1,1]; title="Normalized Error")
    primal_err = sol(t, Val{0}, idxs=foo.x).u .- cos.(t)
    lines!(ax, t, primal_err./maximum(abs.(primal_err)); label="primal")

    d1_err = sol(t, Val{1}, idxs=foo.x).u .+ sin.(t)
    lines!(ax, t, d1_err./maximum(abs.(d1_err)); label="first-derivative")

    d2_err = sol(t, Val{2}, idxs=foo.x).u .+ cos.(t)
    lines!(ax, t, d2_err./maximum(abs.(d2_err)); label="second-derivative")
    axislegend(ax)
    fig
end

Plots:

[ChrisRackauckas]

Is specific to FBDF:

using ModelingToolkit

@variables t
D = Differential(t)

@mtkmodel Foo begin
    @parameters begin
        ω₁
        ω₂
        l
    end
    @variables begin
        x(t) # differential variable
        y(t) # algebraic variable
    end
    @equations begin
        D(D(x)) ~ -ω₁^2 * x     # solution: x(t) = cos(ω₁t)
        y^l     ~ sin(ω₂ * t)   # solution: y(t) = sin(ω₂t) when l == 1.0
    end
end

@mtkbuild foo = Foo();

prob = ODEProblem(foo,
                  [foo.x => 1.0, D(foo.x) => 0.0, foo.y => 0.0],
                  (0.0, 2pi),
                  [foo.ω₁ => 1.0, foo.ω₂ => 10.0, foo.l => 1.0]);

using OrdinaryDiffEq
using Plots

sol = solve(prob, RadauIIA5(), abstol=1e-8, reltol=1e-8)
t = 0.01:1e-4:2pi
plot(t, sol(t, Val{2}, idxs=foo.x).u)
savefig("plot1.png")

sol = solve(prob, QNDF(), abstol=1e-8, reltol=1e-8)
t = 0.01:1e-4:2pi
plot(t, sol(t, Val{2}, idxs=foo.x).u)
savefig("plot2.png")

sol = solve(prob, FBDF(), abstol=1e-8, reltol=1e-8)
t = 0.01:1e-4:2pi
plot(t, sol(t, Val{2}, idxs=foo.x).u)
savefig("plot3.png")

[ChrisRackauckas]

f(integrator.fsallast, u, p, tdt)

FBDF is slightly lying about its derivative. Or rather, the first derivative is "correct" for differential variables, but the interpolation polynomial of the FBDF method is different from the interpolation polynomial that we are showing from its solution. We should populate fsallast with the first derivative estimate from its internal polynomial interpolation. This would also make algebraic variable interpolation work for this method in principle.