StackOverflow: Second order ODE giving wrong results #751

LouisIX · 2021-05-12T19:21:20Z

This is from a post I made at stackoverflow. I'll copy/paste it here for ease of reference:

I am trying to use the DifferentialEquations.jl provided by julia, and it's working all right until I try to use it on a second order ODE. Consider for instance the second order ODE

x''(t) = x'(t) + 2* x(t), with initial conditions

x'(0) = 0, x(0) = 1

which has an analytic solution given by: x(t) = 2/3 exp(-t) + 1/3 exp(2t).

To solve it numerically, I run the following code:

`using DifferentialEquations;

function f_simple(ddu, du, u, p, t)
ddu[1] = du[1] + 2*u[1]
end;

du0 = [0.]
u0 = [1.]
tspan = (0.0,5.0)
prob2 = SecondOrderODEProblem(f_simple, du0, u0, tspan)
sol = solve(prob2,reltol=1e-8, abstol=1e-8);
`

With that,
julia> sol(3)[2] 122.57014434362732
whereas the analytic solution yields 134.50945587649028

The text was updated successfully, but these errors were encountered:

ChrisRackauckas · 2021-05-12T22:57:32Z

sol = solve(prob2,reltol=1e-8, abstol=1e-8)
sol(3)[2] # 125.01841715554878

sol = solve(prob2,reltol=1e-12, abstol=1e-12)
sol(3)[2] # 89.64575438673656

sol = solve(prob2,reltol=1e-12, abstol=1e-12, tstops = [3.0], saveat = [3.0])
sol[1][2] # 89.65888006441445

sol = solve(prob2,Tsit5(),reltol=1e-8, abstol=1e-8)
sol(3)[2] # 134.5094559950625

sol = solve(prob2,Vern7(),reltol=1e-8, abstol=1e-8)
sol(3)[2] # 134.5094558872943

Seems like it's a DPRKN6 issue as other ODE solvers work. So it's a duplicate of SciML/OrdinaryDiffEq.jl#1372 . Hmm, @YingboMa we need to dig into this.

ptoche · 2021-05-13T01:59:01Z

I chanced upon this thread. You probably know this: there is no problem solving it as a system of two linear ODEs:

  """
  The ODE system
  
      x''(t) = x'(t) + 2 x(t)
       x'(0) = 0
        x(0) = 1
  
  has solution
  
      x(t) = 2/3 exp(-t) + 1/3 exp(2t)
      
  It may also be written as
  
      x'(t) = y(t)
      y'(t) = y(t) + 2 x(t)
  
  """

  using DifferentialEquations, Plots
  
  x(t) = 2/3 * exp(-t) + 1/3 * exp(2*t)
  function p(t)
      return plot(t, x, color = "red", line = (2), label = "EXACT")
  end
  
  function ode2(ddu, du, u, p, t)
      ddu[1] = du[1] + 2 * u[1]
  end
  
  function ode1(du, u, p, t)
      du[1] = u[2]
      du[2] = u[2] + 2 * u[1]
  end
  
  tspan = (0.0, 5.0)
  
  # First-Order ODE System
  u0 = [1.0; 0.0]
  prob1 = ODEProblem(ode1, u0, tspan)
  sol1 = solve(prob1, reltol = 1e-8, abstol = 1e-8)
  
  # Second-Order ODE Equation
  du0 = [0.0]
   u0 = [1.0]
  prob2 = SecondOrderODEProblem(ode2, du0, u0, tspan)
  sol2 = solve(prob2, reltol = 1e-8, abstol = 1e-8)
  
  t = LinRange(0.0, 5.0, 100)
  
  p(t)
  plot!(sol1.t, sol1[1,:], color = "blue", line = (:dash,2), label = "FIRST-ORDER")
  plot!(sol2.t, sol2[1,:], color = "gray", line = (2), label = "SECOND-ORDER")

ptoche · 2021-05-13T05:05:55Z

As the tolerance is increased, the method switches from DPRKN6 to DPRKN8 and to DPRKN12, with steps something like 1e-6 and 1e-9, if I'm understanding correctly. So I'm guessing at 1e-8 it's DPRKN8 and at 1e-12 it's DPRKN12.

DifferentialEquations.jl/src/ode_default_alg.jl

Line 22 in 975e1a2

if tol_level == :low_tol || tol_level == :med_tol

Without changing the tolerance, performance isn't great:

  p(t)
  prob3 = SecondOrderODEProblem(ode2, du0, u0, tspan)
  sol3a = solve(prob3, DPRKN6(), reltol = 1e-8, abstol = 1e-8)
  plot!(sol3a.t, sol3a[1,:], color = "blue", label = "DPRKN6")
  sol3b = solve(prob3, DPRKN8(), reltol = 1e-8, abstol = 1e-8)
  plot!(sol3b.t, sol3b[1,:], color = "green", label = "DPRKN8")
  sol3c = solve(prob3, DPRKN12(), reltol = 1e-8, abstol = 1e-8)
  plot!(sol3c.t, sol3c[1,:], color = "orange", label = "DPRKN12")

