Choosing suitable options for solver ODE

Specific Domains Numerics
, ,
1

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.

2

If you have a lower error tolerance for the solution, then decrease the error tolerance. And this would run hundreds of times faster if you didn’t use non-constant globals.

1 Like