关于星历表的 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

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

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

这是诊断这些问题的好方法：尝试其他集成商。尝试 Rosenbrock23()、CVODE_BDF()、radau()、dopri5()、Vern9() 等。如果其中 none 有效，那么您将刚刚测试您的算法混合了经过最充分测试的算法（其中一些是 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)

精确到 1e-13，因为给定常量的大小很难精确到 1e5。如果你想在这个问题上达到这种精度，你需要调整你的单位或使用 BigFloat 数字。所以可能发生的事情是微分方程求解器意识到他们没有达到 1e-13 精度，缩小步长，并无限期地重复这个（因为他们永远不会达到 1e-13 精度，因为尺寸浮点数）直到步长太小并退出。如果您更改单位以使常量的大小更合理，则可以解决此问题。

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

Julia program on ephemerides shows inadequate answers

differential-equations

julia

calculation