Vern9 interpolation produces wrong results with save_idxs
#2140
Comments
|
Chaotic is extremely sensitive to numerical details (i.e. 1 ulp changes cause O(1) changes). Can you show this on a non-chaotic system? |
|
Here is an example in a non-chaotic system: using OrdinaryDiffEq
using GLMakie
l = 1.0 # length [m]
m = 1.0 # mass [kg]
g = 9.81 # gravitational acceleration [m/s²]
function pendulum!(du, u, p, t)
du[1] = u[2] # θ'(t) = ω(t)
du[2] = -3g / (2l) * sin(u[1]) + 3 / (m * l^2) * p(t) # ω'(t) = -3g/(2l) sin θ(t) + 3/(ml^2)M(t)
end
θ₀ = 0.01 # initial angular deflection [rad]
ω₀ = 0.0 # initial angular velocity [rad/s]
u₀ = [θ₀, ω₀] # initial state vector
tspan = (0.0, 10.0) # time interval
M = t -> 0.1sin(t) # external torque [Nm]
prob = ODEProblem(pendulum!, u₀, tspan, M)
sol = solve(prob, Vern9(), reltol = 1e-8, abstol = 1e-8)
sol1 = solve(prob, Vern9(), reltol = 1e-8, abstol = 1e-8, save_idxs=[1])
f = Figure(fontsize=18)
ax = Axis(f[1, 1])
plot!(ax, sol1, idxs=1)
plot!(ax, sol, idxs=1)
ax.scene.plots[2*2-1].color = Makie.wong_colors()[2]
|
|
Oh I see. This is what the warning in the docs is about: https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Explicit-Runge-Kutta-Methods
If followed it's fine: using OrdinaryDiffEq
using GLMakie
l = 1.0 # length [m]
m = 1.0 # mass [kg]
g = 9.81 # gravitational acceleration [m/s²]
function pendulum!(du, u, p, t)
du[1] = u[2] # θ'(t) = ω(t)
du[2] = -3g / (2l) * sin(u[1]) + 3 / (m * l^2) * p(t) # ω'(t) = -3g/(2l) sin θ(t) + 3/(ml^2)M(t)
end
θ₀ = 0.01 # initial angular deflection [rad]
ω₀ = 0.0 # initial angular velocity [rad/s]
u₀ = [θ₀, ω₀] # initial state vector
tspan = (0.0, 10.0) # time interval
M = t -> 0.1sin(t) # external torque [Nm]
prob = ODEProblem(pendulum!, u₀, tspan, M)
sol = solve(prob, Vern9(lazy=false), reltol = 1e-8, abstol = 1e-8)
sol1 = solve(prob, Vern9(lazy=false), reltol = 1e-8, abstol = 1e-8, save_idxs=[1])
f = Figure(fontsize=18)
ax = Axis(f[1, 1])
plot!(ax, sol1, idxs=1)
plot!(ax, sol, idxs=1)
ax.scene.plots[2*2-1].color = Makie.wong_colors()[2]This is why in the global defaulting system, With SciMLStructures.jl interface many cases of this will get fixed by having it save and reconstruct the Discrete AbstractPortion of the |
Describe the bug 🐞
When performing interpolation of the solution solved with
Vern9(), if the user specifiessave_idxs, it outputs different more discontinous interpolation results.Expected behavior
Whether or not
save_idxsis set, the interpolation results for eachidxsshould be identical.Minimal Reproducible Example 👇
Environment (please complete the following information):
using Pkg; Pkg.status()versioninfo()Additional context
Tsit5()does not have this issue. I haven't check other algorithms.The text was updated successfully, but these errors were encountered: