Homoclinic loop accuracy DifferentialEquations.jl

Sergey_Novak · April 14, 2023, 7:10am

Hello!
I try get homoclinic loop to saddle focus for system of Rossler. With the given parameters, the homoclinic definitely exists, according to some articles. To get homoclinic loop, I integrate in backward time. When I do this I get an instability message. Why I get such a message I understand. I suppose the problem is with accuracy. In the same articles, people managed to get a homoclinic loop using the 4-order Runge Kutta. I tried to use Verner9, Feigin12. I used big float and it didn’t help. I also made the time step quite small. What could be the problem? Is it all about accuracy, or am I doing something wrong even before the integration?

Code:

using StaticArrays, DifferentialEquations, DynamicalSystems

x, y, z  = -10..10, -10..10,  -10..10
box = x × y × z

using ForwardDiff
using LinearAlgebra

@inbounds function Res(u, p, t)
    du1 = -u[2]-u[3]
    du2 = u[1]+p[1]*u[2]
    du3 = p[2]*u[1]-p[3]*u[3]+u[1]*u[3]
    return SVector(du1, du2, du3)
end
@inbounds function jac_res(u, p, t)
    SMatrix{3, 3}(0.0, 1.0, p[2]+u[3],
                -1.0, p[1], 0.0,
                -1.0, 0.0, -p[3]+u[1])
end

time_set = (t = 500, Ttr = 250, tstep = 0.001)
integ_set = (alg = Vern9(), adaptive = false, dt = time_set.tstep)

u0 = SA[0.13060105886751308, 0.4782615699018184, 0.7435866656901606]
b = 0.3;  c = 4.89; a = 0.35;
p = [a, b, c]
ds = CoupledODEs(Res, u0, p, diffeq = integ_set)

fp, ei, _ = fixedpoints(ds, box, jac_res)

Jacobian = jac_res(fp[1], p, 0)
eivecs = eigvecs(Jacobian)
vec_stable = real(eivecs[:, 1])
ϵ = 1e-23
shift = SA[0.0, 0.0, 0.0] + vec_stable*ϵ

Jacobi(u0, p, t) = ForwardDiff.jacobian((x) -> Res(x, p, 0), u0)

prob = ODEProblem(Res, shift, (0.0, -1000.0),p)

sol = solve(prob,alg = RK4(), adaptive = false, dt = 0.001);

cmarcotte · April 14, 2023, 8:52am

A couple things:
You didn’t define your a,b,c so we can’t run your code as you have (though I just used a=0.1,b=0.1,c=14.0).

Additionally, if you plot [my] solution, the divergence begins in earnest around t=-16.0, just before which the z-component becomes negative; this is consistent for adaptive stepping and fixed (though the precise time changes). This seems like a problem with solving the system in backward time – \dot{z} < 0 if z=0 and p_2 u_1 > 0, which is eminently possible. In fact, using Roots.jl, we can find the precise times that z=0, t0 = find_roots((t)->sol(t)[3], tspan), and then check whether the (negative) flow has a negative component along z.

Sergey_Novak · April 14, 2023, 9:03am

Sorry, I added parameter values
I don’t quite understand what you mean
I think the problem is the accuracy. If I did everything correctly, then the starting point corresponds to the displacement along the stable manifold from the saddle focus.

empet · April 14, 2023, 9:34am

I recommend reading this recent paper by @ranocha: https://arxiv.org/abs/2304.02365, to understand why you are surprised by the computed homoclinic orbit.

empet · April 14, 2023, 9:56am

Running your code I noticed that you computed the eigenvalues and eigenvectors of Jacobian matrix evaluated at u0, but u0 is not an equilibrium point for your system. Or you are searching for a homoclinic orbit to a particular equilibrium.
To check it I defined the vector field, having as components the RHS of your system, ie.:

F(u::SVector{3,Float64}; p = [0.35, 0.3, 4.89])=(-u[2]-u[3], u[1]+p[1]*u[2], p[2]*u[1]-p[3]*u[3]+u[1]*u[3])

