Vern9 from differentialequations.jl and ode89 from matlab

Sergey_Novak · May 22, 2024, 5:03pm

Hello, I am calculating a solution for the system using Vern9() in Julia and ode89 in Matlab. I get a message about instability in Julia while there is no such thing in matlab.
If I change the abstol, reltol in julia to 1e-10 and count the trajectory at a time of 250,000 or more, then the trajectory begins to look discrete, this is not observed in matlab.

What could be the reason for this and can it be fixed?
Thanks for the help

Code in Julia:

using StaticArrays, DifferentialEquations

function FHN2_try3(u, p ,t)
    x1, y1, x2, y2, z= u
    ϵ, a, g, k, σ, α, k1, k2 = p

    I(ϕ_i) = g * (1.0/(1.0 + exp(k*(cos(σ/2) - cos(ϕ_i - α - σ/2)))))
    ρz = k1 + k2 * z ^ 2

    ϕ2 = atan(y2, x2)
    ϕ1 = atan(y1, x1)

    dx1dt = (x1 - x1 ^ 3 / 3 - y1 + I(ϕ2) + ρz * (x2 - x1) ) / ϵ
    dy1dt = x1 - a
    dx2dt = (x2 - x2 ^ 3 / 3 - y2 + I(ϕ1) + ρz * (x1 - x2) ) / ϵ
    dy2dt = x2 - a
    dzdt = x1 - x2
    return SVector(dx1dt, dy1dt, dx2dt, dy2dt, dzdt)
end

function FHN2_try3_params()
    ϵ = 0.01; a = -1.01;
    g = 0.1; k = 50.0; σ = 50.0 * pi / 180; α = 160.0 * pi / 180;
    k1 = 0.0; k2 = 0.0
    return [ ϵ, a, g, k, σ, α, k1, k2]
end

parameters = FHN2_try3_params()
tspan = (0.0, 50_000.0)
parameters[7] = 0.09
parameters[8] = 75.74
u0 = [-0.9859005363852416, -0.635253572091177, -1.0345181027025165, -0.636382088705782, 0.0011285166148596525] 
u0 = SVector{5}(u0)
prob = ODEProblem(FHN2_try3, u0, tspan, parameters)

alg = Vern9();
abs_tol = 1e-6;
rel_tol = 1e-3;
max_iters = 1e8;
integrator_setting = (alg = alg, abs_tol = abs_tol, rel_tol = rel_tol, max_iters = max_iters)

sol = solve(prob, integrator_setting.alg, adaptive = true, abstol = integrator_setting.abs_tol, reltol = integrator_setting.rel_tol,
maxiters = integrator_setting.max_iters);

Code in Matlab:

eps = 0.01; a = -1.01;
g = 0.1; k = 50.0; sigma = 50.0 * pi / 180; alpha = 160.0 * pi / 180;
k1 = 0.09; k2 = 75.74;

tspan = [0 50000];
y0 = [-0.9859005363852416, -0.635253572091177, -1.0345181027025165, -0.636382088705782, 0.0011285166148596525];
opts = odeset('Reltol',1e-3,'AbsTol',1e-6);

[t_solve, solve] = ode89(@(t,yvec) fhn2(t, yvec, eps, a, g, k, sigma, alpha, k1, k2), tspan, y0, opts);

plot(t_solve(length(solve)-10000:length(solve)), solve(length(solve)-10000:length(solve), 1))
xlabel('time')
ylabel('y(t)')

function dx = fhn2(~, v, eps, a, g, k, sigma, alpha, k1, k2)
    len = length(v);
    x = v(1:2:len - 1);
    y = v(2:2:len);
    dx = zeros(len, 1);

    dx1 = zeros(len, 1);
    rho = k1 + k2*v(5)^2;
    dx(1:2:len - 1) = (x - x.^3/3 - y + I_vect(atan2(y, x), g, sigma, alpha, k) + rho*[-1, 1; 1, -1]*x)/eps;
    dx(2:2:len) = (x - a);
    dx(end) = (x(1) - x(2));
