Calculating equilibrium states takes too long IntervalRootFinding

Specific Domains Modelling & Simulations
1

Hello, i calculate equilibriums with IntervalRootFinding.Newton
The calculation time takes too long
Could it be some kind of error or a problem in the method itself?

I call the method of Newton using fixedpoints from ChaosTools.jl

Code:

using StaticArrays, DifferentialEquations, DynamicalSystems, LinearAlgebra, IntervalRootFinding 
@inbounds function TM(u, p, t)
    U(y, p) = p[8] + p[9] / ( 1.0 + exp( -50.0 * (y - p[7]) ) )
    σ(x, p) = 1.0 / ( 1.0 + exp( -20.0 * (x-p[6]) ) )
    g(E, x, y, p, U_) = log( 1.0 + exp( (p[5] * U_ * x * E + p[11]  ) / (p[1]) ) )
    
    U_ = U(u[3], p)
    du1 = (-u[1] + p[1] * g(u[1], u[2], u[3], p, U_) ) / p[2]
    du2 = (1.0 - u[2]) / p[3] - U_*u[2]*u[1]
    du3 = (-u[3])/p[4] + p[10] * σ(u[2], p)
    
    return SVector(du1, du2, du3)
end

@inbounds function jacob_TM_(u, p, t)
    
    U(y, p, exp50) = p[8] + p[9] / ( 1.0 + exp50 )
    U_y(y, p, exp50) = (50.0 * p[9] * exp50) / (1.0 + exp50)^2
    g(E, x, y, p, U_) = exp((p[5]  * U_ * x * E + p[11]) / p[1])
    σ_der(x, p) = exp( (-20.0) * (x - p[6]) )
    exp50 = exp(-50.0 * (u[3] - p[7]))
    
    U_ = U(u[3], p, exp50)
    Uy = U_y(u[3], p, exp50)
    g_ = g(u[1], u[2], u[3], p, U_)
    σ_deri = σ_der(u[2], p)
    
    g_plus = 1.0 + g_
    g_mult = g_ * U_
    g_plus_mult = p[2] * (g_plus)
    u1p5 = p[5] * u[1]
    Uyu2 = Uy * u[2]
    
    E_E = (-1.0 + ((J * u[2] * g_mult)) / (g_plus) ) / p[2]
    E_x = (u1p5 * g_mult) / (g_plus_mult)
    E_y = (u1p5 * Uyu2 * g_) / (g_plus_mult)
    
    x_E = -U_ * u[2]
    x_x = -1.0 / p[3] - U_ * u[1]
    x_y = -Uyu2 * u[1]
    
    y_x = 20.0 * p[10] * σ_deri / (1.0 + σ_deri)^2
    y_y = -1.0/p[4]
    
    SMatrix{3,3}(E_E, x_E, 0.0,
        E_x, x_x, y_x,
        E_y, x_y, y_y)
end

t = 500.0; tt = 600.0; tstep = 0.001
trange = range(0.0, t, step = tstep)
integ_set = (alg = RK4(), adaptive = false, dt = tstep)

E = interval(-40, 40); x = interval(0, 1); y = interval(0, 1)
box = IntervalBox(E, x, y);

const τ = 0.013;  const τD = 0.07993;  const τy = 3.3;  const J = 3.07;  const β = 0.300
const xthr = 0.75; const ythr = 0.4
const α = 1.58; const ΔU0 = 0.305;

U0, I0 = 0.182490142150041, -1.33122064887577
p = SA[α, τ, τD, τy, J, xthr, ythr, U0, ΔU0, β, I0];

u0 = [7.485044731286254, 0.771932479798953, 0.5235934863993938]
ds = CoupledODEs(TM, u0, p, diffeq = integ_set)

fp, _, _ = fixedpoints(ds, box, jacob_TM_, method = IntervalRootFinding.Newton)
2

Could you explain what “too long” means? why are you expecting the calculation to be shorter, i.e. are you comparing to an implementation in another programming language, does it scale badly from small to large inputs, or do you have some theoretical complexity calculations that tell you what the expected time for an efficient algorithm here should be?

3

With other values parameters of the system , the calculation takes about 10 seconds.
But at the moment the bill lasts more than an hour.

4

what’s the trajectorties of the system at these parametrers where it takes too long?

5

Here I have built three trajectories. Green dots marked the initial conditions. It can be seen from the phase space that they all go to a stable state of equilibrium.

image
image800×600 29.7 KB

6

Well if the code works for other system parameters then the code is fine and it is the system that is the problem.

For some reason at this parameters the newton’s method fails to converge. Unfortuantely the IntervalRootsFinding doesn’t provide any API to tweak the method. Try using a different solver instead of Newton: Krawczyk or Bisection

7

I tried the Krawczyk method, it didn’t help either. I don’t think the problem is that Newton’s method doesn’t converge, because the parameters correspond to the zero hopf bifurcation, which was find matcont

In addition, the bifurcationkit finds an equilibrium state with these parameters. It is not clear why there are problems with finding all the equilibrium states in the interval