为什么这个太阳系的数值积分保持运行？ （麻省理工学院计划 SCMUTILS）
Why does this numerical integration of the solar system keep running? (MIT-Scheme SCMUTILS)
我正在尝试对太阳系进行数值积分。我以前在简单的 Scheme 中做过这个，现在我想用非常有趣的 SCMUTILS-library from MIT 来做。我做了什么：
- 我从喷气推进实验室获取了太阳系数据；即：太阳、水星、金星和地球在重心坐标系中的质量、位置和速度。
- 我为微分方程编写了代码，这样系统中的每个物体（太阳、水星、金星、地球）都会以正确的方式被其他 3 个物体吸引。
- 我运行使用SCMUTILS通过数值积分进行模拟。
如果我 运行 模拟太阳 + 1 个其他行星，就可以了。如果我尝试拍摄太阳 + 其他 2 个行星，它似乎会挂起。这是 st运行ge，因为我 运行 几年前用我自己自制的 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
)
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 的魔力）。我认为这个特殊问题将是一个艰难的攀登。
[edit] SICM（麻省理工学院慷慨地开放获取）提供了两个例子，我想你会发现： restricted three-body problem which economizes on coordinates by focusing on the 3rd body which is assumed to be much smaller than the other two, and having the coordinate system rotate such that the first two bodies lie on the x axis. Another example is spin orbit coupling 其中有两个物体，但卫星的不圆度不是微不足道。这两种方法都给出了坐标比 3 *（物体数量）
少得多的公式
我正在尝试对太阳系进行数值积分。我以前在简单的 Scheme 中做过这个，现在我想用非常有趣的 SCMUTILS-library from MIT 来做。我做了什么：
- 我从喷气推进实验室获取了太阳系数据；即：太阳、水星、金星和地球在重心坐标系中的质量、位置和速度。
- 我为微分方程编写了代码，这样系统中的每个物体（太阳、水星、金星、地球）都会以正确的方式被其他 3 个物体吸引。
- 我运行使用SCMUTILS通过数值积分进行模拟。
如果我 运行 模拟太阳 + 1 个其他行星，就可以了。如果我尝试拍摄太阳 + 其他 2 个行星，它似乎会挂起。这是 st运行ge，因为我 运行 几年前用我自己自制的 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
)
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 的魔力）。我认为这个特殊问题将是一个艰难的攀登。
[edit] SICM（麻省理工学院慷慨地开放获取）提供了两个例子，我想你会发现： restricted three-body problem which economizes on coordinates by focusing on the 3rd body which is assumed to be much smaller than the other two, and having the coordinate system rotate such that the first two bodies lie on the x axis. Another example is spin orbit coupling 其中有两个物体，但卫星的不圆度不是微不足道。这两种方法都给出了坐标比 3 *（物体数量）
少得多的公式