It was developed specifically to solve geodesic equations in 5 dimensions, but is completely general in application, as long as the system is of the form¹

x_{1}^{(n)}
=
f(x_{1}, ẋ_{1}, …, x_{n}, …, t),
…

Advantages over other Mathematics packages include customizability of the integrator and of stopping conditions, and the natural introspective nature of lisp.

;;;; -*- Mode: Lisp; Syntax: ANSI-Common-Lisp; Base: 10; Package: integrate -*-
(in-package #:cl-user)
(

*defpackage*#:integrate (:use #:common-lisp) (:export #:*step-size* #:*tolerance* #:*integrator* #:*target-t* #:*log* #:*stop-function* #:precision-loss #:integrate #:integrated-1-function #:euler #:make-adaptive-integrator #:runge-kutta )) (in-package :integrate) ;;; trap for the unwary: *step-size* gets changed through the course ;;; of an integration. It might be wise to wrap calls to integrate ;;; inside (let ((*step-size* 0.001)) (integrate ...)) (*defvar**step-size* 1.0d-3) (*defvar**tolerance* 1.0d-3) ;;; BUG you moron. *epsilon* appears only to be used as a smoothing ;;; parameter in adaptive integrators. This is probably not clever. (*defvar**epsilon* single-float-epsilon) ;;; *integrator* api: ;;; takes a state, and integrates it from t (the first thingy in the ;;; state, usually 0) to whenever (funcall *stop-function* state) ;;; returns t, logging to the stream defined by *log*. Note that with ;;; higher-order adaptive integrators, it is perfectly possible that ;;; polynomial interpolation be necessary between the points produced. (*defvar**integrator* 'euler "An integrator takes a state (defined as a list, currently) and integrates it from time = (first state) to whenever (funcall *stop-function* state) does not return nil, with logging to the stream defined by *log*.") (*defvar**target-t* 1.0) (*defvar**stop-function* (*lambda*(state) (<= *target-t* (car state))) "The *stop-function* must return nil for the integration to continue.") (*defvar**log* nil "Stream (as interpreted by format) to log to. nil means no logging") ;;; we might eventually want integrate to return a function of the ;;; variable that can be evaluated as normal. (*define-condition*precision-loss (arithmetic-error) ()) (*defun*integrate (initial-time &rest system) "This is just a wrapper, to convert the system into something that will be understood by integrate-1. This is actually not ideal, as the logging then becomes not immediately transparent, but hey. We could always return a function that'll take the output log and return n functions representing the solutions :-)" ;; do this in two passes -- it's just easier. (handler-case (*let**((state (cons initial-time (*loop*for x in (cdr system) by #'cddr appending x)))) ; (nargs (length state))) (*let*((functions (cons (*lambda*(&rest x) (declare (ignore x)) 1) (*loop*for x on system by #'cddr for offset = 0 then (+ offset len) for len = (length (cadr x)) collect (car x) append (*loop*for y from 1 to (1- len) ;; with gensyms = (loop for i from 1 to nargs collect (gensym)) collect (*let*((n (+ y offset))) (*lambda*(&rest stuff) (nth n stuff)))))))) ;; (coerce `(lambda ,gensyms ;; (declare (ignorable ,@gensyms)) ;; ,(nth (+ y offset) gensyms)) ;; 'function)))))) (apply #'integrate-1 state functions))) (precision-loss () (warn "Before finishing, we've encountered loss of precision.") nil) (arithmetic-error (condition) (warn "~@(~a~) error signalled -- we probably hit a divergence." (type-of condition)) nil) )) ;;; document this. It is *so* cute you'll forget what on earth is ;;; going on. (*defun*integrate-1 (state &rest functions) ; this is probably too "cute" "State is a list of numbers, corresponding to (time y^{(n-1)} y^{(n-2)} ...), where the equation to solve is y^{(n)} = f(time, y^{(n-1)}, ..., y). This function performs a reduction of this nth order equation to n first order equations, adds an explicit equation for time (dt/dt = 1) as a check, and then applies the *integrator* function to the result" (*if*(not (= (length state) (length functions))) (*let*((gensyms (*loop*for i from 1 to (length state) collect (gensym)))) (apply #'integrate-1 state (apply #'list (*lambda*(&rest x) (declare (ignore x)) 1) (car functions) (*loop*for x from 3 to (length state) collect (coerce `(*lambda*,gensyms (declare (ignorable ,@gensyms)) ,(nth (- x 2) gensyms)) '*function*))))) (apply *integrator* state functions))) ;;; there's a macro that is just asking to be written somewhere in ;;; here... (*defun*linear-interpolate (x x0 y0 x1 y1) (*flet*((frob (x a b by) (* (/ (* (- x a)) (* (- b a))) by))) (+ (frob x x0 x1 y1) (frob x x1 x0 y0)))) (*defun*quadratic-interpolate (x x0 y0 x1 y1 x2 y2) (*flet*((frob (x a b c cy) (* (/ (* (- x a) (- x b)) (* (- c a) (- c b))) cy))) (+ (frob x x0 x1 x2 y2) (frob x x1 x2 x0 y0) (frob x x2 x0 x1 y1)))) (*defun*cubic-interpolate (x x0 y0 x1 y1 x2 y2 x3 y3) (*flet*((frob (x a b c d dy) (* (/ (* (- x a) (- x b) (- x c)) (* (- d a) (- d b) (- d c))) dy))) (+ (frob x x0 x1 x2 x3 y3) (frob x x1 x2 x3 x0 y0) (frob x x2 x3 x0 x1 y1) (frob x x3 x0 x1 x2 y2)))) ;;; A function that returns a function representing the solution to ;;; the single-variable differential equation. (*defun*integrated-1-function (state*function*) (*let*((alist nil)) (*lambda*(arg) (unless (assoc arg alist :test #'<) (*let*((*log* (make-string-output-stream)) ; should be a broadcast-stream? ;; this is wrong -- it should be something like (- arg last-value) (*step-size* (min *step-size* (/ arg pi))) (*target-t* arg)) (*if*alist ;; this depends on ordering (integrate-1 (car (last alist))*function*) (integrate-1 state*function*)) (*loop*with stream = (make-string-input-stream (get-output-stream-string *log*)) for line = (read-line stream nil 'no-more-lines) until (eq line 'no-more-lines) ;; preserve ordering -- this is inefficient. do (setf alist (sort (pushnew (read-from-string (format nil "(~a)" line)) alist :test #'equal) #'< :key #'car))))) (when (< (length alist) 2) (error "Need more points. Rethink.")) (*loop*for ooos = state then oos for oos = state then os for os = state then s for s in alist*if*(< (car os) arg (car s)) do (*cond*((< (car ooos) (car oos)) (return (apply #'cubic-interpolate (append (list arg) ooos oos os s)))) ((< (car oos) (car os)) (return (apply #'quadratic-interpolate (append (list arg) oos os s)))) (t (return (apply #'linear-interpolate (append (list arg) os s))))))))) (*defun*euler (state &rest functions) (*let*((current-state (copy-list state))) (*loop*until (funcall *stop-function* current-state)*if**log* ;; the benefit of this notation is that re-reading the ;; state is as simple as ;; (read-from-string (concatenate 'string "(" line ")")) ;; oh, and it's gnuplot friendly format, too. do (format *log* "~{~d~^ ~}~%" current-state) do (*let*((derivs (*loop*for fn in functions collect (apply fn current-state)))) (setf current-state (*loop*for y in current-state for d in derivs collect (+ y (* d *step-size*))))) finally (when *log* (format *log* "~{~d~^ ~}" current-state)) (return current-state)))) (*defun*make-adaptive-integrator (fn order) "takes two arguments -- an integrator to make adaptive and a number describing the order of the error (2 for euler, 5 for fourth-order runge-kutta)" (*lambda*(state &rest functions) (*let*((current-state (copy-list state))) (*loop*until (funcall *stop-function* current-state)*if**log* do (format *log* "~{~d~^ ~}~%" current-state) do (*let*((*log* nil)) (*loop*until (*let*((*target-t* (+ (first current-state) (* 1.9 *step-size*)))) (*let*((*stop-function* (*lambda*(state) (<= *target-t* (first state))))) (*let*((new-state-1 (*let*((*step-size* (* 2 *step-size*))) (apply fn current-state functions))) (new-state-2 (apply fn current-state functions))) ;; FIXME (*let*((deltas (mapcar (*lambda*(x y) (abs (- x y))) new-state-1 new-state-2))) (*if*(every (*lambda*(x) (< x *tolerance*)) deltas) (*progn*(setf *step-size* (* 0.99 *step-size* (expt (/ *tolerance* (+ *epsilon* (apply #'max deltas))) (/ 1 order)))) (setf current-state new-state-2)) (*progn*;; we need to retry with a smaller ;; step size, so return nil after ;; adjusting the step. (setf *step-size* (* 0.99 *step-size* (expt (/ *tolerance* (apply #'max deltas)) (/ 1 (1- order))))) nil)))))))) finally (when *log* (format *log* "~{~d~^ ~}" current-state)) (return current-state))))) (*defun*runge-kutta (state &rest functions) (*let*((current-state (copy-list state))) (*loop*until (funcall *stop-function* current-state)*if**log* do (format *log* "~{~d~^ ~}~%" current-state) do (setf current-state ;; those mapcars could be turned into map-intos, ;; not that this seems too slow in the first place. (*let**((derivs1 (*loop*for fn in functions collect (apply fn current-state))) (k1 (mapcar (*lambda*(x) (* *step-size* x)) derivs1)) (tmp-state (mapcar #'+ current-state (mapcar #'(*lambda*(x) (/ x 2)) k1))) (derivs2 (*loop*for fn in functions collect (apply fn tmp-state))) (k2 (mapcar (*lambda*(x) (* *step-size* x)) derivs2)) (tmp-state (mapcar #'+ current-state (mapcar #'(*lambda*(x) (/ x 2)) k2))) (derivs3 (*loop*for fn in functions collect (apply fn tmp-state))) (k3 (mapcar (*lambda*(x) (* *step-size* x)) derivs3)) (tmp-state (mapcar #'+ current-state k3)) (derivs4 (*loop*for fn in functions collect (apply fn tmp-state))) (k4 (mapcar (*lambda*(x) (* *step-size* x)) derivs4))) (*let*((steps (mapcar #'+ (mapcar (*lambda*(x) (/ x 6)) k1) (mapcar (*lambda*(x) (/ x 3)) k2) (mapcar (*lambda*(x) (/ x 3)) k3) (mapcar (*lambda*(x) (/ x 6)) k4)))) ; (format t "~a" (mapcar #'(lambda (x y) (if (= y 0) 1.0 (/ x y))) steps current-state)) (*if*(some (*lambda*(x) (< (abs x) double-float-epsilon)) (mapcar #'(*lambda*(x y) (*if*(= y 0) 1.0 (/ x y))) steps current-state)) (error 'precision-loss) (mapcar #'+ current-state steps)))) ) finally (when *log* (format *log* "~{~d~^ ~}" current-state)) (return current-state))))¹ *While this has to be some generalisation of geodesics equations, it's not clear which generalisation it is exactly. (ed.)*

Public Domain