As noted in the comments, when I ran your code (using LSODE in SciPy) I get
>>> result = odeint(f,y0,t)
lsoda-- warning..internal t (=r1) and h (=r2) are
such that in the machine, t + h = t on the next step
(h = step size). solver will continue anyway
in above, r1 = 0.2036188800000D+10 r2 = 0.1292689215500D-07
lsoda-- warning..internal t (=r1) and h (=r2) are
such that in the machine, t + h = t on the next step
(h = step size). solver will continue anyway
in above, r1 = 0.2036188800000D+10 r2 = 0.1292689215500D-07
lsoda-- at current t (=r1), mxstep (=i1) steps
taken on this call before reaching tout
in above message, i1 = 500
in above message, r1 = 0.2036188801501D+10
/usr/lib64/python2.7/site-packages/scipy/integrate/odepack.py:218: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
warnings.warn(warning_msg, ODEintWarning)
>>>
>>> for i in range(0,6):
... print(result[-1][i])
...
9.07638086329e+223
9.48797814966e-81
1.04657495534e-312
2.52059148034e-306
9.94646703738e+86
8.48798317088e-313
This is what I assume you must mean by the errors. I am not sure why it would do this, but there must be some bad default values for the tolerances, causing the timesteps to go really low with lsode. Multistep methods have to kick start using a low order method like Euler, and so their error estimators are really bad in the first few steps. I think it's going awry there.
I investigated this a bit using DifferentialEquations.jl. It has a bunch of algorithms ready, so it makes it a easier to do this kind of analysis. Setting up the problem is very similar to your Python code:
myu=398600.4418E+9
J2=1.08262668E-3
req=6378137
t0=86400*2.3567000000000000E+04
tN= 86400*2.3567250000000000E+04
y0 =[-9.0656779992979474E+05, -4.8397431127596954E+06, -5.0408120071376814E+06, -7.6805804020022015E+02, 5.4710987127502622E+03, -5.1022193482389539E+03]
function f(t,y,dy)
r2=(y[1]^2 + y[2]^2 + y[3]^2)
r3=r2^(3/2)
w=1+1.5J2*(req*req/r2)*(1-5y[3]*y[3]/r2)
dy[1] = y[4]
dy[2] = y[5]
dy[3] = y[6]
dy[4] = -myu*y[1]*w/r3
dy[5] = -myu*y[2]*w/r3
dy[6] = -myu*y[3]*w/r3
end
t = linspace(t0,tN)
using DifferentialEquations
prob = ODEProblem(f,y0,(t0,tN))
I don't think much of a comment needs to be said there, since it's almost a straight copy plus changing 0-based indexing to 1-based indexing. First I tried this with LSODA:
using LSODA
@time sol = solve(prob,lsoda(),abstol=1e-12)
#0.000690 seconds (17.74 k allocations: 373.578 KB)
println(sol[end])
#[1.20141e6,1.14517e6,7.13797e6,74.8469,-7231.03,1079.12]
When I run this with Sundials:
@time sol = solve(prob,CVODE_Adams(),abstol=1e-14)
#0.000730 seconds (8.87 k allocations: 187.641 KB)
println(sol[end])
#[1.1016e6,3.68894e6,6.3157e6,513.044,-6351.17,3642.44]
It seems like it works out alright. When I use the FORTRAN dopri
method (this is just a wrapper to the Hairer method):
using ODEInterfaceDiffEq
@time sol = solve(prob,dopri5(),abstol=1e-12)
#0.000885 seconds (7.83 k allocations: 180.469 KB)
println(sol[end])
[-965604.0,-1.84395e6,-5.65673e6,-292.419,7785.53,-2437.98]
You can see that errors build up more. I checked using a more modern Julia re-implementation of the Dormand-Prince 4/5 algorithm
@time sol = solve(prob,DP5(),abstol=1e-12)
#0.000467 seconds (6.82 k allocations: 155.828 KB)
println(sol[end])
[-9.60193e5,-1.96886e6,-5.61328e6,-319.081,7737.2,-2594.4]
It seems like this method is just very error prone. I would suggest using newer algorithms, like Tsit5()
and Vern7()
. You can see from the benchmarks that I posted that they are more efficient than the old Fortran codes, and I discuss this more here.
@time sol = solve(prob,Tsit5(),abstol=1e-12)
#0.000485 seconds (8.97 k allocations: 218.261 KB)
println(sol[end])
[1.03967e6,-2.41448e6,6.47286e6,-527.843,-7082.73,-2564.45]
@time sol = solve(prob,Vern7(),abstol=1e-13)
#0.001131 seconds (11.19 k allocations: 266.861 KB)
println(sol[end])
[1.10658e6,-1.5227e6,6.79846e6,-369.784,-7333.39,-1590.5]
So it looks like the higher order methods (CVODE_Adams
and Vern7
) are in the same direction. To see what the correct solution should be, we can re-solve it using higher-precision numbers. Here we will use BigFloat
s in the algorithms. We can do this just by changing the number types that we build the problem type with:
y0_big =big([-9.0656779992979474E+05, -4.8397431127596954E+06, -5.0408120071376814E+06, -7.6805804020022015E+02, 5.4710987127502622E+03, -5.1022193482389539E+03])
t0_big=big(86400*2.3567000000000000E+04)
tN_big=big(86400*2.3567250000000000E+04)
prob_big = ODEProblem(f,y0_big,(t0_big,tN_big))
And then solve the problem with higher precision numbers:
@time sol = solve(prob_big,Vern7(),abstol=1e-20,reltol=1e-20)
#5.124532 seconds (20.51 M allocations: 1020.807 MB, 29.95% gc time)
println(sol[end])
BigFloat[1.1053381662331820146336046671754933670595930944811486000087489128455927171218624720058718205588449974168956065779438050345e+06,-1.5416202537237067975875829305475944291202837487809289389890156042765665624065710316274224354945044382194625246683296005503e+06,6.7926473294913369302731863142467582412609188710382594274952998303466256241086526296652443665837323225497701006408603155527e+06,-3.7312828821058534895294734492513244053231049180643905427777841548594598586886829969951338572602652254542799098196952577885e+02,-7.3297961245134371967675775894582329860535459110035746501350371674509262650960485899126786425646757794442551813929572586948e+03,-1.61096155619704279518364795928351631960055375601824673571833246176656122912319418152730921855088538873761614569864638284e+03]
@time sol = solve(prob_big,Tsit5(),abstol=1e-20,reltol=1e-20)
#46.843028 seconds (202.64 M allocations: 9.848 GB, 21.61% gc time)
println(sol[end])
BigFloat[1.1053381662331820147067348726438400371709496335752162175134992231756774664804231225139143546847966065233311732755080517428e+06,-1.5416202537237067964033929090489735675629358124880169057038770565130982115037191409720907051171224745683301005104028510218e+06,6.7926473294913369306108500165460363424229267867186882151308990443549948472674738460107391298061753478752197008344571867852e+06,-3.7312828821058534874482439660747198554520591273413619386987734456375288506664244304076979140861220572371658224120877522505e+02,-7.3297961245134371970128642262106938950880875625957763378149307834830636800982964288295248480875424454937258559641758284114e+03,-1.6109615561970427939085757874332499247300266093609788206242850106668114350752364680610553537125616823780476995508356156203e+03]
This is confirmed by two algorithms at high precision, so this must be the value.
Recap
So let's recap a little bit. The timings are really to quick to say much, but at least in Julia there isn't a problem with getting high accuracy. The older methods like lsoda
and dopri5
have higher errors. Even DP5
. Methods and coefficients have been refined since the time these algorithms were made, which is why I wouldn't really recommend them anymore. Sundials and the newer Runge-Kutta methods in DifferentialEquations.jl (OrdinaryDiffEq.jl) are a lot more accurate.
There are a few sources of error that this shows. For one, each method has discretization error, which causes a difference between them. But that should go away with a small enough tolerance. This isn't necessarily the case because of floating point issues, the second source of error. It's a very long timespan, and so you need very low tolerances to not have a substantial buildup of error. This is confirmed by the BigFloat
results which show that we can only get the first few digits right if we go to much higher precision numbers (tolerances of about 1e-13
without getting junk, and we needed to go lower!).
But lastly, even after going to arbitrary precision numbers, putting a low tolerance, and essentially having two algorithms converge to the same result for the first 10 or so decimal places, we are still off from your proposed "true result". This is likely due to the fact that the initial condition was slightly inaccurate as well, which builds up over time to cause an error in the first digit. The error from the initial condition can grow very fast, so you need to make sure that is very accurate!
(However, I was not able to re-create the SciPy problems here... I am not sure why)
So you need to be careful about all 3 sources of errors, as demonstrated here.
Adding a plot for good measure.
plot(sol)

