# Vern9 from differentialequations.jl and ode89 from matlab
**URL:** https://discourse.julialang.org/t/vern9-from-differentialequations-jl-and-ode89-from-matlab/114587
**Category:** Modelling & Simulations
**Created:** [May 22, 2024, 5:03pm UTC](https://discourse.julialang.org/t/vern9-from-differentialequations-jl-and-ode89-from-matlab/114587 "2024-05-22T17:03:08Z")
**Posts on this page:** 16
**Page:** 1
### Author:  [@Sergey\_Novak](https://discourse.julialang.org/u/Sergey_Novak)
#### Post date: [May 22, 2024, 5:03pm UTC](https://discourse.julialang.org/t/vern9-from-differentialequations-jl-and-ode89-from-matlab/114587/1 "2024-05-22T17:03:08Z")
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:
```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:
```julia
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
```
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [May 22, 2024, 5:34pm UTC](https://discourse.julialang.org/t/vern9-from-differentialequations-jl-and-ode89-from-matlab/114587/2 "2024-05-22T17:34:50Z")
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.
---
### Author:  [@Sergey\_Novak](https://discourse.julialang.org/u/Sergey_Novak)
#### Post date: [May 22, 2024, 6:02pm UTC](https://discourse.julialang.org/t/vern9-from-differentialequations-jl-and-ode89-from-matlab/114587/3 "2024-05-22T18:02:12Z")
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:
```julia
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)
```
---
### Author:  [@Sergey\_Novak](https://discourse.julialang.org/u/Sergey_Novak)
#### Post date: [May 22, 2024, 6:12pm UTC](https://discourse.julialang.org/t/vern9-from-differentialequations-jl-and-ode89-from-matlab/114587/4 "2024-05-22T18:12:06Z")
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
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 22, 2024, 7:58pm UTC](https://discourse.julialang.org/t/vern9-from-differentialequations-jl-and-ode89-from-matlab/114587/5 "2024-05-22T19:58:55Z")
Check the values, i.e. no saveat. It could be an instability artifact in the interpolation?
---
### Author:  [@Sergey\_Novak](https://discourse.julialang.org/u/Sergey_Novak)
#### Post date: [May 23, 2024, 8:21am UTC](https://discourse.julialang.org/t/vern9-from-differentialequations-jl-and-ode89-from-matlab/114587/6 "2024-05-23T08:21:54Z")
Do you mean this?
Code:
```julia
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)
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 23, 2024, 11:15am UTC](https://discourse.julialang.org/t/vern9-from-differentialequations-jl-and-ode89-from-matlab/114587/7 "2024-05-23T11:15:37Z")
> [@Sergey\_Novak](#):
>
> `values_interp = sol(trange)`
No you’re still interpolating.
```julia
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.
---
### Author:  [@Sergey\_Novak](https://discourse.julialang.org/u/Sergey_Novak)
#### Post date: [May 23, 2024, 11:40am UTC](https://discourse.julialang.org/t/vern9-from-differentialequations-jl-and-ode89-from-matlab/114587/8 "2024-05-23T11:40:22Z")
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:
```julia
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)
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 23, 2024, 12:18pm UTC](https://discourse.julialang.org/t/vern9-from-differentialequations-jl-and-ode89-from-matlab/114587/9 "2024-05-23T12:18:44Z")
@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.
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [May 23, 2024, 12:54pm UTC](https://discourse.julialang.org/t/vern9-from-differentialequations-jl-and-ode89-from-matlab/114587/10 "2024-05-23T12:54:08Z")
Will investigate.
---
### Author:  [@Sergey\_Novak](https://discourse.julialang.org/u/Sergey_Novak)
#### Post date: [May 24, 2024, 4:50pm UTC](https://discourse.julialang.org/t/vern9-from-differentialequations-jl-and-ode89-from-matlab/114587/11 "2024-05-24T16:50:35Z")
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.
---
### Author:  [@abraemer](https://discourse.julialang.org/u/abraemer)
#### Post date: [May 24, 2024, 4:58pm UTC](https://discourse.julialang.org/t/vern9-from-differentialequations-jl-and-ode89-from-matlab/114587/12 "2024-05-24T16:58:45Z")
Then I have good news: Makie.jl just released an updated version 0.21 ([[ANN] Makie@0.21 - #4 by barucden](https://discourse.julialang.org/t/ann-makie-0-21/114678/4)) that should fix that just by updating 🙂
---
### Author:  [@Sergey\_Novak](https://discourse.julialang.org/u/Sergey_Novak)
#### Post date: [May 24, 2024, 6:01pm UTC](https://discourse.julialang.org/t/vern9-from-differentialequations-jl-and-ode89-from-matlab/114587/13 "2024-05-24T18:01:41Z")
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
---
### Author:  [@jules](https://discourse.julialang.org/u/jules)
#### Post date: [May 24, 2024, 7:32pm UTC](https://discourse.julialang.org/t/vern9-from-differentialequations-jl-and-ode89-from-matlab/114587/14 "2024-05-24T19:32:08Z")
For me that works, maybe your plot pane might be disabled for some reason?
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [May 24, 2024, 9:59pm UTC](https://discourse.julialang.org/t/vern9-from-differentialequations-jl-and-ode89-from-matlab/114587/15 "2024-05-24T21:59:37Z")
> [@Sergey\_Novak](#):
>
> Makie.inline!()
```julia
Makie.inline!(true)
```
---
### Author:  [@Sergey\_Novak](https://discourse.julialang.org/u/Sergey_Novak)
#### Post date: [May 25, 2024, 8:00am UTC](https://discourse.julialang.org/t/vern9-from-differentialequations-jl-and-ode89-from-matlab/114587/16 "2024-05-25T08:00:17Z")
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.