share/contrib/symplectic_ode/symplectic_ode.lisp

   1 ;;;; Common Lisp/Maxima code for symplectic integration of time independent separable hamiltonian
   2 ;;;; systems.
   3
   4 ;;;; Copyright (C) 2017, Barton Willis <willisb@unk.edu>
   5 ;;;; This work is licensed under a Creative Commons Attribution 4.0 International License.
   6
   7 (in-package :maxima)
   8
   9 ;;; Let's have version numbers 1,2,3, ...
  10 ($put "symplectic_ode" 1 '$version)
  11
  12 ;;; Define a mlist msimpind. This definition should be somewhere else.
  13 (setf (get 'mlist 'msimpind) (list 'mlist 'simp))
  14
  15 ;;; A Maxima coerce function. This function is adequate for the needs of symplectic_ode, but
  16 ;;; it might be useful to extend it to additional types.
  17 (defun mcoerce (x ntype)
  18   "Either convert x to type ntype or signal an error. The allowable values of ntype are float, bfloat, fixnum, rational, and any (no coercion)."
  19   (let ((type-test) (xx x))
  20        (cond
  21          ((eql ntype '$float)
  22           (setq xx ($float x))
  23           (setq type-test #'$floatnump)) ;what happens to binary32 numbers?
  24
  25          ((eql ntype '$bfloat)
  26           (setq xx ($bfloat x))
  27           (setq type-test #'$bfloatp))
  28
  29          ((eql ntype 'fixnum)
  30           (setq xx (coerce x 'fixnum))
  31           (setq type-test #'fixnump))
  32
  33          ((eql ntype '$rational)
  34           (setq xx ($rationalize x))
  35           (setq type-test #'$ratnump))
  36
  37          ((eql ntype '$any)
  38           (setq type-test #'(lambda (x) (declare (ignore x)) t)))
  39
  40          (t (merror (intl:gettext "Unknown type ~M") ntype)))
  41        (if (funcall type-test xx) xx
  42            (merror (intl:gettext "Unable to convert ~M to type ~M") x ntype))))
  43
  44 ;;; Possibly coerce-float-fun (defined in plot.lisp) is an alternative to expr-to-compiled-fun.
  45 ;;; But unlike coerce-float-fun, the function expr-to-compiled-fun calls the compiler. After
  46 ;;; compiling, return the function identifier.
  47
  48 (defun expr-to-compiled-fun (e ntype args)
  49   "Return a compiled CL function args |--> e where the args are modedeclared to be type ntype. When ntype isn't one of the types
  50   accepted by modedeclare, the type is changed to the type any. The argument args must be a CL list. Return the identifier for the function."
  51
  52   (unless (member ntype (list '$any '$boolean '$fixnum '$number '$rational '$float))
  53     (mtell  (intl:gettext  "Type ~M is unknown to the compiler; changing to type any.") ntype)
  54     (setq ntype '$any))
  55
  56   (let ((f (gensym)))
  57        (meval `((mdefine) ((,f) ,@args) ((mprogn) (($modedeclare) ((mlist) ,@args) ,ntype) ,e)))
  58        (second (mfuncall '$compile f)))) ;compile returns a Maxima list of function identifiers.
  59
  60 ;;; The Poisson bracket {f,g}.
  61 (defun $poisson_bracket (f g p q)
  62   "Return the poisson bracket {f,g} for the canonical coordinates p & q. Either both p & q must be
  63   equal length Maxima lists or both must be ."
  64   (cond
  65     ((and ($listp p) ($listp q))
  66      (setq p (rest p))
  67      (setq q (rest q)))
  68     ((and ($mapatom p) ($mapatom q))
  69      (setq p (list p))
  70      (setq q (list q)))
  71     (t
  72      (merror (intl:gettext
  73               "Either the third and fourth arguments to poisson_bracket must (i) both be lists or (ii) both be mapatoms"))))
  74
  75   (reduce #'add
  76           (mapcar #'(lambda (pk qk) (sub
  77                                      (mul ($diff f qk) ($diff g pk))
  78                                      (mul ($diff f pk) ($diff g qk)))) p q)))
  79
  80 ;;; Hash table for coefficients of symplectic methods; each value is a list of the
  81 ;;; form ((c1 d1) (c2 d2) ...)
  82
  83 (defvar *symplectic-method-coefficients* (make-hash-table :size 5 :test 'eql))
  84 (defvar *symplectic-method-digits* (make-hash-table :size 5 :test 'eql))
  85
  86 (setf (gethash '$symplectic_euler *symplectic-method-coefficients*)
  87       (list (cons 1 1)))
  88 (setf (gethash '$symplectic_euler *symplectic-method-digits*) '$inf)
  89
  90 (setf (gethash '$verlet *symplectic-method-coefficients*)
  91       (list (cons 0 (div 1 2))
  92             (cons 1 (div 1 2))))
  93 (setf (gethash '$verlet *symplectic-method-digits*) '$inf)
  94
  95 (setf (gethash '$symplectic_third_order *symplectic-method-coefficients*)
  96       (list (cons 1 (div -1 24))
  97             (cons (div -2 3) (div 3 4))
  98             (cons (div  2 3) (div 7 24))))
  99 (setf (gethash '$symplectic_third_order *symplectic-method-digits*) '$inf)
 100
 101 (setf (gethash '$symplectic_fourth_order *symplectic-method-coefficients*)
 102       (let ((k (div 1 (sub 2 (power 2 (div 1 3))))))
 103            (list
 104             (cons (div k 2) k)
 105             (cons (div (mul (sub 1 (power 2 (div 1 3))) k) 2) (mul (mul -1 (power 2 (div 1 3))) k))
 106             (cons (div (mul (sub 1 (power 2 (div 1 3))) k) 2) k)
 107             (cons (div k 2) 0))))
 108 (setf (gethash '$symplectic_fourth_order *symplectic-method-digits*) '$inf)
 109
 110 ;; see Kostas Tselios and T. E. Simos, ``Optimized fifth order symplectic integrators
 111 ;;; for orbital problems,'' Revista Mexicana de Astronom´ıa y Astrof´ısica, 49, 11–24
 112 ;;; (2013)
 113
 114 ;;; cut and paste coefficients from article--the article says something about using 40
 115 ;;; digits, but the coeffficeints are given to 57 digits. Let's be cautious and say these
 116 ;;; are good to 40 digits.
 117
 118 #|
 119 c1 = 0.112569584468347104973189684884327785393840239333314075493 ,
 120 c2 = 0.923805029000837468447500070054064432491178527428114178991 ,
 121 c3 = −1.362064898669775624786044007840908597402026042205084284026 ,
 122 c4 = 0.980926531879316517259793318227431991923428491844523669724 ,
 123 c5 = 0.400962967485371350147918025877657753577504227492190779513 ,
 124 c6 = 0.345821780864741783378055242038676806930765132085822482512 ,
 125 c7 = −0.402020995028838599420412333241250172914690575978880873429 ,
 126 d1 = 0.36953388878114957185081450061701658106775743968995046842 ,
 127 d2 = −0.032120004263046859169923904393901683486678946201463277409 ,
 128 d3 = −0.011978701020553903586622444048386301410473649207894475166 ,
 129 d4 = 0.51263817465269673604202785657395553607442158325539698102 ,
 130 d5 = −0.334948298035883491345320878224434762455516821029015086331 ,
 131 d6 = 0.021856594741098449005512783774683495267598355789295971623 ,
 132 d7 = 0.47501834514453949720351208570106713494289203770372938037
 133 |#
 134
 135 (setf (gethash '$symplectic_fifth_order *symplectic-method-coefficients*)
 136       (mapcar #'(lambda (a b) (cons a b))
 137               (list
 138                (div 5652885872144153264930671803559342469576954521714583915513
 139                     50216813883093446110686315385661331328818843555712276103168)
 140                (div 362426134418756238993387367402013039958979964654753376945
 141                     392318858461667547739736838950479151006397215279002157056)
 142                (div -4274909969574666030908365902533105271208568565239518523207
 143                     3138550867693340381917894711603833208051177722232017256448)
 144                (div 3078687817773247970859746019343925292770783479419914796411
 145                     3138550867693340381917894711603833208051177722232017256448)
 146                (div 629221334757054438917474234802547577988423100782937501865
 147                     1569275433846670190958947355801916604025588861116008628224)
 148                (div 4341517001601166158290432120166041572705763778637534302909
 149                     12554203470773361527671578846415332832204710888928069025792)
 150                (div -5047053371114805865785231327168622346738894852075304064687
 151                     12554203470773361527671578846415332832204710888928069025792))
 152               (list
 153                (div 2319601814552342659136359435206712577059674598536966859091
 154                     6277101735386680763835789423207666416102355444464034512896)
 155                (div -3225928552003184381834780354143261178475501353749927450079
 156                     100433627766186892221372630771322662657637687111424552206336)
 157                (div -4812257597683015175999923204116192762946117963505537572213
 158                     401734511064747568885490523085290650630550748445698208825344)
 159                (div 402235246967237878148042461203282093231351221007491081787
 160                     784637716923335095479473677900958302012794430558004314112)
 161                (div -4205009085731718839023168486933272699269839833464512439795
 162                     12554203470773361527671578846415332832204710888928069025792)
 163                (div 274392137557984949615301132056520558869905595208602668291
 164                     12554203470773361527671578846415332832204710888928069025792)
 165                (div 5963476957294596320417393096022684997935176483307435518905
 166                     12554203470773361527671578846415332832204710888928069025792))))
 167
 168 (setf (gethash '$symplectic_fifth_order *symplectic-method-digits*) 40)
 169
 170 ;;; Usage: symplectic_ode(ham, [p1,..., pn], [q1, ... , qn],
 171 ;;;                    [p01, ... p0n], [q01, ... q0n], dt, N, [method], [ntype])
 172
 173 ;;;   ham = time independent hamiltonian (function of p1 ... pn and q1 ... qn only)
 174 ;;;   p1 ... pn; q1 ... qn = canonical  momenta and position, respectively
 175 ;;;   p01 .. p0n; q01 ... q0n  = initial values for p1 ... pn & q1 ... qn, respectively
 176 ;;;   dt = time step--can be negative
 177 ;;;   N = number of time steps--must be a positive fixnum
 178 ;;;   method = integration method (default symplectic_euler); must be one of symplectic_euler,
 179 ;;;            verlet, symplectic_third_order, symplectic_fourth_order
 180 ;;;   ntype = optional number type (default float)--can be float, rational, or any (no type)
 181
 182 ;;; The function symplectic returns a two member list of two lists. The first list has the form
 183 ;;; [[p00,p10,p20, ...],[p01,p11, p21, ..] ...] where pij = i-th momentum at time step j.
 184
 185 ;;; For a separable hamiltonian of the form ham = F(p) + G(q), symplectic_ode approximately
 186 ;;; integrates the hamiltonian equation using a method that preserves the poisson brackets of p & q.
 187 ;;; For a description, see https://en.wikipedia.org/wiki/Symplectic_integrator#A_first-order_example.
 188 ;;; For non separable hamiltonians, the method does not in general preserve the Poisson brackets.
 189
 190 ;;; Basically, the method is
 191
 192 ;;; for k = 1 ... n
 193 ;;;   q <-  q + c(k)*diff(ham,p)*dt
 194 ;;;   p <-  p - d(k)*diff(ham,q)*dt ;use updated position q
 195
 196 ;;; For example, the symplectic Euler has c(1)=1, d(1)=1, and n=1. This code has built-in methods
 197 ;;; for the symplectic Euler, the Verlet method, and third and fourth order methods due to Ronald Ruth.
 198 ;;; A hash table *symplectic-method-coefficients* contains the coefficients for c(1), d(1), ... for
 199 ;;; each of the methods.
 200
 201 ;;; There is no user level mechanism for appending methods, but a new method can be
 202 ;;; added by including a new entry in the hash table *symplectic-method-coefficients*. The coefficients
 203 ;;; should be Maxima expressions in exact form---depending on the optional value ntype (number type)
 204 ;;; the coefficients are converted at runtime to the correct type.
 205
 206 ;;; This function creates and compiles functions for updating p & q. The arguments to these functions
 207 ;;; are modedeclared to have type ntype (default $float). Since float is a Maxima option variable,
 208 ;;; users need to remember to quote float. Of course, hand coding of these functions could increase
 209 ;;; speed or accuracy, but the automatically generated functions are convenient.
 210
 211 (defun $symplectic_ode (ham p q po qo dt NN &optional (sym-method '$symplectic_euler) (ntype '$float))
 212   "symplectic ode solver: Usage: (non scalar or scalar version)
 213        symplectic_ode(ham, [p1,..., pn], [q1, ... , qn], [p01, ... p0n],
 214              [q01, ... q0n], dt, N, [method], [ntype])
 215        symplectic_ode(ham, p, q, po, qo, dt, [method],[ntype])
 216
 217   where
 218
 219   ham = time independent hamiltonian (function of p1 ... pn and q1 ... qn only)
 220   p1 ... pn; q1 ... qn = canonical  momenta and position, respectively
 221   p01 .. p0n, q01 ... q0n  = initial values for p1 ... pn & q1 ... qn, respectively
 222   dt = time step--can be negative
 223   N = number of time steps--must be a positive fixnum
 224   method = integration method (default symplectic_euler); must be one of symplectic_euler,
 225   verlet, symplectic_third_order, symplectic_fourth_order
 226   ntype = optional number type (default float)--can be float, rational, or any (no type)"
 227
 228
 229   (let* ((update-p) (update-q) (cfs) (ddt (gensym)) (args) (poo) (qoo)
 230          (xmlist (get 'mlist 'msimpind)) (scalar-case) (N (mcoerce NN 'fixnum)))
 231
 232         (declare (type fixnum N))
 233         (cond
 234           ((every #'$listp (list p q po qo)) ;Maxima to CL list conversion
 235            (setq scalar-case nil)
 236            (setq po (list (mapcar #'(lambda (x) (mcoerce x ntype)) (rest po))))
 237            (setq qo (list (mapcar #'(lambda (x) (mcoerce x ntype)) (rest qo))))
 238            (setq p (rest p))
 239            (setq q (rest q)))
 240
 241           ((every #'(lambda (s) (not ($listp s))) (list p q po qo))
 242            (setq scalar-case t)
 243            (setq po (list (list (mcoerce po ntype))))
 244            (setq qo (list (list (mcoerce qo ntype))))
 245            (setq p (list p))
 246            (setq q (list q)))
 247
 248           (t (merror (intl:gettext "Either the second through fifth arguments must all be lists or all must be nonlists."))))
 249
 250
 251         (when (or (some #'(lambda (s) (not ($symbolp s))) p)
 252                   (some #'(lambda (s) (not ($symbolp s))) q))
 253               (merror "The second and third arguments must be either be symbols or lists of symbols."))
 254
 255         (setq args (append p q (list ddt)))
 256         (setq update-p (mapcar #'(lambda (pk qk)
 257                                          (expr-to-compiled-fun (sub pk (mult ddt ($diff ham qk))) ntype args)) p q))
 258
 259         (setq update-q (mapcar #'(lambda (pk qk)
 260                                         (expr-to-compiled-fun (add qk (mult ddt ($diff ham pk))) ntype args)) p q))
 261
 262         (setq cfs (gethash sym-method *symplectic-method-coefficients* ))
 263
 264         (when (not cfs)
 265           (let ((all-methods nil)) ;build a list of all methods
 266                (maphash #'(lambda (a b) (declare (ignore b)) (push a all-methods)) *symplectic-method-coefficients*)
 267                (push xmlist all-methods)
 268                (merror (intl:gettext "The method must be a member of ~M, but found ~M") all-methods sym-method)))
 269
 270         (setq cfs (mapcar #'(lambda (x) (cons (mcoerce (mul dt (car x)) ntype)
 271                                               (mcoerce (mul dt (cdr x)) ntype))) cfs))
 272
 273         (setq poo (car po))
 274         (setq qoo (car qo))
 275         (while (> N 0)
 276                (decf N 1)
 277                (dolist (cx cfs)
 278                  (setq args (append poo qoo (list (car cx))))
 279                  (setq poo (mapcar #'(lambda (fk) (apply fk args)) update-p))
 280                  (setq args (append poo qoo (list (cdr cx))))
 281                  (setq qoo (mapcar #'(lambda (fk) (apply fk args)) update-q)))
 282                (push poo po)
 283                (push qoo qo))
 284
 285         (cond
 286           (scalar-case
 287             (setq po (mapcar #'car (reverse po)))
 288             (setq qo (mapcar #'car (reverse qo))))
 289           (t
 290            (setq po (mapcar #'(lambda (x) (push xmlist x)) (reverse po)))
 291            (setq qo (mapcar #'(lambda (x) (push xmlist x)) (reverse qo)))))
 292
 293         (mapcar #'fmakunbound update-p) ;The functions update-p and update-q are no longer needed.
 294         (mapcar #'fmakunbound update-q)
 295         (push xmlist po) ;Convert po & qo from CL lists to Maxima lists.
 296         (push xmlist qo) ;Arguably circumventing simplification by pushing (mlist simp) is a bug.
 297         (list xmlist po qo)))