If you fence your code with backticks ‘```’ it will format nicely in your message.
Good to know about the backticks haha
Do you have some reason to use the CVODE_BDF solver?
Actually yes… When I tried using other stiff solvers like QNDF, KenCarp etc, the solution didn’t converges, even for smaller tolerances. However, I would like to use others solvers because CVODE_BDF can’t handle complex-valued ODE, that’s why I always open complex functions in polar coordinate, going from a N-dim complex space to a 2N-dim real space where CVODE_BDF do its job.
The code should be modified to pass these in as part of p
I need to understand better how this works…
There are some constants that you have not defined in this example;
Can you show what you tried and how it failed?
The code is a bit too long, that’s why I haven’t copied it here haha. Nevertheless, here is a runnable example
const ħ = 1.054571817e-34 # [J.s]
const c = 299792458 # [m/s]
const kb = 1.380649e-23 # [J/K]
# Inputs
λ = 1560.0e-9 # bare optical cavity ressonance wavelength [m]
P₀ = 425e-6 # optical power that reaches the cavity (not necessarily, and probably not, the same of the pump laser) [W]
ω₀ = 2*π*c/λ # optical angular frequency [rad/s]
Q₀ = 165000 # optical quality factor
κ = ω₀/Q₀ # optical loss [rad/s]
g₀ = -2*π*16200 # vacuum single-photon optomechanical coupling [rad/s]
η = 0.75 # κe/κ
κₑ = η*κ # external coupling [rad/s]
κ₀ = κ - κₑ # intrinsic optical loss [rad/s]
Ωₘ = 2*π*31.86e6 # natural angular mechanical oscillator frequency [rad/s]
Qₘ = 1255 # mechanical quality factor
Γₘ = Ωₘ/Qₘ # mechanical loss [rad/s]
T = 300 # temperature [K]
nth_a = (kb*T)/(ħ*ω₀) # thermal photons number
nth_b = (kb*T)/(ħ*Ωₘ) # thermal phonons number
# Rescaling parameter
ξ = Ωₘ
# Optical Dynamic Parameters
Δi_x = 8.0*κ # initial optical detuning with bistability correction
Δf_x = 0.35*κ # final optical detuning with bistability correction
#-----------------------------------------------------------------------
# Static Corrections (i.e., Optomechanical Bistability)
function Δi_correction!(F, Δi)
F[1] = Δi[1] - Δi_x + (2*g₀*g₀*κₑ/Ωₘ)*(P₀/(ħ*(Δi[1]+ω₀)))/(κ^2/4 + Δi_x^2)
end
function Δf_correction!(F, Δf)
F[1] = Δf[1] - Δf_x + (2*g₀*g₀*κₑ/Ωₘ)*(P₀/(ħ*(Δf[1]+ω₀)))/(κ^2/4 + Δf_x^2)
end
Δi_guess = Δi_x - (2*g₀*g₀*κₑ/Ωₘ)*(P₀/(ħ*(Δi_x+ω₀)))/(κ^2/4 + Δi_x^2)
Δf_guess = Δf_x - (2*g₀*g₀*κₑ/Ωₘ)*(P₀/(ħ*(Δf_x+ω₀)))/(κ^2/4 + Δf_x^2)
Δi = (nlsolve(Δi_correction!, [Δi_guess], show_trace=false, ftol=1E-10, xtol=1E-12)).zero[1]
Δf = (nlsolve(Δf_correction!, [Δf_guess], show_trace=false, ftol=1E-10, xtol=1E-12)).zero[1]
#-----------------------------------------------------------------------
vΔ = 0.01*κ^2/(2*π) # optical detuning sweep velocity
vΔ = ifelse(Δf-Δi>0, abs(vΔ), -abs(vΔ)) # this corrects the signal of the sweep velocity, if necessary
# RF Parameters
ε = 0.03 # modulation depth
ρ = 2 # winding number, i.e., the harmonic that you are are pumping (for example, ρ = 2 means Ωd around 2Ω₀)
Ω₀_guess = (2*π)*3.2054910037748326e7 # what is, approximatly, the self-sustained frequency of the oscillator?
Ωⁱd = ρ*Ω₀_guess - 4.8*Γₘ # initial modulation frequency
Ωᶠd = ρ*Ω₀_guess + 9.1*Γₘ # final modulation frequency
vΩd = 0.275*Γₘ^2/(2*π) # drive frequency sweep velocity
vΩd = ifelse(Ωᶠd-Ωⁱd>0, abs(vΩd), -abs(vΩd)) # this corrects the signal of the sweep velocity, if necessary
# Time/Delay/Waiting Setup
t_vΔ = (Δf-Δi)/vΔ # time until the optical sweep reaches Δf
t_Δf = 5/(Γₘ/(2*π)) # waiting time after t_vΔ to turn on the modulation depth ε
t_ε = 5/(Γₘ/(2*π)) # waiting time after t_Δf to turn on vΩd
t_vΩd = (Ωᶠd - Ωⁱd)/vΩd # waiting time after t_ε to reach Ωᶠd
t_Ωᶠd = 1/(Γₘ/(2*π)) # waiting time after t_vΩd, here the experiment has already finished
t_final = t_vΔ + t_Δf + t_ε + t_vΩd + t_Ωᶠd
# Optical Detuning function
Δ_func(t) = ifelse(t<=t_vΔ, Δi + vΔ*t, Δf)
# Modulation Depth function
ε_func(t) = ifelse(t<=t_vΔ+t_Δf, 0, ε)
# RF drive function
Θ₁ = 0.5*vΩd*(t_vΔ+t_Δf+t_ε)^2 - vΩd*(t_vΔ+t_Δf+t_ε)*(t_vΔ+t_Δf+t_ε)
Θ₂ = Ωⁱd*(t_vΔ+t_Δf+t_ε+t_vΩd) + 0.5*vΩd*(t_vΔ+t_Δf+t_ε+t_vΩd)^2 - vΩd*(t_vΔ+t_Δf+t_ε)*(t_vΔ+t_Δf+t_ε+t_vΩd) - Θ₁
Θd_func(t) = @. ifelse(t<=t_vΔ+t_Δf+t_ε+t_vΩd, ifelse(t<=t_vΔ+t_Δf+t_ε, Ωⁱd*t, Ωⁱd*t + 0.5*vΩd*t^2 - vΩd*(t_vΔ+t_Δf+t_ε)*t - Θ₁), Ωᶠd*(t - (t_vΔ+t_Δf+t_ε+t_vΩd)) + Θ₂)
Ωd_func(t) = @. ForwardDiff.derivative(Θd_func, t)
fd_func(t) = Ωd_func.(t)/(2*π)
# Optical Pump functions
s₀_func(t) = sqrt(P₀/(ħ*(Δ_func(t)+ω₀)))
s_func(t) = s₀_func(t)*sqrt(1+ε_func(t)*sin(Θd_func(t)));
# Amplitude and Phase Model
#--------------------
function amplitude_phase_optomechanics!(du,u,p,τ)
A,ϕa,B,ϕb = u
du[1] = -(1/2)*(κ/ξ)*A + sqrt(κₑ/ξ)*(s_func(τ/ξ)/sqrt(ξ))*cos(ϕa)
du[2] = (Δ_func(τ/ξ)/ξ) - (2*g₀/ξ)*B*cos(ϕb) - sqrt(κₑ/ξ)*(s_func(τ/ξ)/sqrt(ξ))*(sin(ϕa)/A)
du[3] = -(1/2)*(Γₘ/ξ)*B - (g₀/ξ)*(A^2)*(sin(ϕb))
du[4] = -(Ωₘ/ξ) - (g₀/ξ)*(A^2)*(cos(ϕb)/B)
end
function amplitude_phase_optomechanics_jac(J,u,p,τ)
A,ϕa,B,ϕb = u
J[1,1] = -(1/2)*(κ/ξ)
J[2,1] = sqrt(κₑ/ξ)*(s_func(τ/ξ)/sqrt(ξ))*(sin(ϕa)/(A*A))
J[3,1] = -(2*g₀/ξ)*A*(sin(ϕb))
J[4,1] = -(2*g₀/ξ)*A*(cos(ϕb)/B)
J[1,2] = -sqrt(κₑ/ξ)*(s_func(τ/ξ)/sqrt(ξ))*sin(ϕa)
J[2,2] = -sqrt(κₑ/ξ)*(s_func(τ/ξ)/sqrt(ξ))*(cos(ϕa)/A)
J[3,2] = 0
J[4,2] = 0
J[1,3] = 0
J[2,3] = -(2*g₀/ξ)*cos(ϕb)
J[3,3] = -(1/2)*(Γₘ/ξ)
J[4,3] = (g₀/ξ)*(A^2)*(cos(ϕb)/(B*B))
J[1,4] = 0
J[2,4] = (2*g₀/ξ)*B*sin(ϕb)
J[3,4] = -(g₀/ξ)*(A^2)*(cos(ϕb))
J[4,4] = (g₀/ξ)*(A^2)*(sin(ϕb)/B)
nothing
end
A₀ = sqrt(nth_a)
ϕa₀ = 2*π*rand(1)[1]
B₀ = sqrt(nth_b)
ϕb₀ = 2*π*rand(1)[1]
#--------------------
τ_span = (0, ξ*t_final)
u₀ = [A₀,ϕa₀,B₀,ϕb₀]
ODEfunc = ODEFunction(amplitude_phase_optomechanics!)
ODEprob = ODEProblem(ODEfunc, u₀, τ_span)
# Solver method
sol_abstol = 1e-6
sol_reltol = 1e-6
@time sol = solve(ODEprob, CVODE_BDF(), abstol=sol_abstol, reltol=sol_reltol, maxiters=Int(100e6))
I also would appreciate any hint/good practice that I am not doing in my code so far