But even increasing tolerance doesn't seem to do enough:

  p(t)
  sol3g = solve(prob2, DPRKN6(),reltol = 1e-6, abstol = 1e-6)
  plot!(sol3g.t, sol3g[1,:], color = "blue", label = "DPRKN6 + TOL=1e-6")
  sol3h = solve(prob2, DPRKN8(),reltol = 1e-12, abstol = 1e-12)
  plot!(sol3h.t, sol3h[1,:], color = "green", label = "DPRKN8 + TOL=1e-12")
  sol3i = solve(prob2, DPRKN12(),reltol = 1e-20, abstol = 1e-20)
  plot!(sol3i.t, sol3i[1,:], color = "gray", label = "DPRKN12 + TOL=1e-20")

And unless I'm completely misunderstanding, other methods have similar outcome

  p(t)
  sol3e = solve(prob2, Tsit5(),reltol = 1e-8, abstol = 1e-8)
  plot!(sol3e.t, sol3e[1,:], color = "blue", label = "TSIT5")
  sol3f = solve(prob2, Vern7(),reltol = 1e-8, abstol = 1e-8)
  plot!(sol3f.t, sol3f[1,:], color = "orange", line = :dash, label = "VERN7")

While I haven't checked it manually, I'm assuming the exact solution is correct, because the solution for the first-order linear system matches the functional form (see the first graph). Turning off 'adaptive' didn't seem to help.

ChrisRackauckas · 2021-05-13T07:12:49Z

Turning off adaptive seems to do fine:

using DifferentialEquations

function f_simple(ddu, du, u, p, t)
    ddu[1] = du[1] + 2*u[1]
end

du0 = [0.]
u0 = [1.]
tspan = (0.0,5.0)
prob2 = SecondOrderODEProblem(f_simple, du0, u0, tspan)

sol = solve(prob2,DPRKN6(),dt=0.00001,adaptive=false)
sol[end][2] # 134.5067059149833

This makes sense too because it would explain why the convergence tests didn't catch it

And unless I'm completely misunderstanding, other methods have similar outcome

Since it is exponentially increasing, you do get exponentially increasing error. But the test case is reasonably well-behaved at 3.0:

t = 3.0
2/3 * exp(-t) + 1/3 * exp(2t) # 134.50945587649028

and all of the first order methods hit it, so it's a good test. At t=5, not so much because of the exponential growth of the errors though.

ChrisRackauckas · 2021-05-13T07:48:20Z

I hand checked the adaptivity parameters:

sum_b = 329//4212 + (8411_9543+366_727R)/4096_22616 + (8411_9543-366_727R)/4096_22616 + 200//17901 # 0.5
sum_bhat = (2701+23R)/4563 - (9829+131R)/9126 + 5(1798+17R)/9126 # 0.5
sum_bphat = 115//2106 + (8411_9543+366_727R)/2560_14135 + (8411_9543-366_727R)/2560_14135 + 6950//17901 - 1//10 # 1.0

Everything looks fine. And it looks like it converges fine too:

using DifferentialEquations

function f_simple(ddu, du, u, p, t)
    ddu[1] = du[1] + 2*u[1]
end

du0 = [0.]
u0 = [1.]
tspan = (0.0,5.0)
prob2 = SecondOrderODEProblem(f_simple, du0, u0, tspan)

tspan = (0.0,3.0)
prob2 = SecondOrderODEProblem(f_simple, du0, u0, tspan)
sol = solve(prob2,DPRKN6(),dt=0.00001,adaptive=false)
sol[end][2] # 134.5067059149833

sol = solve(prob2,DPRKN6(),reltol=1e-14, abstol=1e-14)
sol[end][2] # 133.8237798051805

sol = solve(prob2,DPRKN6(),reltol=1e-16, abstol=1e-16)
sol[end][2] # 134.23469439807363

tspan = (0.0,5.0)
prob2 = SecondOrderODEProblem(f_simple, du0, u0, tspan)
sol = solve(prob2,DPRKN6(),reltol=1e-16, abstol=1e-16)
sol(3)[2] # 134.23469439807363

Even DPRKN12 displays the same behavior:

using DifferentialEquations

function f_simple(ddu, du, u, p, t)
    ddu[1] = du[1] + 2*u[1]
end

du0 = [0.]
u0 = [1.]
tspan = (0.0,5.0)
prob2 = SecondOrderODEProblem(f_simple, du0, u0, tspan)