end

function res = I_vect(phi, g, sigma, alpha ,k)
    G = [0, g;
         g, 0];
    res = G*(1./(1 + exp(k*(cos(sigma/2) - cos(phi - alpha - sigma/2)))));
end

Oscar_Smith · May 22, 2024, 5:34pm

What package versions are you using? I’m not seeing an instability warning. That said I do see the weird jumps that you are mentioning (but only with Vern9 and other Vern integrators, not Tsit5 or Rodas5P(). This does seem to be a bug.

Sergey_Novak · May 22, 2024, 6:02pm

I use version: DifferentialEquations v7.13.0

Do you mean this? If yes, that is not bug, it is extreme events

When i using time for integrate 500 000 i see next in end time series
The trajectory becomes discrete, and over time this discrete becomes more noticeable

For calculate solve for the time 500 000 i use next code:

using StaticArrays, DifferentialEquations, JLD2 

function FHN2_try3(u, p ,t)
    x1, y1, x2, y2, z= u
    ϵ, a, g, k, σ, α, k1, k2 = p

    I(ϕ_i) = g * (1.0/(1.0 + exp(k*(cos(σ/2) - cos(ϕ_i - α - σ/2)))))
    ρz = k1 + k2 * z ^ 2

    ϕ2 = atan(y2, x2)
    ϕ1 = atan(y1, x1)

    dx1dt = (x1 - x1 ^ 3 / 3 - y1 + I(ϕ2) + ρz * (x2 - x1) ) / ϵ
    dy1dt = x1 - a
    dx2dt = (x2 - x2 ^ 3 / 3 - y2 + I(ϕ1) + ρz * (x1 - x2) ) / ϵ
    dy2dt = x2 - a
    dzdt = x1 - x2
    return SVector(dx1dt, dy1dt, dx2dt, dy2dt, dzdt)
end

function FHN2_try3_params()
    ϵ = 0.01; a = -1.01;
    g = 0.1; k = 50.0; σ = 50.0 * pi / 180; α = 160.0 * pi / 180;
    k1 = 0.0; k2 = 0.0
    return [ ϵ, a, g, k, σ, α, k1, k2]
end

function clear_work_space(sol, len_sol, ttr)
    sol = nothing
    len_sol = nothing
    ttr = nothing
    return sol, len_sol, ttr
end

t_truncate(t) = floor(Int64, t / 2)

function save_data(iteration, data, path_to_save)
    file_name = "$(iteration)_sol_x1.jld"
    full_path = path_to_save * file_name
    jldsave(full_path; data)
end

function first_iteration(u0, parameters, tspan, integrator_setting, path_to_save)

    prob = ODEProblem(FHN2_try3, u0, tspan, parameters)
    sol = solve(prob, integrator_setting.alg, adaptive = true,
    abstol = integrator_setting.abs_tol, reltol = integrator_setting.rel_tol,
    maxiters = integrator_setting.max_iters);
    last_point = sol[end]
    println("last point: $(last_point)"); flush(stdout)
    len_sol = length(sol.t)

    ttr = t_truncate(len_sol); tend = len_sol
    data_x1 = [sol[1, ttr:tend], sol.t[ttr:tend]]
    sol, len_sol, ttr = clear_work_space(sol, len_sol, ttr)
    GC.gc()

    save_data(1, data_x1, path_to_save)
    return last_point
end

function calculate_timeseris(u0_start, parameters, integrator_setting,
    t_point, count_iteration, path_to_save)

    global tspan = (0.0, t_point)
    global u0 = u0_start

    u0 = first_iteration(u0, parameters, tspan, integrator_setting, path_to_save)
    println("----------------------------------------------"); flush(stdout)
    println("");flush(stdout)

    for iteration in range(2, count_iteration)

        println("ITERATION: $(iteration)")

        tspan = (t_point*(iteration-1), t_point*iteration)
        prob = ODEProblem(FHN2_try3, u0, tspan, parameters)
        println("u0: $(u0)"); flush(stdout)
        sol = solve(prob, integrator_setting.alg, adaptive = true,
        abstol = integrator_setting.abs_tol, reltol = integrator_setting.rel_tol,
        maxiters = integrator_setting.max_iters);
        u0 = sol[end]

        data_x1 = [sol[1, :], sol.t]
        sol = nothing; GC.gc();
        save_data(iteration, data_x1, path_to_save)

        println("last point: $(u0)"); flush(stdout)
        println("----------------------------------------------"); flush(stdout)
        println("");flush(stdout)
    end
end

alg = Vern9()
abs_tol = 1e-14;
rel_tol = 1e-14;
max_iters = 1e8;
println("alg: $alg"); println("abstol: $abs_tol; reltol: $(rel_tol)")
integrator_setting = (alg = alg, abs_tol = abs_tol, rel_tol = rel_tol, max_iters = max_iters)

path_to_save = "/home/sergey/MEGA/dynamical-systems/FHN_Korotkov/data/timeseries_k2_75_74_change_tol/"
parameters = FHN2_try3_params()
parameters[7] = 0.09
parameters[8] = 75.74

u0_start = [1.7, 0.7, -1.4, 0.35, 0.7 - 0.35]
u0_start = SVector{5}(u0_start)

t_point = 50000.0
count_iteration = 15

calculate_timeseris(u0_start, parameters, integrator_setting,
    t_point, count_iteration, path_to_save)

Sergey_Novak · May 22, 2024, 6:12pm

End of time series for time = 1 000 000
Matlab does not have such a problem and I do not understand what this is related to

End of time series for time = 1 500 000

ChrisRackauckas · May 22, 2024, 7:58pm

Check the values, i.e. no saveat. It could be an instability artifact in the interpolation?

Sergey_Novak · May 23, 2024, 8:21am

Do you mean this?

Code:

using StaticArrays, DifferentialEquations, CairoMakie, GLMakie

function FHN2_try3(u, p ,t)
    x1, y1, x2, y2, z= u
    ϵ, a, g, k, σ, α, k1, k2 = p

    I(ϕ_i) = g * (1.0/(1.0 + exp(k*(cos(σ/2) - cos(ϕ_i - α - σ/2)))))
    ρz = k1 + k2 * z ^ 2

    ϕ2 = atan(y2, x2)
    ϕ1 = atan(y1, x1)

    dx1dt = (x1 - x1 ^ 3 / 3 - y1 + I(ϕ2) + ρz * (x2 - x1) ) / ϵ
    dy1dt = x1 - a
    dx2dt = (x2 - x2 ^ 3 / 3 - y2 + I(ϕ1) + ρz * (x1 - x2) ) / ϵ
    dy2dt = x2 - a
    dzdt = x1 - x2
    return SVector(dx1dt, dy1dt, dx2dt, dy2dt, dzdt)
end

function FHN2_try3_params()
    ϵ = 0.01; a = -1.01;
    g = 0.1; k = 50.0; σ = 50.0 * pi / 180; α = 160.0 * pi / 180;
    k1 = 0.0; k2 = 0.0
    return [ ϵ, a, g, k, σ, α, k1, k2]
end

parameters = FHN2_try3_params()
tspan = (0.0, 200_000.0)
parameters[7] = 0.09
parameters[8] = 75.74
u0 = [-0.9859005363852416, -0.635253572091177, -1.0345181027025165, -0.636382088705782, 0.0011285166148596525] 
u0 = SVector{5}(u0)
prob = ODEProblem(FHN2_try3, u0, tspan, parameters)

alg = Vern9();
abs_tol = 1e-10;
rel_tol = 1e-10;
max_iters = 1e8;
integrator_setting = (alg = alg, abs_tol = abs_tol, rel_tol = rel_tol, max_iters = max_iters)

sol = solve(prob, integrator_setting.alg, adaptive = true, abstol = integrator_setting.abs_tol, reltol = integrator_setting.rel_tol,
maxiters = integrator_setting.max_iters);

tstart = tspan[2] - 100000; tend = tspan[2]
trange = range(tstart, tend, length = 1000000)
values_interp = sol(trange)

f = Figure()
ax = Axis(f[1, 1], xlabel = L"time", ylabel = L"x_1")
lines!(ax, trange, values_interp[1, :], linewidth = 1.0, color = :blue)
display(GLMakie.Screen(), f)

ChrisRackauckas · May 23, 2024, 11:15am

No you’re still interpolating.

sol = solve(prob, integrator_setting.alg, adaptive = true, abstol = integrator_setting.abs_tol, reltol = integrator_setting.rel_tol,
maxiters = integrator_setting.max_iters, dense=false);

and then just plot the values that come out.

Sergey_Novak · May 23, 2024, 11:40am

I disabled interpolation
The result has improved greatly, but there is still a slight discreteness
Thank you

time 200 000

time 500 000

time 1 000 000

Code:

parameters = FHN2_try3_params()
tspan = (0.0, 200_000.0)
parameters[7] = 0.09
parameters[8] = 75.74
u0 = [-0.9859005363852416, -0.635253572091177, -1.0345181027025165, -0.636382088705782, 0.0011285166148596525] 
u0 = SVector{5}(u0)
prob = ODEProblem(FHN2_try3, u0, tspan, parameters)

alg = Vern9();
abs_tol = 1e-10;
rel_tol = 1e-10;
max_iters = 1e8;
integrator_setting = (alg = alg, abs_tol = abs_tol, rel_tol = rel_tol, max_iters = max_iters)

sol = solve(prob, integrator_setting.alg, adaptive = true, abstol = integrator_setting.abs_tol, reltol = integrator_setting.rel_tol,
maxiters = integrator_setting.max_iters, dense = false);

len_sol = length(sol);
tstart = len_sol - 20000;
tend = len_sol;

f = Figure()
ax = Axis(f[1, 1], xlabel = L"time", ylabel = L"x_1")
lines!(ax, sol.t[tstart:tend], sol[1, tstart:tend], linewidth = 1.0, color = :blue)
display(GLMakie.Screen(), f)

ChrisRackauckas · May 23, 2024, 12:18pm

@Oscar_Smith can you look into this? There are very low tolerance tests on the interpolation but we might want to setup more L2 (big L) convergence tests which would check the interpolation is hitting the required 7th order. It could be an instability in the interpolation definition itself, i.e. the math, which is something I have seen before. But we should check this deeply to make sure there’s not something like a 15th digit coefficient that is causing this.

Oscar_Smith · May 23, 2024, 12:54pm

Will investigate.

Sergey_Novak · May 24, 2024, 4:50pm

I think I found the culprit in the fact that the time series looked discrete. That is CairoMakie. I calculated the solution in matlab and upload it to Julia, then plot it using CairoMakie. The time series also looked discrete, although there was a graph from matlab next to it, where everything was continuous. I checked the time diagram calculated using differentialequation.jl by plotting it in Plots and everything looked fine.

abraemer · May 24, 2024, 4:58pm

Then I have good news: Makie.jl just released an updated version 0.21 ([ANN] Makie@0.21 - #4 by barucden) that should fix that just by updating

Sergey_Novak · May 24, 2024, 6:01pm

I tried it. CairoMakie now does not display the image in vs code. I used CairoMakie.activate!() and Makie.inline!(). It did not help

When i open terminal and using Julia REPL in it CairoMakie show image

jules · May 24, 2024, 7:32pm

For me that works, maybe your plot pane might be disabled for some reason?

rveltz · May 24, 2024, 9:59pm

Makie.inline!(true)

Sergey_Novak · May 25, 2024, 8:00am

The picture was not displayed for a very strange reason. I already had the REPL running, where the calculation was going on, and I opened another one in the terminal via vs-code and try plot in it. When the calculating ran out in the first REPL and I started again REPL, everything worked.