Numerical Integration
Posted: Sun Feb 20, 2011 9:25 am
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:
Then some utility functions for what follows:
And the loop to integrate the entire system in all it's poorly written glory:
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)
But I think it's generally correct, as I get an ellipse out of simple parameters:
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).
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))))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))))))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))))))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))))))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
NILSecond, 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).