I can solve the equations by ode45 in MATLAB, However, failed in julia with Warning: Instability detected. Aborting

i have use the same odefun in both matlab and julia,however, i got a warning in julia,the codes are as follows. I have changed the methods in julia and dimin,maxiters, but failed.

using DifferentialEquations
u0=[2.468249415089429,-0.22594658481728938,0.48341013085174944,-0.05587010826568186,-0.1878112092045923,0.20977961906960021,-8.568597818431394,0.0,-0.0,-0.10227692148167625,0.0,0.04764245128595939]
tspan=(0.0, 1433.00717532166)
function odeXe!(du,u,p,t::Float64)
#u[1:12]x1 vx1 y1 vy1 z1 vz1 …
soft=1.5

Ex=E0*(1.0*(t>=0.0 && t<10.0*T1)+(-t/(3.0*T1)+13.0/3.0)*(t>=10.0*T1 && t<=13.0*T1))*cos(omg1*t)
#Ex=E0*(1.0*(t>=0.0 && t<4.0*T1)+(-t/(2.0*T1)+3.0)*(t>=4.0*T1 && t<=6.0*T1))*cos(omg1*t)
Ey=0.0
Ez=0.0

Bx=-Ey/137.0
By=Ex/137.0
#By=0.0
Bz=0.0

r1=sqrt(u[1]^2+u[3]^2+u[5]^2+soft)
r2=sqrt(u[7]^2+u[9]^2+u[11]^2+soft)

ξ=1.3208
η=5.3128

diff_vr1=2.0/r1^2+52.0*(1-η/ξ+η/ξ*exp(ξ*r1)+η*r1*exp(ξ*r1))/((η/ξ*exp(ξ*r1)+1-η/ξ)^2*r1^2)
diff_vr2=2.0/r2^2+52.0*(1-η/ξ+η/ξ*exp(ξ*r2)+η*r2*exp(ξ*r2))/((η/ξ*exp(ξ*r2)+1-η/ξ)^2*r2^2)

#diff_vr1=diff_MHF(r1)
#diff_vr2=diff_MHF(r2)

ne1=diff_vr1/sqrt(u[1]^2+u[3]^2+u[5]^2+soft)
ne2=diff_vr2/sqrt(u[7]^2+u[9]^2+u[11]^2+soft)

e1e2=1.0/((u[1]-u[7])^2+(u[3]-u[9])^2+(u[5]-u[11])^2)/sqrt((u[1]-u[7])^2+(u[3]-u[9])^2+(u[5]-u[11])^2)

du[1]=u[2]
du[2]=-ne1*u[1]+(u[1]-u[7])*e1e2+(u[6]*By-u[4]*Bz)-Ex
du[3]=u[4]
du[4]=-ne1*u[3]+(u[3]-u[9])*e1e2+(u[2]*Bz-u[6]*Bx)
du[5]=u[6]
du[6]=-ne1*u[5]+(u[5]-u[11])*e1e2+(u[4]*Bx-u[2]*By)

du[7]=u[8]
du[8]=-ne2*u[7]+(u[7]-u[1])*e1e2+(u[12]*By-u[10]*Bz)-Ex
du[9]=u[10]
du[10]=-ne2*u[9]+(u[9]-u[3])*e1e2+(u[8]*Bz-u[12]*Bx)
du[11]=u[12]
du[12]=-ne2*u[11]+(u[11]-u[5])*e1e2+(u[10]*Bx-u[8]*By)

end
tra=ODEProblem(odeXe!,u0,tspan)
sol=solve(tra,Vern8(),reltol=1e-9, abstol=1e-9,save_everystep=false,save_start=false)

Have a look here


It shows you what Julia solver corresponds to your old matlab solver. You can also call the same matlab solver from julia to verify that it computes the same solution, if it doesn’t, you might have made an error in your translation.
5 Likes