Thank you
1 Like
Code:
@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());
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?
I’m sorry, I reduced the code
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())
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
@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];
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