Julia 星历表程序显示出不充分的答案

Question

在求解卫星运动微分方程时遇到此错误:

dt <= dtmin. Aborting. If you would like to force continuation with dt=dtmin, set force_dtmin=true

这是我的代码:

    using JPLEphemeris
    spk = SPK("some-path/de430.bsp")
    jd = Dates.datetime2julian(DateTime(some-date))#date of the calculatinons 
    yyyy/mm/dd hh/mm/ss
    jd2 = Dates.datetime2julian(DateTime(some-date))#date of the calculatinons 
    yyyy/mm/dd hh/mm/ss

   println(jd)
println(jd2)

st_bar_sun = state(spk, 0, 10, jd)
st_bar_moon_earth = state(spk, 0, 3, jd)
st_bar_me_earth = state(spk, 3, 399, jd)
st_bar_me_moon = state(spk, 3, 301, jd)

moon_cord = st_bar_me_moon - st_bar_me_earth

a = st_bar_moon_earth + st_bar_me_earth
sun_cord = st_bar_sun - a


println(sputnik_cord)
sputnik_cord = [8.0,8.0,8.0,8.0,8.0,8.0,8.0]
moon_sputnik = sputnik_cord - moon_cord
sun_sputnic = sputnik_cord - sun_cord

Req = 6378137
J2 = 1.08262668E-3

GMe = 398600.4418E+9
GMm = 4.903E+12
GMs = 1.32712440018E+20

function f(dy,y,p,t)

  re2=(y[1]^2 + y[2]^2 + y[3]^2)
  re3=re2^(3/2)

  rs3 = ((y[1]-sun_cord[1])^2 + (y[2]-sun_cord[2])^2 + (y[3]-sun_cord[3])^2)^(3/2)
  rm3 = ((y[1]-moon_cord[1])^2 + (y[2]-moon_cord[2])^2 + (y[3]-moon_cord[3])^2)^(3/2)

  w = 1 + 1.5J2*(Req*Req/re2)*(1 - 5y[3]*y[3]/re2)
  w2 = 1 + 1.5J2*(Req*Req/re2)*(3 - 5y[3]*y[3]/re2)

  dy[1] = y[4]
  dy[2] = y[5]
  dy[3] = y[6]
  dy[4] = -GMe*y[1]*w/re3 
  dy[5] = -GMe*y[2]*w/re3 
  dy[6] = -GMe*y[3]*w2/re3
end

function f2(dy,y,p,t)

  re2=(y[1]^2 + y[2]^2 + y[3]^2)
  re3=re2^(3/2)

  rs3 = ((y[1]-sun_cord[1])^2 + (y[2]-sun_cord[2])^2 + (y[3]-sun_cord[3])^2)^(3/2)
  rm3 = ((y[1]-moon_cord[1])^2 + (y[2]-moon_cord[2])^2 + (y[3]-moon_cord[3])^2)^(3/2)

  w = 1 + 1.5J2*(Req*Req/re2)*(1 - 5y[3]*y[3]/re2)
  w2 = 1 + 1.5J2*(Req*Req/re2)*(3 - 5y[3]*y[3]/re2)

  dy[1] = y[4]
  dy[2] = y[5]
  dy[3] = y[6]
  dy[4] = -GMe*y[1]*w/re3 - GMm*y[1]/rm3 - GMs*y[1]/rs3
  dy[5] = -GMe*y[2]*w/re3 - GMm*y[2]/rm3 - GMs*y[2]/rs3
  dy[6] = -GMe*y[3]*w2/re3- GMm*y[3]/rm3 - GMs*y[3]/rs3
end



y0 = sputnik_cord

jd=jd*86400
jd2=jd2*86400
using DifferentialEquations
prob = ODEProblem(f,y0,(jd,jd2))
sol = solve(prob,DP5(),abstol=1e-13,reltol=1e-13)


prob2 = ODEProblem(f2,y0,(jd,jd2))
sol2 = solve(prob2,DP5(),abstol=1e-13,reltol=1e-13)
println("Without SUN and MOON")
println(sol[end])
for i = (1:6)
 println(sputnik_cord[i]-sol[end][i])
end
println("With SUN and MOON")
println(sol2[end])

什么(除了值)可能是造成这种情况的原因？在我在函数 f2 的 dy[4]、dy[5]、dy[6] 定义中添加 sun_coords 和 Moon_coords 项之前，它运行良好(我认为函数 f1 工作正常)。

Answer 1

发生这种情况有两个原因。其一，您可能会看到此错误，因为模型由于实现问题而不稳定。如果您不小心犯了错误，解决方案可能会发散到无穷大，并且随着发散，时间步长会缩短，并且存在此错误。

可能发生的另一件事是您的模型可能很僵硬。如果不同组件之间存在较大的时间尺度差异，则可能会发生这种情况。在这种情况下，显式 Runge-Kutta 方法 DP5() 不是解决该问题的合适算法。相反，您会想看一些东西 for stiff equations 。我建议尝试一下 Rosenbrock23():它不是最快的，但它非常稳定，如果问题是可积的，它会处理它。

这是诊断这些问题的一个非常好的方法:尝试其他集成商。尝试Rosenbrock23()、CVODE_BDF()、radau()、dopri5()、Vern9( ) 等。如果这些都不起作用，那么您将刚刚使用经过最充分测试的算法(其中一些是 Julia 实现，但其他只是标准经典 C++ 和Fortran 方法)，这表明问题出在您的模型公式中，而不是该问题的特定求解器的特性。由于我无法运行您的示例(您应该使您的示例可运行，即不需要额外的文件，如果您希望我测试的话)，我无法确定您的模型实现是否正确，这将是找出答案的好方法.

我的猜测是，由于浮点问题，您写下的模型并不是一个好的实现。

GMe = 398600.4418E+9
GMm = 4.903E+12
GMs = 1.32712440018E+20

这些值的精度仅比规定值小 16 位:

> eps(1.32712440018E+20)
16384.0
> 1.32712440018E+20 + 16383
1.3271244001800002e20

> 1.32712440018E+20 + 16380
1.3271244001800002e20

> 1.32712440018E+20 + 16000
1.3271244001800002e20

请注意，该值的精度低于机器 epsilon。好吧，你要求的是

sol = solve(prob,DP5(),abstol=1e-13,reltol=1e-13)

当给定常量的大小很难精确到 1e5 时，精度为 1e-13。如果您想在这个问题上达到这种精度，您需要调整单位或使用 BigFloat 数字。因此，可能发生的情况是，微分方程求解器意识到它们没有达到 1e-13 精度，因此缩小了步长，并无限地重复此操作(因为它们永远无法达到 1e-13 由于 float 的大小而导致精度)，直到步长太小并退出。如果您更改单位，使常量的大小更加合理，则可以解决此问题。