为什么这个太阳系的数值积分一直在运行？(MIT-scme SCMUTILS)

Question

我想做一个太阳系的数值积分。我以前用普通的Scheme做过，现在我想用非常有趣的 "太阳系数值积分 "来做。麻省理工学院的SCMUTILS-库。. 我做了什么。

1. 我从喷气推进实验室获取了太阳系数据，即：太阳、墨丘利斯、金星和地球的质量、位置和速度的双中心坐标。
2. 我写了微分方程的代码，使系统中的每个物体（太阳、水银星、金星、地球）都能以正确的方式被其他三个物体吸引。
3. 我使用SCMUTILS通过数值积分运行模拟。

如果我用太阳+1个其他星球来运行模拟，它就会工作。如果我试着把太阳+2颗其他行星，它似乎挂了。这很奇怪，因为几年前我用同样的数据，用自己自制的Runge-Kutta积分器运行模拟，结果运行良好。

请注意，我对MIT-Scheme和数值积分很熟悉，我只想学习SCMUTILS。我显然是做错了什么，如果这个问题不能用SCMUTILS解决，我会很惊讶。

另外，我的代码并不固定：如果有人能给我提供SCMUTILS中的工作实现，那也可以，只要我明白我的程序中做错了什么。 我只是想用一种习惯的方式来使用SCMUTILS...。

我的代码在下面（大约60行文件）。感谢任何意见或改进，以实现工作模拟。

;JPL-DE - Ephemerides from Jet Propulsion Laboratory http://ssd.jpl.nasa.gov
(define solar-system
  (up (+ 0.0 (decoded-time->universal-time (make-decoded-time 0 0 0 1 1 2000 0))) ; January 1st 2000 at 0h0m0s UT
      (up (up 1.3271244004193937e20                                               ; Sun mass * gravitational constant                    
          (up -1068000648.30182 -417680212.56849295 30844670.2068709)             ; Sun position (x,y,z) in meters in barycentric coordinates
          (up 9.305300847631916 -12.83176670344807 -.1631528028381386))           ; Sun velocity (vx,vy,vz) in meters per second
      (up 22031780000000.                                                         ; Mercurius
          (up -22120621758.62221 -66824318296.10253 -3461601353.17608) 
          (up 36662.29236478603 -12302.66986781422 -4368.33605178479)) 
      (up 324858592000000.                                                        ; Venus
          (up -108573550917.8141 -3784200933.160055 6190064472.97799)
          (up 898.4651054838754 -35172.03950794635 -532.0225582712421))
;     (up 398600435436000.                                                        ; Earth
;         (up -26278929286.8248 144510239358.6391 30228181.35935813)
;         (up -29830.52803283506 -5220.465685407924 -.1014621798034465))
      )))

(define (ss-time state)       ; Select time from solar system state
  (ref state 0))
(define (ss-planets state)    ; Select up-tuple with all planets
  (ref state 1))
(define (ss-planet state i)   ; Select i-th planet in system (0: sun, 1: mercurius, 2: venus, 3: earth) (Note: the sun is also a "planet")
  (ref (ss-planets state) i))
(define (GM planet)           ; Select GM from planet (GM is gravitational constant times mass of planet)
  (ref planet 0))
(define (position planet)     ; Select position up-tuple (x,y,z) from planet
  (ref planet 1))
(define (velocity planet)     ; Select velocity up-tuple (vx,vy,vz) from planet
  (ref planet 2))

