Why does this numerical integration of the solar system keep running? (MIT-Scheme SCMUTILS)

Question

I'm trying to do a numerical integration of the solar system. I did this before in plain Scheme, now I want to do it using the very interesting SCMUTILS-library from MIT. What I did:

1. I took solar system data from the Jet Propulsion Laboratory; that is: the mass, the position and the velocity of the sun, mercurius, venus and earth in barycentric coordinates.
2. I wrote code for the differential equation, such that every object in the system (sun, mercurius, venus, earth) gets attracted by the 3 others in the correct way.
3. I ran the simulation through numerical integration using SCMUTILS.

If I run the simulation with the sun + 1 other planet, it works. If I try to take the sun + 2 other planets, it seems to hang. This is strange as I ran the simulation with the same data a few years ago with my own home-built Runge-Kutta integrator, and that worked fine.

Note that I'm well-known with MIT-Scheme and numerical integration, and that I only want to learn SCMUTILS. I'm doing something wrong clearly, and it would surprise me if this problem couldn't be tackled with SCMUTILS.

Also, I'm not fixed on my code: if somebody can provide me with a working implementation in SCMUTILS, then that's fine as well, as long as I understand what I'm doing wrong in my program. I just want to use SCMUTILS in an idiomatic way...

My code is below (about 60 well-documented lines). Thanks for any comments or improvements towards a working simulation.

;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

The way that scmutils handles the integration is interesting. The state derivative function works with a local tuple as described in SICM, but the integrator wants to work with a function that takes an array of floats as input and produces an array of equal size as output. To do this, scmutils takes the initial state data and replaces the values in it with symbols and passes this to your derivative. This produces symbolic output, which can be used to prepare a function with the right signature for the integrator. (I can describe this process in greater detail if you would like).

Your problem is in Cartesian coordinates, however, and the resulting symbolic expression is hairy. You can see this process in action by creating your own symbolic state and passing it to the derivative function, and simplifying the output (by passing the result through pe (print expression)):

(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))

The result is immense, so I haven't pasted it here. If you now add another planet, by uncommenting the lines with subscript 2, you will find the print hangs, which means the expression simplifier has bogged down. The numerical integrator hasn't even run yet.

What to do? You might be able to get some capacity back by eliminating the z coordinate. You could move the GM parameters, which are constant, to the argument list of the state derivative constructor function, leaving only things that will change in the state tuple itself. You could flatten the state tuple a little; its structure is entirely up to you.

Ultimately, though, the function that gets integrated will be much more complicated than the function you would have written yourself, and a lot of that has to do with the sqrt(x^2 + y^2 + ...) terms that you get from the Cartesian coordinates. Scmutils was designed for problems where the use of generalized coordinates produces compact Lagrangians, from which simpler state derivative functions can be derived (automatically, which is the magic of scmutils). I think this particular problem would be an uphill climb.

[edit] SICM (which MIT has graciously made open-access) provides two examples I think you will find illustrative: the 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 in which there are two bodies, but the satellite's out of roundness is not negligible. Both of these approaches give formulations where there are many fewer coordinates than 3 * (number of bodies)