1

我正在尝试对太阳系进行数值积分。我以前用普通的 Scheme 做过这个，现在我想用MIT 的非常有趣的 SCMUTILS-library来做这个。我做了什么：

  1. 我从喷气推进实验室获取了太阳系数据；即：太阳、水星、金星和地球在重心坐标中的质量、位置和速度。
  2. 我为微分方程编写了代码，这样系统中的每个物体（太阳、水星、金星、地球）都会以正确的方式被其他 3 个物体吸引。
  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
)
4

1 回答 1

3

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（麻省理工学院慷慨地开放获取）提供了两个例子，我认为你会发现：受限三体问题，它通过关注假设比另一个小得多的第三个物体来节省坐标二，并让坐标系旋转，使前两个物体位于 x 轴上。另一个例子是自旋轨道耦合，其中有两个物体，但卫星的不圆度是不可忽略的。这两种方法都给出了坐标比 3 *（物体数量）少得多的公式

于 2020-05-13T02:59:19.967 回答