Numerical Integration

[jstoddard]

Because doing database-driven programming for businesses gets boring (but hey, it's a living).

This solves the equation of motion r'' = GM/r^2 for any number of point masses in two dimensions. These actually do each step of the integration:

Code: Select all

(defparameter *G* 6.67438d-11)

(defun leapfrog-accel (x y masses x-positions y-positions)
  "Calculates acceleration of a particle at x,y due to gravity from
particles of masses and positions. a = -Gm2q2/r^2"
  (list
   (loop for other-mass in masses
      for particle-index upfrom 0
      when (not (= (aref x-positions particle-index) x))
      ;; Gm2cos(theta)/((x2-x1)^2 + (y2-y1)^2)
      sum (/ (* *G* other-mass (- (aref x-positions particle-index) x))
	     (expt (+ (expt (- (aref x-positions particle-index) x) 2)
		      (expt (- (aref y-positions particle-index) y) 2))
		   1.5d0)))
   (loop for other-mass in masses
      for particle-index upfrom 0
      when (not (= (aref y-positions particle-index) y))
      ;; Gm2sin(theta)/((x2-x1)^2 + (y2-y1)^2)
      sum (/ (* *G* other-mass (- (aref y-positions particle-index) y))
             (expt (+ (expt (- (aref x-positions particle-index) x) 2)
                      (expt (- (aref y-positions particle-index) y) 2))
                   1.5d0)))))

(defun leapfrog-velocity (v-x v-y a-x a-y time-step)
  "Calculates velocity update of one-half step."
  (list
   (+ v-x (/ (* a-x time-step) 2))
   (+ v-y (/ (* a-y time-step) 2))))

(defun leapfrog-position (x y v-x v-y time-step)
  "Calculates position update of one time step."
  (list
   (+ x (* v-x time-step))
   (+ y (* v-y time-step))))

Then some utility functions for what follows:

Code: Select all

(defun first-elements (list)
  "Return a list of the car of each sublist of list."
  (let ((outlist nil))
    (dolist (element list outlist) (setf outlist (nconc outlist (list (car element)))))))

(defun second-elements (list)
  "Return a list of the second element of each sublist of list."
  (let ((outlist (nconc)))
  (dolist (element list outlist) (setf outlist (nconc outlist (list(second element)))))))

(defmacro new-axis-array (num-elements array-type array-axis)
  "Create an array for a set of coordinates, velocities or accelerations."
  `(make-array ,num-elements
	       :initial-contents (,@(cond ((= array-axis 2) `(second-elements ,array-type))
					   (t `(first-elements ,array-type))))))

And the loop to integrate the entire system in all it's poorly written glory:

Code: Select all

(defun leapfrog-integrate (mass position velocity
			   &key (time-step 1/32) (number-of-steps 256) (output-modulus 4))
  "Leapfrog integrate a system of particles. mass is a list of masses for each particle,
position velocity and acceleration is a list of lists in the form (x y) for each particle,
time-step is the time interval for integration, number of steps is the number of time
steps integrated, and output-modulus should be set to n in order to output results every
n steps."
  (let* ((num-particles (length mass))
	 (acceleration (loop for i from 1 to num-particles collecting '(0.0d0 0.0d0)))
	 (x-positions (new-axis-array num-particles position 1))
	 (y-positions (new-axis-array num-particles position 2))
	 (x-velocities (new-axis-array num-particles velocity 1))
	 (y-velocities (new-axis-array num-particles velocity 2))
	 (x-accels (new-axis-array num-particles acceleration 1))
	 (y-accels (new-axis-array num-particles acceleration 2)))
    (dotimes (step number-of-steps)
      (dotimes (particle num-particles)
	(let ((vect (leapfrog-accel
		     (aref x-positions particle) (aref y-positions particle) mass
		     x-positions y-positions)))
	  (setf (aref x-accels particle) (first vect))
	  (setf (aref y-accels particle) (second vect)))
	(let ((vect (leapfrog-velocity (aref x-velocities particle) (aref y-velocities particle)
				       (aref x-accels particle) (aref y-accels particle) time-step)))
	  (setf (aref x-velocities particle) (first vect))
	  (setf (aref y-velocities particle) (second vect)))
	(let ((vect (leapfrog-position (aref x-positions particle) (aref y-positions particle)
				       (aref x-velocities particle) (aref y-velocities particle)
				       time-step)))
	  (setf (aref x-positions particle) (first vect))
	  (setf (aref y-positions particle) (second vect)))
	(cond ((= 0 (mod step output-modulus))
	       (format t "~a ~a ~a ~a ~a ~a~%"
		       (aref x-positions particle) (aref y-positions particle)
		       (aref x-velocities particle) (aref y-velocities particle)
		       (aref x-accels particle) (aref y-accels particle))))))
      (dotimes (particle num-particles) ; update accels, bump up velocities second half-step
	(let ((vect (leapfrog-accel
		     (aref x-positions particle) (aref y-positions particle) mass
		     x-positions y-positions)))
	  (setf (aref x-accels particle) (first vect))
	  (setf (aref y-accels particle) (second vect)))
	(let ((vect (leapfrog-velocity (aref x-velocities particle) (aref y-velocities particle)
				       (aref x-accels particle) (aref y-accels particle) time-step)))
	  (setf (aref x-velocities particle) (first vect))
	  (setf (aref y-velocities particle) (second vect))))))

This is a bit complicated, so to test out the algorithm, I did something simpler (which doesn't quite work w/ the default parameters -- you'll need adjust *G* a few orders of magnitude, I think)

Code: Select all

(defun planet-sun (&key (solar-mass 1.9891d30) (x 149598261) (y 0) (x-velocity 0) (y-velocity 29780)
		   (time-step 1/2) (number-of-steps 63115200) (output-modulus 1209600))
  "Calculate a planet's trajectory around the sun via leapfrog integration. The values
default to earth's, calculated every half second for one year, and output once per week."
  (let ((x-positions (make-array 2 :initial-contents (list 0.0d0 x)))
	(y-positions (make-array 2 :initial-contents (list 0.0d0 y)))
	(masses (list solar-mass 0.0d0))
	(x-accel 0) (y-accel 0))
    (dotimes (step number-of-steps)
      (let ((vect (leapfrog-accel (aref x-positions 1) (aref y-positions 1) masses
				  x-positions y-positions)))
	(setf x-accel (first vect))
	(setf y-accel (second vect)))
      (let ((vect (leapfrog-velocity x-velocity y-velocity x-accel y-accel time-step)))
	(setf x-velocity (first vect))
	(setf y-velocity (second vect)))
      (let ((vect (leapfrog-position (aref x-positions 1) (aref y-positions 1)
				     x-velocity y-velocity time-step)))
	(setf (aref x-positions 1) (first vect))
	(setf (aref y-positions 1) (second vect)))
      (let ((vect (leapfrog-accel (aref x-positions 1) (aref y-positions 1) masses
				  x-positions y-positions)))
	(setf x-accel (first vect))
	(setf y-accel (second vect)))
      (let ((vect (leapfrog-velocity x-velocity y-velocity x-accel y-accel time-step)))
	(setf x-velocity (first vect))
	(setf y-velocity (second vect)))
      (cond ((= 0 (mod step output-modulus))
	     (format t "~a ~a ~a ~a ~a ~a~%" (aref x-positions 1) (aref y-positions 1)
		     x-velocity y-velocity x-accel y-accel))))))

But I think it's generally correct, as I get an ellipse out of simple parameters:

Code: Select all

INTEGRATION> (setf *G* 1)
1
INTEGRATION> (planet-sun :solar-mass 10 :x 10 :y 0 :x-velocity 0 :y-velocity 1 :time-step 1/16 :number-of-steps 640 :output-modulus 64)
9.9998046875d0 1/16 -0.006249938963055646d0 0.9999804687500112d0 -0.09999804681778067d0 -6.249999996423722d-4
9.186089263033907d0 3.951698750984734d0 -0.3951657018967911d0 0.9186089910251384d0 -0.09186067358807011d0 -0.03951689328160332d0
6.922091396907244d0 7.2170112821823444d0 -0.7216919111046535d0 0.6922100649290486d0 -0.06922028978702764d0 -0.0721694620461256d0
3.565253368118854d0 9.342923804059755d0 -0.9342773855923646d0 0.35652938047010857d0 -0.03565186157694772d0 -0.09342747675800077d0
-0.35445209922215987d0 9.99381739199921d0 -0.9993618744828243d0 -0.03543473975799068d0 0.0035444134707291216d0 -0.09993514233924003d0
-4.218200083339988d0 9.066949048541678d0 -0.9066734615521503d0 -0.42180026676992666d0 0.04218024380824546d0 -0.09066571378985058d0
-7.416013147836746d0 6.7086623517255655d0 -0.6708482742985552d0 -0.7415716947720666d0 0.07415634773903634d0 -0.06708320067683003d0
-9.443056196378434d0 3.2912743920025047d0 -0.32911797599226766d0 -0.9442687037157195d0 0.09442518328703713d0 -0.0329108692408245d0
-9.979326898219654d0 -0.6457064570371644d0 0.06456875569459729d0 -0.9978937095749154d0 0.09978742811230788d0 0.006456686640332846d0
-8.940166639624119d0 -4.480749896857714d0 0.4480623272953247d0 -0.89398163315686d0 0.08939670592949832d0 0.04480501281683719d0
NIL

First, I'd like to Lisp-ify the code some. Any suggestions? Should I be using arrays to hold the current state of the system, or just stick to lists? If I stick to arrays I realize I should probably use 2-dimensional arrays instead of separate for x & y, but somehow the obvious slipped past me while I was coding...

Second, are there any n-body toolkits for Common Lisp to do this kind of stuff? A quick Google search didn't really turn up anything. Computational physics seems to be dominated by C and Fortran, and Java for teaching. C and Fortran are probably a bit more efficient, and this stuff is computationally intensive, but it still seems like fun to do it in Lisp. If not much is available, I want to look into writing some Common Lisp tools similar to a few of the tools found in NEMO http://carma.astro.umd.edu/nemo/, and maybe hack together some very simple visualization routines with LispbuilderSDL. Just enough to simulate some galaxy collisions because, hey, it's fun. We'll see where it goes from there (and what I can actually get done in my free time).

[FAU]

Well I don't know anything about that but maybe the MAXIMA system is interesting to you as apparently you want to do lots of math stuff.

http://maxima.sourceforge.net/

[Warren Wilkinson]

You say integration -- but your integration function seems to be simply a stepping from one state to the next.

Do you mean integration as in the inverse of a derivative?

[jstoddard]

You say integration -- but your integration function seems to be simply a stepping from one state to the next.

Do you mean integration as in the inverse of a derivative?

Yes and yes. Numerical integration, meaning using a discrete timestep (delta-t) to estimate the results where integrating analytically would be difficult or impossible. We're solving the differential equation r'' = GM/r^2 by multiplying r'' (acceleration) by our timestep to estimate r' (velocity), then again for speed. Keeping the timestep small and iterating over and over again gets us an approximation of the r' = integral (r''*dt).

It's an estimate, sure, but with a proper integration algorithm and a well-chosen timestep we can get very good results; good enough to run physics simulations.

[bit-player]

You might want to have a look at this book:

Structure and Interpretation of Classical Mechanics. Gerald Jay Sussman and Jack Wisdom with Meinhard E. Mayer. xxi + 534 pp. MIT Press. 2001.

The programs are written in Scheme rather than Common Lisp, but that's a lot closer than Fortran or Java.

[nuntius]

Sundials is an industrial-strength integration tool.

There are a couple simple integration tools in CL; but I recommend a CFFI binding to sundials.

[jstoddard]

I also found this today: The GNU Scientific Library for Lisp http://common-lisp.net/project/gsll/.