tspan = (0.0,3.0)
prob2 = SecondOrderODEProblem(f_simple, du0, u0, tspan)
sol = solve(prob2,DPRKN12(),dt=0.00001,adaptive=false)
sol[end][2] # 134.5067059149833

sol = solve(prob2,DPRKN12(),reltol=1e-14, abstol=1e-14)
sol[end][2] # 133.8237798051805

sol = solve(prob2,DPRKN12(),reltol=1e-16, abstol=1e-16)
sol[end][2] # 123.79520275557033

sol = solve(prob2,DPRKN12(),reltol=1e-18, abstol=1e-18)
sol[end][2] # 134.15180697251196

sol = solve(prob2,DPRKN12(),reltol=1e-20, abstol=1e-20)
sol[end][2] # 134.5071611034424

tspan = (0.0,5.0)
prob2 = SecondOrderODEProblem(f_simple, du0, u0, tspan)
sol = solve(prob2,DPRKN12(),reltol=1e-20, abstol=1e-20)
sol(3)[2] # 134.50716110344243

The first order methods seem to just be more stable:

using DifferentialEquations

function f_simple(ddu, du, u, p, t)
    ddu[1] = du[1] + 2*u[1]
end

du0 = [0.]
u0 = [1.]
tspan = (0.0,5.0)
prob2 = SecondOrderODEProblem(f_simple, du0, u0, tspan)

tspan = (0.0,3.0)
prob2 = SecondOrderODEProblem(f_simple, du0, u0, tspan)
sol = solve(prob2,Tsit5(),dt=0.00001,adaptive=false)
sol[end][2] # 134.50945587348434

sol = solve(prob2,Tsit5(),reltol=1e-8, abstol=1e-8)
sol[end][2] # 134.5094559698895

sol = solve(prob2,Tsit5(),reltol=1e-14, abstol=1e-14)
sol[end][2] # 134.50945587649144

tspan = (0.0,5.0)
prob2 = SecondOrderODEProblem(f_simple, du0, u0, tspan)
sol = solve(prob2,Tsit5(),reltol=1e-14, abstol=1e-14)
sol(3)[2] # 134.50945587649144

tspan = (0.0,3.0)
prob2 = SecondOrderODEProblem(f_simple, du0, u0, tspan)
sol = solve(prob2,Vern9(),dt=0.001,adaptive=false)
sol[end][2] # 134.50945587654974

sol = solve(prob2,Vern9(),reltol=1e-8, abstol=1e-8)
sol[end][2] # 134.50945587935112

sol = solve(prob2,Vern9(),reltol=1e-14, abstol=1e-14)
sol[end][2] # 134.50945587649028

tspan = (0.0,5.0)
prob2 = SecondOrderODEProblem(f_simple, du0, u0, tspan)
sol = solve(prob2,Vern9(),reltol=1e-14, abstol=1e-14)
sol(3)[2] # 134.50945587648263

ChrisRackauckas · 2021-05-13T08:00:03Z

Looking again at , it looks like it's just unstable at 1e-8:

using DifferentialEquations
using Plots

b = range(1, 3, length=6)

t0 = -100.
t1 = 150.
tspan = (t0, t1)
i = 1

function hill!(ddu, du, u, p, t)
    f = -3/(u[1]^2 + u[2]^2)^1.5

    ddu[1] = u[1]*(f + 3) + 2du[2]
    ddu[2] = u[2]*f - 2du[1]
end
u0 = zeros(2)
du0 = zeros(2)

u0[1] = b[i] - 8/(3b[i]^2 * (t0^2 + 4/9)^0.5)
u0[2] = -1.5*u0[1]*t0
du0[1] = 0
du0[2] = -1.5*u0[1]

prob = SecondOrderODEProblem(hill!, du0, u0, tspan, abstol=1e-8, reltol=1e-8)

sol = solve(prob,DPRKN6())
plot(sol)

sol = solve(prob,Tsit5())
plot(sol)

and at a lower tolerance you can see how it manifests from instability:

prob = SecondOrderODEProblem(hill!, du0, u0, tspan, abstol=1e-12, reltol=1e-12)

sol = solve(prob,DPRKN6())
plot(sol)
savefig("plot1.png")

sol = solve(prob,Tsit5())
plot(sol)
savefig("plot2.png")

Before it fully dampens:

prob = SecondOrderODEProblem(hill!, du0, u0, tspan, abstol=1e-14, reltol=1e-14)

sol = solve(prob,DPRKN6())
plot(sol)
savefig("plot1.png")

sol = solve(prob,Tsit5())
plot(sol)
savefig("plot2.png")

