# Homoclinic loop accuracy DifferentialEquations.jl

**URL:** <https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456>\
**Category:** Modelling & Simulations\
**Created:** [April 14, 2023, 7:10am UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456 "2023-04-14T07:10:59Z")\
**Posts on this page:** 6\
**Page:** 2

<div class="post-metadata">

**Author:** ![Sergey\_Novak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sergey_novak/32/37716_2.png) [@Sergey\_Novak](https://discourse.julialang.org/u/Sergey_Novak)\
**Post date:** [July 18, 2023, 11:57am UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/21 "2023-07-18T11:57:25Z")

</div>

Thank you

---

<div class="post-metadata">

**Author:** ![Sergey\_Novak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sergey_novak/32/37716_2.png) [@Sergey\_Novak](https://discourse.julialang.org/u/Sergey_Novak)\
**Post date:** [July 18, 2023, 12:00pm UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/22 "2023-07-18T12:00:18Z")

</div>

Code:

```julia
@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

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

I0 = -1.706317723300771; U0 = 0.26500158941163476
p = SA[α, τ, τD, τy, J, xthr, ythr, U0, ΔU0, β, I0]
u0 = [11.325905642223786, 0.6594706953104683, 0.4863175026548461];

ds = CoupledODEs(TM, u0, p)
fp, _, _ = fixedpoints(ds, box, jacob_TM_)

ϵ_sep = 1e-6

Ju0 = jacob_TM_(fp[1], p, 0);
eigen_val_vec = eigen(Ju0)
eigen_vectors = eigen_val_vec.vectors

v1 = real(eigen_vectors[:, 1])

dot_on_sep_stable = fp[1] + v1*ϵ_sep
t = 1000.0; tspan = [0.0, t]
prob = ODEProblem(TM, dot_on_sep_stable, tspan, p)
sol = solve(prob, ODEInterfaceDiffEq.radau());

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 18, 2023, 1:50pm UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/23 "2023-07-18T13:50:33Z")

</div>

I copy pasted it and it said `CoupledODEs` is not defined. Please post a runnable MWE. Also, “minimal” seems to not be taken to heart here. Are you sure you need `CoupledODEs` for this error to occur, or are you getting it with only the ODE?

---

<div class="post-metadata">

**Author:** ![Sergey\_Novak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sergey_novak/32/37716_2.png) [@Sergey\_Novak](https://discourse.julialang.org/u/Sergey_Novak)\
**Post date:** [July 18, 2023, 2:04pm UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/24 "2023-07-18T14:04:47Z")

</div>

I’m sorry, I reduced the code

```julia
using StaticArrays, DifferentialEquations, ODEInterfaceDiffEq

@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

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

I0 = -1.706317723300771; U0 = 0.26500158941163476
p = SA[α, τ, τD, τy, J, xthr, ythr, U0, ΔU0, β, I0]
u0 = SA[8.345789408165038,0.7384945878731853, 0.43829843650314865];

t = 1000.0; tspan = [0.0, t]
prob = ODEProblem(TM, u0, tspan, p)
sol = solve(prob, ODEInterfaceDiffEq.radau())

```

---

<div class="post-metadata">

**Author:** ![Sergey\_Novak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sergey_novak/32/37716_2.png) [@Sergey\_Novak](https://discourse.julialang.org/u/Sergey_Novak)\
**Post date:** [July 19, 2023, 1:58pm UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/25 "2023-07-19T13:58:25Z")

</div>

I managed to fix the error. I have removed static arrays and vectors. Most likely, I haven’t read the documentation carefully enough, sorry  
The changes were in the following code snippets

```julia
@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 [du1, du2, du3]
end

p = [α, τ, τD, τy, J, xthr, ythr, U0, ΔU0, β, I0]
u0 = [8.345789408165038,0.7384945878731853, 0.43829843650314865];

```

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [July 19, 2023, 2:19pm UTC](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456/26 "2023-07-19T14:19:24Z")

</div>

This isn’t really related to your question but your `U` and `σ` should probably be using `tanh` since it has some better numerical properties. Also `g` is a `softplus` function which isn’t in Base, but I think we have in LogExpFunctions.jl

[Previous page](https://discourse.julialang.org/t/homoclinic-loop-accuracy-differentialequations-jl/97456.md?page=1)