(define ((dy/dt) state)
  (define (gravitational-force on-planet by-planet)              ; Calculate gravitational force on planet "on-planet" by "by-planet"
    (if (equal? on-planet by-planet)                             ; Compare planets
        (up 0.0 0.0 0.0)                                         ; There is no force of a planet on itself
        (let* ((r (- (position on-planet) (position by-planet))) ; Position of "on-planet" seen from "by-planet"
               (d (abs r)))                                      ; Distance between the two
          (* -1.0 (GM by-planet) (/ r (* d d d))))))             ; Gravitational force is negatively directed, we divide by d^3 to cancel out r in nominator
  (define (dy/dt-planet planet)                                                                                         ; Calculate dy/dt for a planet
    (up 0.0                                                                                                             ; No change in GM
        (velocity planet)                                                                                               ; Change in position is velocity
        (* (s:generate (s:length (ss-planets state)) 'up (lambda (i) (gravitational-force planet (ss-planet state i)))) ; Calculate gravitation forces from
           (s:generate (s:length (ss-planets state)) 'down (lambda (i) 1.0)))))           ; all other planets, and sum them via inner product with down-tuple
  (up 1.0                                                                                              ; Timestep: dt/dt=1.0
      (s:generate (s:length (ss-planets state)) 'up (lambda (i) (dy/dt-planet (ss-planet state i)))))) ; Calculate dy/dt for all planets

(define win (frame -150e9 150e9 -150e9 150e9 512 512)) ; Viewport: a square of 300e9 meters by 300e9 meters plotted in a 512 by 512 window

(define ((monitor-xy) state)
  ((s:elementwise (lambda (planet) (plot-point win (ref (position planet) 0) (ref (position planet) 1)))) (ss-planets state))) ; Plot X,Y (no Z) for planets 

(define end                                                             ; Define end state
  ((evolve dy/dt)                                                       ; Run simulation
   solar-system                                                         ; Starting state, Jan. 1st 2000 0h0m0s
   (monitor-xy)                                                         ; Plot positions throughout simulation
   (* 1.0 60 60)                                                        ; Timestep: 1 hour
   (decoded-time->universal-time (make-decoded-time 0 0 0 1 1 2005 0))) ; Run for 5 years
)

Answer 1

scmutils处理积分的方式很有趣。正如SICM中描述的那样，状态导数函数使用一个局部元组，但是积分器希望使用一个函数，将一个浮动数组作为输入，并产生一个大小相等的数组作为输出。要做到这一点，scmutils会接收初始状态数据，并用符号替换其中的值，并将其传递给你的导数。这就产生了符号输出，它可以用来为积分器准备一个具有正确签名的函数。如果你愿意，我可以更详细地描述这个过程）。

然而，你的问题是在笛卡尔坐标中，由此产生的符号表达式是毛毛雨。您可以通过创建您自己的符号状态并将其传递给导数函数，然后简化输出（通过将结果通过 pe 打印表达式））。

(define symbolic-system
  (up 't
      (up (up 'g_0
          (up 'x_0 'y_0 'z_0)             ; Sun position (x,y,z) in meters in barycentric coordinates
          (up 'vx_0 'vy_0 'vz_0))           ; Sun velocity (vx,vy,vz) in meters per second
      (up 'g_1
          (up 'x_1 'y_1 'z_1)
          (up 'vx_1 'vy_1 'vz_1))
;     (up 'g_2
;          (up 'x_2 'y_2 'z_2)
;          (up 'vx_2 'vy_2 'vz_2))
;     (up 'g_3
;          (up 'x_3 'y_3 'z_3)
;          (up 'vx_3 'vy_3 'vz_3))
      )))

(pe ((dy/dt) symbolic-system))

结果是巨大的，所以我没有把它粘贴在这里。如果你现在再加一个星球，把带下标2的那几行去掉，你会发现打印挂了，这说明表达式简化器已经卡壳了。数值积分器还没运行呢。

怎么办呢？你也许可以通过消除z坐标来找回一些容量。你可以把GM参数，也就是常数，移到状态导数构造函数的参数列表中，只留下状态元组本身会变化的东西。你可以将状态元组扁平化一点，它的结构完全由你决定。

不过最终，被整合的函数会比你自己写的函数要复杂得多，这其中有很多是与你的 sqrt(x^2 + y^2 + ...) 术语，你可以从笛卡尔坐标中得到。Scmutils是为使用广义坐标产生紧凑拉格朗日的问题而设计的，从这些问题中可以推导出更简单的状态导数函数（自动的，这就是scmutils的神奇之处）。我认为这个特殊的问题将是一个上坡。

[编辑]SICM（麻省理工学院很慷慨地开放了它的访问权限）提供了两个例子，我想你会发现这两个例子很有说明意义：SICM中的 局部三体问题 它通过将焦点放在第3个体上，假设第3个体比其他两个体小得多，并使坐标系旋转，使前两个体位于x轴上，从而节省了坐标。另一个例子是 自旋轨道耦合 其中有两个体，但卫星的出圆度是不可忽视的。这两种方法都给出了坐标数比3*（天体数）少很多的公式。