Strange solution with DifferentialEquations.jl

[Sh_166]

I have run into a weird problem with DifferentialEquations.jl.

MWE:

Pkg.add("Symbolics")
Pkg.add("Plots")
Pkg.add("DifferentialEquations")

using Symbolics, Plots, DifferentialEquations

function trial(dx, x, p, t)
    a = p 
    dx[1] = -im*a*x[1]
    dx[2] = -im*a*x[2]
end

u_0 = [1.0+0.0*im, 0.0]
P = 1e5

PR = ODEProblem(trial, u_0, (0.0, 100.0), P, saveat = 0.1)
S = solve(PR, maxiters = 1e8)

f(t, x, y) = (t, (abs(x)^2 + abs(y)^2)^2)
plot(S, idxs = (f, 0, 1, 2))

This system should just plot a straight line at 1 - but what plot (and the solve() output) show is this:

This is strange. I do not need to set P = 1e5 to get this, even 1 shows the decrease (but not as steep - so this value enhances the effect). The saveat portion is also not the culprit because I replicated this even without it. Finally, I put the maxiters there because the output showed Error: Large maxiters needed (or some variation of it).

It probably is due to some syntax issue, so any help will be appreciated!

[Sh_166]

Update: I replicated this issue on both Julia 1.9.2 and Julia 1.8.2 (on both my usual JupyterLab interface and the basic Julia IDE). I also included a stiff solver Rodas4P(autodiff=false), but it persists. It seems to be a simple differential equation and the code throws up no errors either, so I am a bit confused as to what is going wrong here.

The weirder thing is that the following code replicates the issue:

Pkg.add("Symbolics")
Pkg.add("Plots")
Pkg.add("DifferentialEquations")

using Symbolics, Plots, DifferentialEquations

function trial(dx, x, p, t)
    a = p 
    dx[1] = -im*a*x[1]
end

u_0 = [1.0+0.0*im]
P = 1e5

PR = ODEProblem(trial, u_0, (0.0, 100.0), P, saveat = 0.1)
S = solve(PR, maxiters = 1e8)

f(t, x) = (t, abs(x)^2)
plot(S, idxs = (f, 0, 1))

But this one does not:

Pkg.add("Symbolics")
Pkg.add("Plots")
Pkg.add("DifferentialEquations")

using Symbolics, Plots, DifferentialEquations

function trial(dx, x, p, t)
    a = p 
    dx = -im*a*x
end

u_0 = [1.0+0.0*im]
P = 1e5

PR = ODEProblem(trial, u_0, (0.0, 100.0), P, saveat = 0.1)
S = solve(PR, maxiters = 1e8)

f(t, x) = (t, abs(x)^2)
plot(S, idxs = (f, 0, 1))

I played around a bit and this one also gets the right solution:

Pkg.add("Symbolics")
Pkg.add("Plots")
Pkg.add("DifferentialEquations")

using Symbolics, Plots, DifferentialEquations

function trial(du, u, p, t)
    a = p 
    x, y = u
    dx = -im*a*x
    dy = -im*a*y
end

u_0 = [1.0+0.0*im, 0.0]
P = 1e5

PR = ODEProblem(trial, u_0, (0.0, 100.0), P, saveat = 0.1)
S = solve(PR, Rodas4P(autodiff=false))

f(t, x, y) = (t, (abs(x)^2 + abs(y)^2)^2)
plot(S, idxs = (f, 0, 1, 2))

However, not specifying the solver leads to a different kind of wrong solution.

Update 2: I tried doing this for my original system of equations (which was more complicated and led to this system of equations in a limit) - it now shows a constant line at 1 (but with a random assortment of small spikes). I would like to remove those spikes as well, but they’re too small for the tick resolution and so, for now, I am working around the issue by specifying ylims.

[abraemer]

Your last two versions don’t overwrite the contents of the du so don’t do anything.

I wonder whether the issue here is, that your solution oscillates very quickly. I’ll if I can investigate.

[Sh_166]

Thanks for the reply - so, do you mean I should not use those last 2 versions?

[abraemer]

Yes they are both equivalent to:

