Hello everyone, I want about the most suitable ODE solver for the following problem:
using DifferentialEquations
using Plots; plotly()
#Parameters
omega_m = 1.6e6
N = 1
Delta_mw = -0.7
g0 = 3.6e3#/omega_m
g_eff = (g0/sqrt(N))
Gamma = 60000#/omega_m
Gamma_tilde = 0.0
Omega_mw = 0.7
Q = 2*(10^5)
gamma_m = (omega_m/Q) #/ omega_m
N_th = 6.05
Angle_th = 0;
tau = (Gamma/(g_eff^2));
#omega_m = 1
function mean_field(du,u,p,t)
#p[1]=Delta_mw
#p[2]=Omega_mw
Sx =u[1]
Sy =u[2]
Sz =u[3]
alpha_x=u[4]
alpha_y=u[5]
du[1] = (Delta_mw-4*g_eff*u[4])*u[2] - (Gamma+Gamma_tilde)*u[1]
du[2] = -(Delta_mw-4*g_eff*u[4])*u[1] -Omega_mw*u[3] - (Gamma+Gamma_tilde)*u[2]
du[3] = Omega_mw*u[2] - 2*Gamma*(u[3]+N)
du[4] = omega_m*u[5] - (gamma_m /2)*u[4]
du[5] = -omega_m*u[4] -g_eff*u[3] - (gamma_m /2)*u[5]
#du[1] = (p[1]-4*g_eff*alpha_x)*Sy - (Gamma+Gamma_tilde)*Sx
#du[2] = -(p[1]-4*g_eff*alpha_x)*Sx -p[2]*Sz - (Gamma+Gamma_tilde)*Sy
#du[3] = p[2]*Sy - 2*Gamma*(Sz+N)
#du[4] = omega_m*alpha_y - (gamma_m /2)*alpha_x
#du[5] = -omega_m*alpha_x - g_eff*Sz - (gamma_m /2)*alpha_y
end
function Result_lim(Array_sol)
print((Array_sol[2][4])^2+(Array_sol[2][5])^2," ")
end
#Initial state
u0 = [0;0;-N;sqrt(N_th);0];
time= 100*tau
function Quanta_fin(Delta_mw,Omega_mw)
parameters=[Delta_mw ;Omega_mw ]
tspan = (0.0,time)
prob = ODEProblem(mean_field,u0,tspan,parameters) #,parameters)
Val1=1e-6
Val2=1e-6
Iter=1e14
dt_size=1e-7
sol1 = solve(prob,RK4(),dt=dt_size,save_everystep=false,reltol=Val1,abstol=Val2,MaxIters=Iter)
sol2 = solve(prob,alg_hints=[:stiff],dt=dt_size,reltol=Val1,abstol=Val2,save_everystep=false,MaxIters=Iter)
sol3 = solve(prob,Rodas4() ,dt=dt_size,reltol=Val1,abstol=Val2,save_everystep=false,MaxIters=Iter)
sol4 = solve(prob,EulerHeun() ,dt=dt_size,reltol=Val1,abstol=Val2,save_everystep=false,MaxIters=Iter)
sol5 = solve(prob,OwrenZen5() ,dt=dt_size,reltol=Val1,abstol=Val2,save_everystep=false,MaxIters=Iter)
sol6 = solve(prob,TanYam7() ,dt=dt_size,reltol=Val1,abstol=Val2,save_everystep=false,MaxIters=Iter)
sol7 = solve(prob,alg_hints=[:nonstiff] ,dt=dt_size,reltol=Val1,abstol=Val2,save_everystep=false,MaxIters=Iter)
Array_sol1=sol1.u
Array_sol2=sol2.u
Array_sol3=sol3.u
Array_sol4=sol4.u
Array_sol5=sol5.u
Array_sol6=sol6.u
Array_sol7=sol7.u
Result_lim(Array_sol1)
Result_lim(Array_sol2)
Result_lim(Array_sol3)
Result_lim(Array_sol4)
Result_lim(Array_sol5)
Result_lim(Array_sol6)
Result_lim(Array_sol7)
end
Quanta_fin(-Delta_mw,Omega_mw)
I tried to use multiple solvers and even I decreased the step size and increased the number of iterations and as you can see the results are different (Rodas4() and alg_hints=[:stiff] have similar values), the main problem relies long time needed to achieve the steady state dynamics. That is why I want to ask how can I change the options of the solver to obtain a reliable solution and predict properly the steady state of the system. Thanks.
