Hello, I got this error after two hours of calculation
ArgumentError: The interval [a,b] is not a bracketing interval.
You need f(a) and f(b) to have different signs (f(a) * f(b) < 0).
Consider a different bracket or try fzero(f, c) with an initial guess c.
Please help me with fix this error. I don’t understand the reason why it occurred. Code:
using DifferentialEquations, DynamicalSystems, StaticArrays, GLMakie, 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 plot_bifurcation_diagram(output, rangek1)
markersize = 0.5;
lbsize = 30;
ticksize = 30;
fig = Figure(size = (1200, 350))
axis = Axis(fig[1,1],
xlabel = L"k_1", ylabel = L"x_1",
xlabelsize = lbsize, ylabelsize = lbsize,
xticklabelsize = ticksize,yticklabelsize = ticksize)
for (j, p) in enumerate(rangek1)
scatter!(axis, fill(p, length(output[j])), output[j];
color = ("black", 0.5), markersize = markersize)
end
display(GLMakie.Screen(), fig)
end
function bifurcation_diagram()
#u0 = [1.0, 0.0, 0.01, -1.0, 0.0]
u0 = [-1.0, 0.0, 0.01, -1.0, 0.0]
params = FHN2_try3_params()
integ_set = (alg = RK4(), adaptive = false, dt = 0.001)
ds = CoupledODEs(FHN2_try3, u0, params, diffeq = integ_set)
t = 1000
ttr = 250
k1_start = 0.0
k1_end = 0.094
len = 5000
rangek1 = range(k1_start, k1_end, length = len)
index_control_param = 7
index_saving_var = 1
index_fixed_var = 4
value_fixed_var = 0.0
surface = (index_fixed_var, value_fixed_var)
setting_root = (xrtol = 1e-11, atol = 1e-11)
pmap = PoincareMap(ds, surface, rootkw = setting_root)
output = orbitdiagram(pmap, index_saving_var, index_control_param, rangek1;
n = t, Ttr = ttr, show_progress = true)
return output
end
k1_start = 0.0
k1_end = 0.094
len = 5000
rangek1 = range(k1_start, k1_end, length = len)
output = bifurcation_diagram()
plot_bifurcation_diagram(output, rangek1)
Stacktrace:
[1] assert_bracket
@ ~/.julia/packages/Roots/BMiNe/src/Bracketing/bracketing.jl:52 [inlined]
[2] init_state(::Roots.A42, F::Roots.Callable_Function{…}, x₀::Float64, x₁::Float64, fx₀::Float64, fx₁::Float64; c::Nothing, fc::Nothing)
@ Roots ~/.julia/packages/Roots/BMiNe/src/Bracketing/alefeld_potra_shi.jl:59
[3] init_state
@ ~/.julia/packages/Roots/BMiNe/src/Bracketing/alefeld_potra_shi.jl:49 [inlined]
[4] init_state(M::Roots.A42, F::Roots.Callable_Function{…}, x::Tuple{…})
@ Roots ~/.julia/packages/Roots/BMiNe/src/Bracketing/bracketing.jl:6
[5] #init#42
@ ~/.julia/packages/Roots/BMiNe/src/find_zero.jl:299 [inlined]
[6] init
@ ~/.julia/packages/Roots/BMiNe/src/find_zero.jl:289 [inlined]
[7] solve(𝑭𝑿::Roots.ZeroProblem{…}, M::Roots.A42, p::Nothing; verbose::Bool, kwargs::@Kwargs{…})
@ Roots ~/.julia/packages/Roots/BMiNe/src/find_zero.jl:491
[8] find_zero(f::Function, x0::Tuple{…}, M::Roots.A42, p′::Nothing; p::Nothing, verbose::Bool, tracks::Roots.NullTracks, kwargs::@Kwargs{…})
@ Roots ~/.julia/packages/Roots/BMiNe/src/find_zero.jl:220
[9] find_zero (repeats 2 times)
@ ~/.julia/packages/Roots/BMiNe/src/find_zero.jl:210 [inlined]
[10] poincare_step!(ds::CoupledODEs{…}, plane_distance::Function, planecrossing::PlaneCrossing{…}, Tmax::Float64, rootkw::@NamedTuple{…})
@ DynamicalSystemsBase ~/.julia/packages/DynamicalSystemsBase/N5A3S/src/derived_systems/poincare/poincaremap.jl:238
[11] step!(pmap::PoincareMap{…})
@ DynamicalSystemsBase ~/.julia/packages/DynamicalSystemsBase/N5A3S/src/derived_systems/poincare/poincaremap.jl:155
[12] step!(pmap::PoincareMap{…}) (repeats 2 times)
@ DynamicalSystemsBase ~/.julia/packages/DynamicalSystemsBase/N5A3S/src/derived_systems/poincare/poincaremap.jl:165 [inlined]
[13] fill_orbitdiagram!(output::Vector{…}, ds::PoincareMap{…}, i::Int64, pvalues::StepRangeLen{…}, p_index::Int64; n::Int64, Ttr::Int64, u0::Nothing, Δt::Int64, ulims::Nothing, show_progress::Bool, periods::Nothing)
@ ChaosTools ~/.julia/packages/ChaosTools/49ZtY/src/timeevolution/orbitdiagram.jl:83
[14] fill_orbitdiagram!
@ ChaosTools ~/.julia/packages/ChaosTools/49ZtY/src/timeevolution/orbitdiagram.jl:56 [inlined]
[15] #orbitdiagram#1
@ ChaosTools ~/.julia/packages/ChaosTools/49ZtY/src/timeevolution/orbitdiagram.jl:44 [inlined]
[16] bifurcation_diagram()
@ Main ~/work/repo/dynamical-systems/FHN_Korotkov/poincare_diagram.jl:81
[17] top-level scope
@ ~/work/repo/dynamical-systems/FHN_Korotkov/poincare_diagram.jl:92
Some type information was truncated. Use `show(err)` to see complete types.
Thank you for your help