A package to aid in solving differential equations, currently slightly limited in scope.

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¹

x1(n) = f(x1, ẋ1, …, xn, …, 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