I am trying to solve a 2nd order differential equation with ModelingToolkit.jl in Julia:
using ModelingToolkit
using DifferentialEquations: solve
using Plots: plot
@variables s ρ(s) Dρ(s)
@parameters N,M,Λ
D = Differential(s)
ω = sqrt(N - 1)*N/(N*(4*N*s^2 - 4*N*s + N - 4*s^2 + 4*s))
eqn = D(D(ρ)) ~ -4*ω*ρ
@named eqnODE = ODESystem([eqn],s)
I seem to get very different answers depending on whether I run
sol1 = solve(ODEProblem(eqnODE, [ρ => 0.5, D(ρ) => 0], (0.0,1.0), [N => 10000]))
or
lowered_eqnODE = ode_order_lowering(eqnODE)
sol2 = solve(ODEProblem(lowered_eqnODE, [ρ => 0.5, D(ρ) => 0], (0.0,1.0), [N => 10000]))
How is this possible? My first thought would be that I messed up the parameters, but I think I'm doing the same as here: https://mtk.sciml.ai/dev/tutorials/higher_order/
Here are the two plots:
Plot of sol1
Plot of sol2 (1st order system)
Thanks for any help!