share/dynamics/complex_dynamics.lisp

   1 ;; complex_dynamics.lisp - functions julia, mandelbrot and rk
   2 ;;
   3 ;; Copyright (C) 2006-2021 Jaime E. Villate <villate@fe.up.pt>
   4 ;;
   5 ;; This program is free software; you can redistribute it and/or modify
   6 ;; it under the terms of the GNU General Public License as published by
   7 ;; the Free Software Foundation; either version 2 of the License, or
   8 ;; (at your option) any later version.
   9 ;;
  10 ;; This program is distributed in the hope that it will be useful,
  11 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 ;; GNU General Public License for more details.
  14 ;;
  15 ;; You should have received a copy of the GNU General Public License
  16 ;; along with this program; if not, write to the Free Software
  17 ;; Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
  18 ;; MA 02110-1301, USA.
  19 ;;
  20
  21 (in-package :maxima)
  22 (macsyma-module dynamics)
  23
  24 ;; Function $mandelbrot displays the Mandelbrot set
  25 ;; on the region: xmin < x <xmax, ymin < y <ymax.
  26
  27 (defmfun $mandelbrot (&rest extra-options)
  28   (let (plot output-file (options (copy-tree *plot-options*))
  29         a b c d e num  dx dy x xmax xmin y ymax ymin m nx ny)
  30     (setf (getf options '$type) "plot2d")
  31     (unless (getf options '$x) (setf (getf options '$x) '(-2 2)))
  32     (unless (getf options '$y) (setf (getf options '$y) '(-2 2)))
  33     (unless (getf options '$iterations) (setf (getf options '$iterations) 9))
  34     (unless (getf options '$xlabel) (setf (getf options '$xlabel) "x"))
  35     (unless (getf options '$ylabel) (setf (getf options '$ylabel) "y"))
  36     (unless (member '$color_bar options) (setf (getf options '$color_bar) t))
  37     (setf (getf options '$palette)
  38           '(((mlist) $gradient $magenta $violet $blue $cyan $green $yellow
  39             $orange $red $brown $black)))
  40
  41     ;; Parses the options given in the command line
  42     (setq options (plot-options-parser extra-options options))
  43     (unless (getf options '$yx_ratio) (setf (getf options '$same_xy) t))
  44     (setq xmin (car (getf options '$x)))
  45     (setq xmax (cadr (getf options '$x)))
  46     (setq ymin (car (getf options '$y)))
  47     (setq ymax (cadr (getf options '$y)))
  48     (setq m (getf options '$iterations))
  49     (setq nx (or (car (getf options '$grid)) 30))
  50     (setq ny (or (cadr (getf options '$grid)) 30))
  51     (setq dx (/ (rationalize (- xmax xmin)) nx)) ; x incr. per pixel
  52     (setq dy (/ (rationalize (- ymax ymin)) ny)) ; y incr. per pixel
  53
  54     ;; Creates the object that will be passed to the external graphic program
  55     (setq plot (make-instance 'gnuplot-plot))
  56     (if (eq (getf options '$plot_format) '$gnuplot_pipes)
  57         (setf (slot-value plot 'pipe) T)
  58         (setf (getf options '$plot_format) '$gnuplot))
  59
  60     (setq output-file (plot-preamble plot options))
  61     (setf
  62      (slot-value plot 'data)
  63      (concatenate
  64       'string
  65       (slot-value plot 'data)
  66       (with-output-to-string (st)
  67         (format st "set palette ~a~%"
  68                 (gnuplot-palette (rest (first (getf options '$palette)))))
  69         (format st "unset key~%")
  70         (format st "plot '-' with image~%")
  71         ;; iterates through all grid points
  72         (dotimes (i ny)
  73           (setq y (+ ymin (/ dy 2) (* i dy)))
  74           (dotimes (j nx)
  75             (setq x (+ xmin (/ dx 2) (* j dx)))
  76             (setq a 0) (setq b 0)
  77             (setq num m)
  78             (dotimes (l m)
  79               (setq c (* a a))
  80               (setq d (* b b))
  81               (setq e (* 2 a b))
  82               (when (> (+ c d) 4) (progn (setq num l) (return)))
  83               (setq a (+ (float x) (- c d)))
  84               (setq b (+ (float y) e)))
  85             (format st "~f ~f ~d~%" x y num)))
  86         (format st "e~%"))))
  87     (plot-shipout plot options output-file)))
  88
  89 ;; Function $julia(x,y) displays the Julia set for the
  90 ;; point (x,y) of the complex plane, on the region: xmin < x <xmax,
  91 ;; ymin < y <ymax.
  92
  93 (defmfun $julia (x y &rest extra-options)
  94   (let (plot output-file (options (copy-tree *plot-options*))
  95         num dx dy xmax xmin ymax ymin a b c d e a0 b0 m nx ny)
  96     (setf (getf options '$type) "plot2d")
  97     (unless (getf options '$x) (setf (getf options '$x) '(-2 2)))
  98     (unless (getf options '$y) (setf (getf options '$y) '(-2 2)))
  99     (unless (getf options '$iterations) (setf (getf options '$iterations) 9))
 100     (unless (getf options '$xlabel) (setf (getf options '$xlabel) "x"))
 101     (unless (getf options '$ylabel) (setf (getf options '$ylabel) "y"))
 102     (unless (member '$color_bar options) (setf (getf options '$color_bar) t))
 103     (setf (getf options '$palette)
 104           '(((mlist) $gradient $magenta $violet $blue $cyan $green $yellow
 105             $orange $red $brown $black)))
 106
 107     ;; Parses the options given in the command line
 108     (setq options (plot-options-parser extra-options options))
 109     (unless (getf options '$yx_ratio) (setf (getf options '$same_xy) t))
 110     (setq xmin (car (getf options '$x)))
 111     (setq xmax (cadr (getf options '$x)))
 112     (setq ymin (car (getf options '$y)))
 113     (setq ymax (cadr (getf options '$y)))
 114     (setq m (getf options '$iterations))
 115     (setq nx (or (car (getf options '$grid)) 30))
 116     (setq ny (or (cadr (getf options '$grid)) 30))
 117     (setq dx (/ (rationalize (- xmax xmin)) nx)) ; x incr. per pixel
 118     (setq dy (/ (rationalize (- ymax ymin)) ny)) ; y incr. per pixel
 119
 120     ;; Creates the object that will be passed to the external graphic program
 121     (setq plot (make-instance 'gnuplot-plot))
 122     (if (eq (getf options '$plot_format) '$gnuplot_pipes)
 123         (setf (slot-value plot 'pipe) T)
 124         (setf (getf options '$plot_format) '$gnuplot))
 125
 126     (setq output-file (plot-preamble plot options))
 127     (setf
 128      (slot-value plot 'data)
 129      (concatenate
 130       'string
 131       (slot-value plot 'data)
 132       (with-output-to-string (st)
 133         (format st "set palette ~a~%"
 134                 (gnuplot-palette (rest (first (getf options '$palette)))))
 135         (format st "unset key~%")
 136         (format st "plot '-' with image~%")
 137
 138         ;; iterates through all grid points
 139         (dotimes (i ny)
 140           (setq b0 (+ ymin (/ dy 2) (* i dy)))
 141           (dotimes (j nx)
 142             (setq a0 (+ xmin (/ dx 2) (* j dx)))
 143             (setq a (float a0))
 144             (setq b (float b0))
 145             (setq num m)
 146             (dotimes (l m)
 147               (setq c (* a a))
 148               (setq d (* b b))
 149               (setq e (* 2 a b))
 150               (when (> (+ c d) 4) (progn (setq num l) (return)))
 151               (setq a (+ x (- c d)))
 152               (setq b (+ y e)))
 153             (format st "~f ~f ~d~%" a0 b0 num)))
 154         (format st "e~%"))))
 155     (plot-shipout plot options output-file)))
 156
 157 ;; Function $rk implements the 4th order Runge-Kutta numerical method.
 158 ;;  exprs:   an expression or maxima list with n expressions
 159 ;;  vars:    name of (or list of names of) the independent variable(s)
 160 ;;  initial: initial value for the variable (or list of initial values)
 161 ;;  domain:  maxima list with four elements (name of the independent
 162 ;;           variable, its initial and final values and its increment)
 163 (defmfun $rk (exprs vars initial domain
 164             &aux d u fun k1 k2 k3 k4 r1 r2 r3 traj r
 165               (it (mapcar #'coerce-float (cddr domain))))
 166   (unless ($listp exprs) (setq exprs `((mlist) ,exprs)))
 167   (unless ($listp initial) (setq initial `((mlist) ,initial)))
 168   (unless ($listp vars) (setq vars `((mlist) ,vars)))
 169   (dolist (var (cdr vars))
 170     (unless (symbolp var)
 171       (merror (intl:gettext "rk: variable name expected; found: ~M") var)))
 172   (unless (symbolp (cadr domain))
 173     (merror (intl:gettext "rk: variable name expected; found: ~M")
 174             (cadr domain)))
 175   (setq vars (concatenate 'list '((mlist)) (list (cadr domain)) (cdr vars)))
 176   (setq r (concatenate 'list `(,(car it)) (mapcar #'coerce-float (cdr initial))))
 177   (setq fun (mapcar #'(lambda (x) (coerce-float-fun x vars)) (cdr exprs)))
 178   (setq d (/ (- (cadr it) (car it)) (caddr it)))
 179   (setq traj (list (cons '(mlist) r)))
 180   (do ((m 1 (1+ m))) ((> m d))
 181     (ignore-errors
 182       (setq k1 (mapcar #'(lambda (x) (apply x r)) fun))
 183       (setq r1 (map 'list #'+ (cdr r) (mapcar #'(lambda (x) (* (/ (caddr it) 2) x)) k1)))
 184       (push (+ (car r) (/ (caddr it) 2)) r1)
 185       (setq k2 (mapcar #'(lambda (x) (apply x r1)) fun))
 186       (setq r2 (map 'list #'+ (cdr r) (mapcar #'(lambda (x) (* (/ (caddr it) 2) x)) k2)))
 187       (push (+ (car r) (/ (caddr it) 2)) r2)
 188       (setq k3 (mapcar #'(lambda (x) (apply x r2)) fun))
 189       (setq r3 (map 'list #'+ (cdr r) (mapcar #'(lambda (x) (* (caddr it) x)) k3)))
 190       (push (+ (car r) (caddr it)) r3)
 191       (setq k4 (mapcar #'(lambda (x) (apply x r3)) fun))
 192       (setq u (map 'list #'+
 193                    (mapcar #'(lambda (x) (* 1/6 x)) k1)
 194                    (mapcar #'(lambda (x) (* 1/3 x)) k2)
 195                    (mapcar #'(lambda (x) (* 1/3 x)) k3)
 196                    (mapcar #'(lambda (x) (* 1/6 x)) k4)))
 197       (setq r
 198             (concatenate 'list
 199                          `(,(+ (car it) (* m (caddr it))))
 200                          (map 'list #'+ (cdr r) (mapcar #'(lambda (x) (* (caddr it) x)) u))))
 201       (push (cons '(mlist) r) traj)))
 202   (when (< (car r) (cadr it))
 203     (let ((s (- (cadr it) (car r))))
 204       (ignore-errors
 205         (setq k1 (mapcar #'(lambda (x) (apply x r)) fun))
 206         (setq r1 (map 'list #'+ (cdr r) (mapcar #'(lambda (x) (* (/ s 2) x)) k1)))
 207         (push (+ (car r) (/ s 2)) r1)
 208         (setq k2 (mapcar #'(lambda (x) (apply x r1)) fun))
 209         (setq r2 (map 'list #'+ (cdr r) (mapcar #'(lambda (x) (* (/ s 2) x)) k2)))
 210         (push (+ (car r) (/ s 2)) r2)
 211         (setq k3 (mapcar #'(lambda (x) (apply x r2)) fun))
 212         (setq r3 (map 'list #'+ (cdr r) (mapcar #'(lambda (x) (* s x)) k3)))
 213         (push (+ (car r) s) r3)
 214         (setq k4 (mapcar #'(lambda (x) (apply x r3)) fun))
 215         (setq u (map 'list #'+
 216                      (mapcar #'(lambda (x) (* 1/6 x)) k1)
 217                      (mapcar #'(lambda (x) (* 1/3 x)) k2)
 218                      (mapcar #'(lambda (x) (* 1/3 x)) k3)
 219                      (mapcar #'(lambda (x) (* 1/6 x)) k4)))
 220         (setq r
 221               (concatenate 'list
 222                            `(,(cadr it))
 223                            (map 'list #'+ (cdr r) (mapcar #'(lambda (x) (* s x)) u))))
 224         (push (cons '(mlist) r) traj))))
 225   (cons '(mlist) (nreverse traj)))
 226