F evaluated at u0 gives:

F(u0)
(-1.221848235591979, 0.2979926083331495, -3.4998452716657327)

not (0,0,0).

Sergey_Novak · April 14, 2023, 10:30am

Why u0?
Eigenvalues and eigenvectors calculate for fp[1]
Fragment code:

fp, ei, _ = fixedpoints(ds, box, jac_res)

Jacobian = jac_res(fp[1], p, 0)
eivecs = eigvecs(Jacobian)
vec_stable = real(eivecs[:, 1])
ϵ = 1e-23
shift = SA[0.0, 0.0, 0.0] + vec_stable*ϵ

Thank you for link on article

ChrisRackauckas · April 14, 2023, 10:42am

Well I mean, you turned off the accuracy handling. Don’t use adaptive=false if you care about solving to tolerance or high accuracy. Set the abstol and reltol to the desired local tolerance and let the integrator solve to tolerance.

Sergey_Novak · April 14, 2023, 10:58am

With next tol i get message about unstable solution

sol = solve(prob,alg = RK4(), abstol=1e-13, reltol=1e-13, maxiters = 10000000)

dt(-1.1368683772161603e-13) <= dtmin(1.1368683772161603e-13) at t=-13.637725296644767. Aborting. There is either an error in your model specification or the true solution is unstable.

ChrisRackauckas · April 14, 2023, 11:03am

That means RK4 cannot hit that accuracy within the chosen dt, which isn’t too surprising because it’s low order. What about Vern9?

Sergey_Novak · April 14, 2023, 11:05am

I tried the Vern9. It also displays a message about an unstable solution

cmarcotte · April 14, 2023, 11:08am

One issue is it’s not vanilla Rössler:
\dot{x} = -y-z, \dot{y} = x + a y, and \dot{z} = b + z(x-c), but instead
\dot{x} = -y -z, \dot{y} = x + a y, and \dot{z} = b x + z(x-c).

