# Strange solution with DifferentialEquations.jl

**URL:** <https://discourse.julialang.org/t/strange-solution-with-differentialequations-jl/109579>\
**Category:** General Usage\
**Tags:** plotting, symbolics, differentialequation\
**Created:** [February 1, 2024, 3:45pm UTC](https://discourse.julialang.org/t/strange-solution-with-differentialequations-jl/109579 "2024-02-01T15:45:32Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![Sh\_166](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sh_166/32/43965_2.png) [@Sh\_166](https://discourse.julialang.org/u/Sh_166)\
**Post date:** [February 1, 2024, 3:45pm UTC](https://discourse.julialang.org/t/strange-solution-with-differentialequations-jl/109579/1 "2024-02-01T15:45:32Z")

</div>

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

MWE:

```import
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:

 ![Screenshot 2024-02-01 210707](https://global.discourse-cdn.com/julialang/original/3X/4/5/456971b76ca09fe2a99b569df1bbc877515c07ed.jpeg)

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!

---

<div class="post-metadata">

**Author:** ![Sh\_166](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sh_166/32/43965_2.png) [@Sh\_166](https://discourse.julialang.org/u/Sh_166)\
**Post date:** [February 2, 2024, 5:05am UTC](https://discourse.julialang.org/t/strange-solution-with-differentialequations-jl/109579/2 "2024-02-02T05:05:13Z")

</div>

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:

```julia
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:

```julia
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:

```julia
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`.

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [February 2, 2024, 5:52am UTC](https://discourse.julialang.org/t/strange-solution-with-differentialequations-jl/109579/3 "2024-02-02T05:52:51Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Sh\_166](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sh_166/32/43965_2.png) [@Sh\_166](https://discourse.julialang.org/u/Sh_166)\
**Post date:** [February 2, 2024, 6:01am UTC](https://discourse.julialang.org/t/strange-solution-with-differentialequations-jl/109579/4 "2024-02-02T06:01:16Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [February 2, 2024, 6:11am UTC](https://discourse.julialang.org/t/strange-solution-with-differentialequations-jl/109579/5 "2024-02-02T06:11:43Z")

</div>

Yes they are both equivalent to:

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

```

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [February 2, 2024, 6:30am UTC](https://discourse.julialang.org/t/strange-solution-with-differentialequations-jl/109579/6 "2024-02-02T06:30:13Z")

</div>

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.

```julia
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.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [February 2, 2024, 7:33am UTC](https://discourse.julialang.org/t/strange-solution-with-differentialequations-jl/109579/7 "2024-02-02T07:33:47Z")

</div>

You can probably use this solver [GitHub - pnavaro/HOODESolver.jl: High Oscillatory Ordinary Differential Equation Solver in Julia](https://github.com/pnavaro/HOODESolver.jl)

---

<div class="post-metadata">

**Author:** ![Sh\_166](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sh_166/32/43965_2.png) [@Sh\_166](https://discourse.julialang.org/u/Sh_166)\
**Post date:** [February 2, 2024, 7:44am UTC](https://discourse.julialang.org/t/strange-solution-with-differentialequations-jl/109579/8 "2024-02-02T07:44:48Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [February 2, 2024, 8:01am UTC](https://discourse.julialang.org/t/strange-solution-with-differentialequations-jl/109579/9 "2024-02-02T08:01:07Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Sh\_166](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sh_166/32/43965_2.png) [@Sh\_166](https://discourse.julialang.org/u/Sh_166)\
**Post date:** [February 2, 2024, 9:24am UTC](https://discourse.julialang.org/t/strange-solution-with-differentialequations-jl/109579/10 "2024-02-02T09:24:46Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Sh\_166](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sh_166/32/43965_2.png) [@Sh\_166](https://discourse.julialang.org/u/Sh_166)\
**Post date:** [February 2, 2024, 12:39pm UTC](https://discourse.julialang.org/t/strange-solution-with-differentialequations-jl/109579/11 "2024-02-02T12:39:12Z")

</div>

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

```julia
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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 2, 2024, 12:50pm UTC](https://discourse.julialang.org/t/strange-solution-with-differentialequations-jl/109579/12 "2024-02-02T12:50:20Z")

</div>

> <https://scicomp.stackexchange.com/questions/29149/what-does-symplectic-mean-in-reference-to-numerical-integrators-and-does-scip>

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.

---

<div class="post-metadata">

**Author:** ![Sh\_166](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sh_166/32/43965_2.png) [@Sh\_166](https://discourse.julialang.org/u/Sh_166)\
**Post date:** [February 2, 2024, 5:36pm UTC](https://discourse.julialang.org/t/strange-solution-with-differentialequations-jl/109579/13 "2024-02-02T17:36:47Z")

</div>

@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.