share/dynamics/plotdf.lisp

   1 ;; plotdf.mac - Adds a function plotdf() to Maxima, which draws a Direction
   2 ;;              Field for an ordinary 1st order differential equation,
   3 ;;              or for a system of two autonomous 1st order equations.
   4 ;;
   5 ;; Copyright (C) 2004, 2008, 2011 Jaime E. Villate <villate@fe.up.pt>
   6 ;;
   7 ;; This program is free software; you can redistribute it and/or modify
   8 ;; it under the terms of the GNU General Public License as published by
   9 ;; the Free Software Foundation; either version 2 of the License, or
  10 ;; (at your option) any later version.
  11 ;;
  12 ;; This program is distributed in the hope that it will be useful,
  13 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 ;; GNU General Public License for more details.
  16 ;;
  17 ;; You should have received a copy of the GNU General Public License
  18 ;; along with this program; if not, write to the Free Software
  19 ;; Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
  20 ;; MA  02110-1301, USA.
  21 ;;
  22 ;; See plotdf.usg (which should come together with this program) for
  23 ;; a usage summary
  24 ;;
  25 ;; $Id: plotdf.lisp,v 1.12 2011-03-09 11:33:46 villate Exp $
  26
  27 (in-package :maxima)
  28
  29 ;; parses a plotdf option into a command-line option for tcl scripts
  30 (defun plotdf-option-to-tcl (value s1 s2)
  31   (let (v vv)
  32     (unless (and  ($listp value)
  33                   (symbolp (setq name (second value))))
  34       (merror
  35        (intl:gettext
  36         "plotdf-option-to-tcl: Expecting a symbol for the option name, found: \"~M\"") value))
  37     (case name
  38       (($xradius $yradius $xcenter $ycenter $tinitial $tstep)
  39        (setq v (check-option (cdr value) #'realp "a real number" 1))
  40        (setq value (list '(mlist) name v)))
  41       (($width $height $nsteps $versus_t)
  42        (setq v (check-option (cdr value) #'naturalp "a natural number" 1))
  43        (setq value (list '(mlist) name v)))
  44       ($trajectory_at
  45        (setq v (check-option (cdr value) #'realp "a real number" 2))
  46        (setq value (cons '(mlist) (cons name v))))
  47       ($bbox
  48        (setq v (check-option (cdr value) #'realp "a real number" 4))
  49        (setq value (cons '(mlist) (cons name v))))
  50       (($xfun $parameters $sliders $vectors $fieldlines $curves
  51               $windowtitle $xaxislabel $yaxislabel $psfile)
  52        (setq v (check-option (cdr value) #'stringp "a string" 1))
  53        (setq value (list '(mlist) name v)))
  54       ($axes
  55        (if (not (third value))
  56            (setq value '((mlist) $axes 0))
  57            (case (third value)
  58              ($x (setq value '((mlist) $axes "x")))
  59              ($y (setq value '((mlist) $axes "y")))
  60              (t (setq value '((mlist) $axes "xy"))))))
  61       ($box
  62        (if (not (third value))
  63            (setq value '((mlist) $nobox 1))
  64            (setq value '((mlist) $nobox 0))))
  65       ($direction
  66        (or
  67         (member (ensure-string (third value)) '("forward" "backward" "both") :test #'equal)
  68         (merror
  69          (intl:gettext
  70           "plotdf-option-to-tcl: direction should be forward, backward or both."))))
  71       (t (cond
  72            ((eql name s1)
  73             (setq value (cons '(mlist) (cons '$x (cddr (check-range value))))))
  74            ((eql name s2)
  75             (setq value (cons '(mlist) (cons '$y (cddr (check-range value))))))
  76            (t (merror (intl:gettext "plotdf-option-to-tcl: unknown option ~M") name)))))
  77     (setq vv (mapcar #'(lambda (a) (if (symbolp a) (ensure-string a) a)) (cdr value)))
  78     (with-output-to-string (st)
  79       (cond ((or (equal (first vv) "x") (equal (first vv) "y"))
  80              (format st "-~(~a~)center " (first vv))
  81              (format st "{~a} " (/ (+ (third vv) (second vv)) 2))
  82              (format st "-~(~a~)radius " (first vv))
  83              (format st "{~a}" (/ (- (third vv) (second vv)) 2)))
  84             (t
  85              (format st "-~(~a~) " (first vv))
  86              (format st "{~{~a~^ ~}}" (rest vv)))))))
  87
  88 ;; applies float(ev(expression, numer)) to an expression, and returns a string
  89
  90 (defun expr_to_str (fun)
  91   (mstring (mfuncall '$float (mfuncall '$ev fun '$numer))))
  92
  93 ;; plots the direction field for an ODE  dy/dx = f(x,y), or for an autonomous
  94 ;; system of 2 equations dx/dt = f(x,y), dy/dt = g(x,y)
  95 ;;
  96 (defun $plotdf (ode &rest options)
  97   (let (file cmd (opts " ") (s1 '$x) (s2 '$y) vars)
  98     (unless ($listp ode) (setf ode `((mlist) ,ode)))
  99
 100     ;; if the variables are not x and y, their names must be given right
 101     ;; after the expression for the ode's
 102     (when
 103       (and (listp (car options)) (= (length (car options)) 3)
 104             (or (symbolp (cadar options)) ($subvarp (cadar options)))
 105             (or (symbolp (caddar options)) ($subvarp (caddar options))))
 106       (setq s1 (cadar options))
 107       (setq s2 (caddar options))
 108       (setq options (cdr options)))
 109
 110     ;; parse options and copy them to string opts
 111     (cond (options
 112            (dolist (v options)
 113              (setq opts (concatenate 'string opts " "
 114                                   (plotdf-option-to-tcl v s1 s2))))))
 115
 116     (unless (search "-xaxislabel " opts)
 117       (setq opts (concatenate 'string opts " -xaxislabel " (ensure-string s1))))
 118     (unless (search "-yaxislabel " opts)
 119       (setq opts (concatenate 'string opts " -yaxislabel " (ensure-string s2))))
 120
 121     ;; check that the expressions given contain only the axis variables
 122     (setq vars (cdr (mfuncall '$listofvars ode)))
 123     (unless (search "-parameters " opts)
 124       (dolist (var vars)
 125         (unless (or (eq var s1) (eq var s2))
 126           (merror
 127            (intl:gettext "plotdf: expression(s) given can only depend on ~M and ~M~%Found extra variable ~M") s1 s2 var))))
 128
 129     ;; substitute $x by s1 and $y by s2
 130     (defun subxy (expr)
 131       (if (listp expr)
 132           (mapcar #'subxy expr)
 133         (cond ((eq expr s1) '$x) ((eq expr s2) '$y) (t expr))))
 134     (setf ode (mapcar #'subxy ode))
 135
 136     ;; parse the differential equation expressions
 137     (case (length ode)
 138           (3 (setq cmd (concatenate 'string " -dxdt \""
 139                                     (expr_to_str (second ode)) "\" -dydt \""
 140                                     (expr_to_str (third ode)) "\"")))
 141           (2 (setq cmd (concatenate 'string " -dydx \""
 142                                     (expr_to_str (second ode)) "\"")))
 143           (t (merror
 144               (intl:gettext "plotdf: first argument must be either an expression or a list with two expressions."))))
 145     (setq file (plot-temp-file (format nil "maxout~d.xmaxima" (getpid))))
 146     (show-open-plot
 147      (with-output-to-string
 148          (st)
 149        (cond
 150          ($show_openplot (format st "plotdf ~a ~a~%" cmd opts))
 151          (t (format st "{plotdf ~a ~a} " cmd opts))))
 152      file)
 153     file))
 154
 155 ;; plot equipotential curves for a scalar field f(x,y)
 156 (defun $ploteq (fun &rest options)
 157
 158   (let (file cmd mfun (opts " ") (s1 '$x) (s2 '$y))
 159     (setf mfun `((mtimes) -1 ,fun))
 160     ;; if the variables are not x and y, their names must be given right
 161     ;; after the expression for the ode's
 162     (when
 163       (and (listp (car options)) (= (length (car options)) 3)
 164             (or (symbolp (cadar options)) ($subvarp (cadar options)))
 165             (or (symbolp (caddar options)) ($subvarp (caddar options))))
 166       (setq s1 (cadar options))
 167       (setq s2 (caddar options))
 168       (setq options (cdr options)))
 169     (defun subxy (expr)
 170       (if (listp expr)
 171           (mapcar #'subxy expr)
 172           (cond ((eq expr s1) '$x) ((eq expr s2) '$y) (t expr))))
 173     (setf mfun (mapcar #'subxy mfun))
 174 ;; the next two lines should take into account parameters given in the options
 175 ;;    (if (delete '$y (delete '$x (rest (mfuncall '$listofvars ode))))
 176 ;;        (merror "The equation(s) can depend only on 2 variable which must be specified!"))
 177     (setq cmd (concatenate 'string " -dxdt \""
 178                            (expr_to_str (mfuncall '$diff mfun '$x))
 179                            "\" -dydt \""
 180                            (expr_to_str (mfuncall '$diff mfun '$y))
 181                            "\" "))
 182
 183     ;; parse options and copy them to string opts
 184     (cond (options
 185            (dolist (v options)
 186              (setq opts (concatenate 'string opts " "
 187                                   (plotdf-option-to-tcl v s1 s2))))))
 188
 189     (unless (search "-vectors " opts)
 190       (setq opts (concatenate 'string opts " -vectors {}")))
 191     (unless (search "-fieldlines " opts)
 192       (setq opts (concatenate 'string opts " -fieldlines {}")))
 193     (unless (search "-curves " opts)
 194       (setq opts (concatenate 'string opts " -curves {red}")))
 195     (unless (search "-windowtitle " opts)
 196       (setq opts (concatenate 'string opts " -windowtitle {Ploteq}")))
 197     (unless (search "-xaxislabel " opts)
 198       (setq opts (concatenate 'string opts " -xaxislabel " (ensure-string s1))))
 199     (unless (search "-yaxislabel " opts)
 200       (setq opts (concatenate 'string opts " -yaxislabel " (ensure-string s2))))
 201
 202     (setq file (plot-temp-file (format nil "maxout~d.xmaxima" (getpid))))
 203     (show-open-plot
 204      (with-output-to-string
 205          (st)
 206        (cond ($show_openplot (format st "plotdf ~a ~a~%" cmd opts))
 207              (t (format st "{plotdf ~a ~a}" cmd opts))))
 208      file)
 209     file))