I don't know your problem, but you should always double-check the plot!
Can we get around the floating point errors?
The next question to address is, can we get around the floating point errors? Yes, we can rescale the problem. There is no reason to have such large numbers for every number! We can instead move every number down by something like 10^3 by changing units. I would suggest doing a change of units like that in both the dependent variables and time, and see if it can get better results with just Float64
's.
But most importantly, I hope this shows a good way to investigate the results.
Edit: More Verification
Because it was mentioned in the comments, I ran this with
@time sol = solve(prob_big,Vern9(),abstol=1e-25,reltol=1e-25)
println(sol[end])
BigFloat[1.1053381662331822596169776849641239592819519884844738459016726262029644062128757651516073415939693720405423642205909058064e+06,-1.5416202537237025194193402796238592848398830150314947716991093949659223769788623589030795132080240499488785780018224381598e+06,6.7926473294913380453538931349439465820220085181755818261687088166926884907951721678701419759674780893612240340931187099615e+06,-3.7312828821058461338614861968657625567908036695262270829525438975640521184553269025483381748908274236026477510068778050456e+02,-7.3297961245134380914206076512518039416894076161354630688132265052545303764963243830516528946565213622293265857749536691384e+03,-1.6109615561970383265735660519068721854582438082504069296606070807069042413835561362598331583440076713255077056006112545402e+03]
@time sol = solve(prob_big,Feagin14(),abstol=1e-25,reltol=1e-25)
println(sol[end])
BigFloat[1.1053381662331822596169777141870567234610511965557086835077483653533437799111462386523440877691562284278052381367094500505e+06,-1.5416202537237025194193397999916272100650863657026496400673463469591815591302649671026419236684086191564320676752769500241e+06,6.7926473294913380453538932693163486282736145849617871379561032756359556363271724242234264204395645390979862595520421289454e+06,-3.7312828821058461338614853616747368391649624257444346156270400632544621072526444156009615705008474094021665763109262437141e+02,-7.3297961245134380914206077511660934608408482056577325775582350370158926687300526721134971356964860345847978865296553096107e+03,-1.6109615561970383265735655400955954012952369987169718357788394884494500522058249689736956400471634474979232165731626140887e+03]
So many different algorithms are giving, at very high precision and accuracy, the same result to the first 14 decimal places. Thus for these inputs, this is a value you can trust. If this value is off, then it's the input numbers (the initial conditions, the constants, or the times) which are off.
Another Edit
We worked it out in the JuliaDiffEq Gitter chat channel. There was an error in the z-acceleration component. The new script is:
myu=398600.4418E+9
J2=1.08262668E-3
req=6378137
t0=86400*2.3567000000000000E+04
tN= 86400*2.3567250000000000E+04
y0 =[-9.0656779992979474E+05, -4.8397431127596954E+06, -5.0408120071376814E+06, -7.6805804020022015E+02, 5.4710987127502622E+03, -5.1022193482389539E+03]
function f(t,y,dy)
r2=(y[1]^2 + y[2]^2 + y[3]^2)
r3=r2^(3/2)
w=1+1.5J2*(req*req/r2)*(1-5y[3]*y[3]/r2)
w2=1+1.5J2*(req*req/r2)*(3-5y[3]*y[3]/r2)
dy[1] = y[4]
dy[2] = y[5]
dy[3] = y[6]
dy[4] = -myu*y[1]*w/r3
dy[5] = -myu*y[2]*w/r3
dy[6] = -myu*y[3]*w2/r3
end
t = linspace(t0,tN)
using DifferentialEquations
prob = ODEProblem(f,y0,(t0,tN))
sol = solve(prob,DP5(),abstol=1e-13,reltol=1e-13)
println(sol[end])
[1.09206e6,-1.97193e6,6.68529e6,-408.338,-7208.31,-2057.99]
which hits the true solution dead on. Credit to @crbinz