All of this plus the manual tableau checks seem to indicate that it's not an implementation issue, but rather it seems to be a fatal flaw in the method. Both the stepper and the embedded method seem to have small stability regions and this issue is amplified on some problems. Trying to dig into it, I see https://www.sciencedirect.com/science/article/pii/0898122187900666 the paper doesn't even reference stability (😱), so it may have been a fatal overlooked issue in the method. RKN methods may still be a ripe area for researching better methods.

I'll let @YingboMa take a dig into this and see if he finds anything I didn't see, but my conclusion is that the DPRKN series of methods with local extrapolation, they mention in the paper that you should do local extrapolation, but I don't think it was well-tested as to how stable it truly is, and I don't think that was a mistake. So while still useful on many problems, I can see how users are finding some problematic cases, and am thinking of changing the default over to using the first order methods.

ChrisRackauckas · 2021-05-13T08:02:31Z

A property that can give rise to this behavior is if the stepping method is much more stable than the embedded method, making the error estimates unstable even when the method is okay, making them a bit wild and unpredictable on the numerically difficult problems. Some stiff ODE solvers have displayed this behavior (where the embedded method isn't even A-stable), so it could be a similar phenomena.

ptoche · 2021-05-13T20:10:57Z

Thank you for sharing this super detailed analysis. Is there a case for having a check on linearity of first- and second-order problems and quite simply hard-coding exact solutions in these cases? It would add a bit of overhead, but then you don't need to worry about what the default solver is in these cases. Users will expect linear systems to "just work". Just a thought.

ChrisRackauckas · 2021-05-14T14:40:43Z

I recently added that to Pumas. Making it work automatically on all of DifferentialEquations.jl would require a few more assumptions, which we might want better compiler tools for. But basically the way to do it is:

https://openreview.net/pdf?id=rJlPdcY38B

calculate the sparsity of the Hessian, and then that the sparsity pattern with respect to t is empty. In Symbolics.jl terms:

https://github.com/JuliaSymbolics/Symbolics.jl/blob/7cfe405b7def6cde138c0bb19ef3effaa1209e45/src/diff.jl#L481

Using the same ModelingToolkit tricks shown in https://arxiv.org/abs/2103.05244, we can automatically apply it to any quasi-static code, but fully imperative codes would need a fallback.

ptoche · 2021-05-14T22:53:18Z

Fantastic project!

You just taught me the difference between affine and linear in English. For some reason, the meanings are different in France (where linear means proportional, affine means linear!!): wiki/Fonction_affine

ChrisRackauckas · 2021-05-25T09:38:53Z

Defaults were changed to avoid DPRKN

HenryLangner · 2023-05-01T15:46:36Z

Another possibility for the wrong results may be due to the construction of the DPRKN methods. They are designed for 2nd order ODEs which are not dependent on the first derivative. This could have lead to the wrong results. More information and some testing can be found here:

SciML/OrdinaryDiffEq.jl#1938

ChrisRackauckas closed this as completed May 25, 2021

HenryLangner mentioned this issue Apr 27, 2023

Repeated occurrence of Runge-Kutta-Nyström issues concerning their accuracy SciML/OrdinaryDiffEq.jl#1938

Open

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

StackOverflow: Second order ODE giving wrong results #751

StackOverflow: Second order ODE giving wrong results #751

LouisIX commented May 12, 2021

ChrisRackauckas commented May 12, 2021

ptoche commented May 13, 2021 •

edited

Loading

ptoche commented May 13, 2021 •

edited

Loading

ChrisRackauckas commented May 13, 2021

ChrisRackauckas commented May 13, 2021

ChrisRackauckas commented May 13, 2021

ChrisRackauckas commented May 13, 2021

ptoche commented May 13, 2021 •

edited

Loading

ChrisRackauckas commented May 14, 2021

ptoche commented May 14, 2021

ChrisRackauckas commented May 25, 2021

HenryLangner commented May 1, 2023

StackOverflow: Second order ODE giving wrong results #751

StackOverflow: Second order ODE giving wrong results #751

Comments

LouisIX commented May 12, 2021

ChrisRackauckas commented May 12, 2021

ptoche commented May 13, 2021 • edited Loading

ptoche commented May 13, 2021 • edited Loading

ChrisRackauckas commented May 13, 2021

ChrisRackauckas commented May 13, 2021

ChrisRackauckas commented May 13, 2021

ChrisRackauckas commented May 13, 2021

ptoche commented May 13, 2021 • edited Loading

ChrisRackauckas commented May 14, 2021

ptoche commented May 14, 2021

ChrisRackauckas commented May 25, 2021

HenryLangner commented May 1, 2023

ptoche commented May 13, 2021 •

edited

Loading

ptoche commented May 13, 2021 •

edited

Loading

ptoche commented May 13, 2021 •

edited

Loading