Identical differential equations give different solutions when an independent equation is added

General Usage
, ,
1

I define the Lorenz system twice, once the standard way and once with an extra equation that depends only on u_1.

using DifferentialEquations
using Makie, CairoMakie

function lorenz!(du, u, p, t)
    du[1] = 10.0 * (u[2] - u[1])
    du[2] = u[1] * (28.0 - u[3]) - u[2]
    du[3] = u[1] * u[2] - (8 / 3) * u[3]
end

u0 = [1.0; 0.0; 0.0]
tspan = (0.0, 100.0)
prob = ODEProblem(lorenz!, u0, tspan)
sol = solve(prob, Tsit5(), abstol = 1e-12, reltol = 1e-12, dt = 1e-6)
# sol = solve(prob, Tsit5(), adaptive = false, dt = 0.01)

function lorenz2!(du, u, p, t)
    du[1] = 10.0 * (u[2] - u[1])
    du[2] = u[1] * (28.0 - u[3]) - u[2]
    du[3] = u[1] * u[2] - (8 / 3) * u[3]
    du[4] = u[1]
end

u02 = [1.0; 0.0; 0.0; 0.0]
prob2 = ODEProblem(lorenz2!, u02, tspan)
sol2 = solve(prob2, Tsit5(), abstol = 1e-12, reltol = 1e-12, dt = 1e-6)
# sol2 = solve(prob2, Tsit5(), adaptive = false, dt = 0.01)

fig = Figure()
ax = Axis(fig[1, 1])
lines!(ax, sol.t, [sol.u[i][1] for i in 1:length(sol.u)])
lines!(ax, sol2.t, [sol2.u[i][1] for i in 1:length(sol2.u)])
fig

adap
adap800×600 142 KB

It gives two radically different solutions. The question is then: can I trust the result of lorenz2, and if not, how could I? This seems to have something to do with the adaptive algorithms, if you try the non-adaptive version (commented lines) you get this:

no-adap
no-adap800×600 117 KB

The two solutions totally overlap, which is what I naively would have expected from the first case. Of course, u'_4 as defined might be changing much faster than any of the other components, so it’s understandable that the adaptive algo might have to compensate for that, but I still don’t understand why u_1 itself wouldn’t be the same. I would also note that the same behavior is observed even if du[4] = 0.0.

2

I’m pretty sure the answer here is that when you solve chaotic systems any error will very quickly put you on a completely uncorrelated trajectory.

2 Likes
3

I’m with Oscar. If you do anything different at all, including even roundoff errors, the chaotic nature of the system will lead to exponential growth of that deviation according to the lyapunov exponent of the system. So if you use a different algorithm in any way you are essentially guaranteed diverging solutions.

4

That’s fair. Integrating 100 steps for this system is somewhat ambitious anyway.

5

The solution lives on a shadow trajectory, which essentially means you get a pretty picture that is on the attractor (and is the exact solution for some value with an epsilon perturbed initial condition) but you never get the real solution in the first place. Any small change will give you a different random shadow trajectory. Chaotic systems are best understood as almost random.