Changing it to Rossler original edition, you immediately blow up (even with Feagin14() and abstol=reltol=1e-13. The issue is that with your initial condition at (approximately) z=0 immediately encounters \dot{z} < 0 (in backward time), and after that, \dot{z} \propto z, so the solution blows up. I suspect the same happens with your system – it steps over z=0 because the p_2u_1 term is negative at that time for the backwards flow. If you use a ContinuousCallback to ensure that z\geq 0, then perhaps it will be more reliable. I’m not convinced this isn’t a mathematical issue, at heart, rather than a numerical one. You’re essentially hoping that the attractor is stable and bounded under the reverse flow, but thats not obvious to me.

So my recommendations are:

change the flow to Rössler, if that’s what you actually want,
Try a different initial condition with low tolerances and high order and adaptivity to see if the flow is actually bounded in backward time,
implement a ContinuousCallback to keep z \geq 0 (or z > 0, probably).

Sergey_Novak · April 14, 2023, 11:24am

The system for which I actually want to get a loop is different. I took the Rössler system to get a loop in a simpler system. How i can correct use Callback? More precisely, what should he do after checking for the condition?

empet · April 14, 2023, 11:33am

You’re right, I don’t know how I read in a hurry on the REPL

cmarcotte · April 14, 2023, 11:43am

For the callback I can’t tell you what should be done – you have the ability to modify the state values to prevent z from becoming negative, but what particular form that takes is up to you and how you wish to perturb* your system.

*state modification is not guaranteed to be ‘perturbative’ in the mathematical sense.

Sergey_Novak · April 28, 2023, 12:30pm

I am experimenting with the Feagin 12. Is it possible to use type of data float64 for this method?
With next code i don’t get message with warnings about maxiter.

prob_for = ODEProblem(TM, shift, (0.0, 500), p)
sol_for = solve(prob_for, alg = Feagin12(), abstol = 1e-18, reltol = 1e-18);

I tried use Feagin14 and i get message about Instability, when using abstol = 1e-16, reltol = 1e-16

ChrisRackauckas · April 30, 2023, 2:23pm

Yes, as always it’s just done by how you choose your initial values.

julia> eps(Float64)
2.220446049250313e-16

so a relative tolerance that low is generally impossible with Float64. Anything in Float64 range I would recommend Vern9 anyways.

Sergey_Novak · July 18, 2023, 10:03am

I decided to try implicit methods to solve stiff problems. With them, the maxiters option is not required, and there is no need to calculate the trajectory for a long time to get an instability message.

The radau method does not work and the following error appears. Did I forget any options for this method?

sol = solve(prob, ODEInterfaceDiffEq.radau())

MethodError: no method matching radau(::ODEInterfaceDiffEq.var"#12#18"{ODEProblem{SVector{3, Float64}, Tuple{Float64, Float64}, false, SVector{11, Float64}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(TM), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tuple{Int64}}, ::Float64, ::Float64, ::SVector{3, Float64}, ::ODEInterface.OptionsODE)

Closest candidates are:
  radau(::Any, ::Real, ::Real, !Matched::Vector, ::ODEInterface.AbstractOptionsODE)
   @ ODEInterface C:\Users\Alex\.julia\packages\ODEInterface\RwRLn\src\Radau.jl:742


Stacktrace:
  [1] __solve(prob::ODEProblem{SVector{3, Float64}, Tuple{Float64, Float64}, false, SVector{11, Float64}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(TM), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::radau{Nothing}, timeseries::Vector{Any}, ts::Vector{Any}, ks::Vector{Any}; saveat::Vector{Float64}, verbose::Bool, save_everystep::Bool, save_on::Bool, save_start::Bool, timeseries_errors::Bool, dense_errors::Bool, callback::Nothing, alias_u0::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ODEInterfaceDiffEq C:\Users\Alex\.julia\packages\ODEInterfaceDiffEq\lGV1B\src\solve.jl:144
  [2] __solve(prob::ODEProblem{SVector{3, Float64}, Tuple{Float64, Float64}, false, SVector{11, Float64}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(TM), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::radau{Nothing}, timeseries::Vector{Any}, ts::Vector{Any}, ks::Vector{Any})
    @ ODEInterfaceDiffEq C:\Users\Alex\.julia\packages\ODEInterfaceDiffEq\lGV1B\src\solve.jl:1
  [3] __solve
    @ C:\Users\Alex\.julia\packages\ODEInterfaceDiffEq\lGV1B\src\solve.jl:1 [inlined]
  [4] #solve_call#22
    @ C:\Users\Alex\.julia\packages\DiffEqBase\G15op\src\solve.jl:511 [inlined]
  [5] solve_call
    @ C:\Users\Alex\.julia\packages\DiffEqBase\G15op\src\solve.jl:481 [inlined]
  [6] #solve_up#30
    @ C:\Users\Alex\.julia\packages\DiffEqBase\G15op\src\solve.jl:972 [inlined]
  [7] solve_up
    @ C:\Users\Alex\.julia\packages\DiffEqBase\G15op\src\solve.jl:945 [inlined]
  [8] #solve#28
    @ C:\Users\Alex\.julia\packages\DiffEqBase\G15op\src\solve.jl:882 [inlined]
  [9] solve(prob::ODEProblem{SVector{3, Float64}, Tuple{Float64, Float64}, false, SVector{11, Float64}, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(TM), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::radau{Nothing})
    @ DiffEqBase C:\Users\Alex\.julia\packages\DiffEqBase\G15op\src\solve.jl:872
 [10] top-level scope
    @ c:\Users\Alex\Desktop\dynamical-systems\Tsodyks Markram\

ChrisRackauckas · July 18, 2023, 10:52am

Make an mwe

Sergey_Novak · July 18, 2023, 11:04am

What is you mean by mwe?

John_Gibson · July 18, 2023, 11:47am

MWE == Minimum Working Example, i.e. a small, self-contained piece of code that exhibits the problem. You should be able to copy-paste a MWE into the REPL, run it, and trigger the error.