function trial(du, u, p, t)
    nothing
end

[abraemer]

Ok so the main problem is as I expected: Your problem oscillates rapidly and the default tolerances are quite loose. If you set e.g. P=1 it will work but slowly and you will be off by half a percent at t=100 (I used Tsit5 as solver). If you additionally set some lower tolerance, e.g.

solve(PR, maxiters = 1e8, abstol=1e-8, reltol=1e-8)

The integration runs almost instantly and the values stay very close to 1. Still P=1e5 does not really work but is that a realistic scenario for you? You probably need different methods then.

[rveltz]

You can probably use this solver GitHub - pnavaro/HOODESolver.jl: High Oscillatory Ordinary Differential Equation Solver in Julia

[Sh_166]

Unfortunately, I do need the 1e5 choice, it is actually a physical parameter. I’ll check out the link @rveltz, thanks a lot to both of you!

[abraemer]

And also a time of 100? That are ~3e6 oscillation periods in total. To resolve that you will a lot of time steps in any case. And play with the tolerances. Lower tolerance does not automatically mean longer solve time. There is usually a sweet spot somewhere.

[Sh_166]

Yeah, in my actual problem its not time, but distance - so, I am easily going up to larger length scales than 100. Never had an issue with it till now, since I was more concerned about very small parameter values.

[Sh_166]

@rveltz Could you help me out a bit with HOODESolver.jl? I basically copied this from the documentation:

using Symbolics, DifferentialEquations, Plots, HOODESolver

ϵ = 0.0001

A = [ -im 0 ;
      0 -im ]

f1 = LinearHOODEOperator(ϵ, A)

f2 = (u,p,t) ->  [0, 0]

tspan = (0.0, 100.0)
u0 = [1.0+0.0*im, 0.0]

prob = SplitODEProblem(f1, f2, u0, tspan) # seems to run fine till this point

S = solve(prob, HOODEAB(), dt=0.1)

This last step throws up an error:

MethodError: no method matching HOODEProblem(::var"#9#10", ::Vector{ComplexF64}, ::Tuple{Float64, Float64}, ::Missing, ::Matrix{Complex{Int64}}, ::Float64, ::Missing)
Closest candidates are:
  HOODEProblem(::Any, ::Any, ::Any, ::Any, ::Any, ::Any) at C:\Users\shiha\.julia\packages\HOODESolver\5dr9W\src\interface.jl:70
  HOODEProblem(::Any, ::Vector{T}, ::Tuple{T, T}, ::Any, ::AbstractMatrix, ::T, ::Union{Missing, Matrix}) where T at C:\Users\shiha\.julia\packages\HOODESolver\5dr9W\src\interface.jl:47

Stacktrace:
 [1] HOODEProblem(f::Function, u0::Vector{ComplexF64}, tspan::Tuple{Float64, Float64}, p::Missing, A::Matrix{Complex{Int64}}, epsilon::Float64)
   @ HOODESolver C:\Users\shiha\.julia\packages\HOODESolver\5dr9W\src\interface.jl:71
 [2] solve(prob::ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, SplitFunction{false, SciMLBase.FullSpecialize, ODEFunction{false, SciMLBase.FullSpecialize, LinearHOODEOperator{Float64}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, var"#9#10", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SplitODEProblem{false}}, alg::HOODEAB{4, 32}; dt::Int64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ HOODESolver C:\Users\shiha\.julia\packages\HOODESolver\5dr9W\src\common_interface.jl:107
 [3] top-level scope
   @ In[49]:1

I tried a couple of things, but can’t seem to figure it out.

[ChrisRackauckas]

What does "symplectic" mean in reference to numerical integrators, and does SciPy's odeint use them?

Basically, standard ODE solvers have drift on highly oscillatory problems and you need to treat the problem differently if you have very long integrations. Or set the tolerance really low.

[Sh_166]

@ChrisRackauckas So, this means that I’ll possibly need to adjust the tolerances and see what works - mine isn’t a Hamiltonian system in the sense needed for the symplectic integrators. Thanks for the explanation, though - now I know what the issue was.