share/contrib/diffequations/lazysolver.lisp

   1 ;; Author Barton Willis
   2 ;; University of Nebraska at Kearney
   3 ;; Copyright (C) 2004, Barton Willis
   4
   5 ;; Brief Description: Maxima code for linear homogeneous second order
   6 ;; differential equations.
   7
   8 ;; Maxima odelin is free software; you can redistribute it and/or
   9 ;; modify it under the terms of the GNU General Public License,
  10 ;; http://www.gnu.org/copyleft/gpl.html.
  11
  12 ;; Maxima odelin has NO WARRANTY, not even the implied warranty of
  13 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  14
  15 ($put '$lazysolver 1 '$version)
  16
  17 (eval-when
  18     (:load-toplevel :compile-toplevel :execute)
  19   ($load "odeutils"))
  20
  21 (defmvar $aalgsys_is_loquacious nil)
  22
  23 (defun variablep (e)
  24   (or ($symbolp e) ($subvarp e)))
  25
  26 ;; Check that sol is a solution to eqs and that nz is nonvanishing.
  27
  28 (defun $checksolution (sol eqs &optional (nz `((mlist))))
  29   (let (($gcd '$spmod) ($algebraic t) ($ratfac nil) ($ratprint nil))
  30     (if (not ($listp sol)) (setq sol `((mlist) ,sol)))
  31     (if (not ($listp eqs)) (setq eqs `((mlist) ,eqs)))
  32     (if (not ($listp nz)) (setq nz `((mlist) ,nz)))
  33     (setq eqs (mapcar #'meqhk (cdr eqs)))
  34     (setq eqs `((mlist) ,@eqs))
  35     (and
  36      (every #'mequalp (cdr sol))
  37      (every #'(lambda (s) (variablep ($lhs s))) (cdr sol))
  38      (every #'(lambda (s) (like 0 s)) (cdr (sratsimp ($substitute sol eqs))))
  39      (every #'(lambda (s) (not (like 0 s)))
  40             (cdr (sratsimp ($substitute sol nz)))))))
  41
  42 ;; The function 'checkedalgsys' tries to return the "simplest" solution.
  43 ;; it uses this simple-minded measure of simple.
  44
  45 (defun my-expr-size (e)
  46   (xmy-expr-size ($totaldisrep e)))
  47
  48 (defun xmy-expr-size (e)
  49   (if (consp e) (apply #'+ (mapcar #'my-expr-size (margs e))) 1))
  50
  51 (defun my-freeof (unks q)
  52   (every #'(lambda (x) ($freeof x q)) (cdr unks)))
  53
  54 ;; This function solves 'eqs' for 'vars' and returns a solution such that
  55 ;; no expression in the Maxima list 'nz' vanishes. Since
  56 ;; 'algsys' sometimes returns bogus solutions--checkedalgsys checks
  57 ;; the putative solutions and rejects bogus solutions.  To return
  58 ;; the 'simpliest' solution, we sort the putative solutions using
  59 ;; 'my-expr-size.
  60
  61 (defun $checkedalgsys (eqs vars &optional (nz `((mlist))))
  62
  63   (let ((sol) ($ratfac nil) ($ratprint nil) ($realonly nil) ($algexact t)
  64         ($gcd '$spmod) ($algebraic t)
  65         ($programmode t) ($globalsolve nil) ($solveexplicit t)
  66         ($listconstvars t) ($solveradcan nil) ($ratvars nil))
  67
  68     (if $aalgsys_is_loquacious
  69         (mtell "...solving ~:M equations in ~:M variables~%"
  70                ($length eqs) ($length vars)))
  71
  72     (setq sol (cdr ($algsys eqs vars)))
  73     (setq sol (sort sol #'(lambda (a b) (< (my-expr-size a) (my-expr-size b)))))
  74     (dolist (si sol)
  75       (if ($checksolution si eqs nz) (return si)))))
  76
  77 ;; The function 'checkedalgsys' tries to return the "simplest" solution.
  78 ;; it uses this simple-minded measure of simple.
  79
  80 (defun nonconstant-factors (e vars)
  81   (let (acc)
  82     (setq e ($factor e))
  83     (setq e (if (mtimesp e) (margs e) (list e)))
  84     (dolist (ei e `(($set) ,@acc))
  85       (if (mexptp ei) (setq ei (car (margs ei))))
  86       (if (not (my-freeof vars ei)) (push ei acc)))))
  87
  88 (defun variablep (e)
  89   (or ($symbolp e) ($subvarp e)))
  90
  91 (defun $checksolution (sol eqs &optional (nz `((mlist))))
  92   (let (($gcd '$spmod) ($algebraic t) ($ratfac nil) ($ratprint nil))
  93     (if (not ($listp sol)) (setq sol `((mlist) ,sol)))
  94     (if (not ($listp eqs)) (setq eqs `((mlist) ,eqs)))
  95     (if (not ($listp nz)) (setq nz `((mlist) ,nz)))
  96     (setq eqs (mapcar #'meqhk (cdr eqs)))
  97     (setq eqs `((mlist) ,@eqs))
  98     (and
  99      (every #'mequalp (cdr sol))
 100      (every #'(lambda (s) (variablep ($lhs s))) (cdr sol))
 101      (every #'(lambda (s) (like 0 s)) (cdr (sratsimp ($substitute sol eqs))))
 102      (every #'(lambda (s) (not (like 0 s)))
 103             (cdr (sratsimp ($substitute sol nz)))))))
 104
 105 (defun unks-in-eq (eq unks)
 106   (let (($listconstvars nil))
 107     ($intersection ($setify ($listofvars eq)) unks)))
 108
 109 (defun $aalgsys (e-eqs eqs unks &optional (nz `((mlist))))
 110   (let ((e-vars) (sol) (e-sol) ($gcd '$spmod) ($algebraic t)
 111         ($ratvars nil) ($radexpand nil) ($ratfac nil) ($ratprint nil))
 112
 113     (setq unks ($setify unks))
 114     (setq unks (mbag-map #'$ratdisrep unks))
 115
 116     (setq eqs ($setify eqs))
 117     (setq eqs (mbag-map #'$ratdisrep eqs))
 118
 119     (setq e-eqs ($setify e-eqs))
 120     (setq e-eqs (mbag-map #'$ratdisrep e-eqs))
 121     (setq e-eqs (mbag-map #'$factor e-eqs))
 122
 123     (setq e-eqs ($union e-eqs
 124                         ($subset eqs #'(lambda (w)
 125                                          (= 1 (number-of-unks w unks))))))
 126
 127     (setq e-vars (unks-in-eq e-eqs unks))
 128     (setq e-eqs ($union e-eqs
 129                         ($subset eqs #'(lambda (w)
 130                                          (like (unks-in-eq w unks) e-vars)))))
 131
 132     (setq e-eqs (mbag-map #'$ratdisrep e-eqs))
 133     (setq e-eqs (mbag-map #'$factor e-eqs))
 134     (setq e-eqs ($disjoin 0 e-eqs))
 135
 136    ;; (displa `((mequal) eeqs ,e-eqs))
 137     (setq e-eqs ($listify e-eqs))
 138     (setq e-vars ($listify e-vars))
 139     (setq eqs ($listify eqs))
 140     (setq unks ($listify unks))
 141
 142     (block bailout
 143       (cond ((not ($emptyp e-vars))
 144              (setq e-sol ($algsys e-eqs e-vars))
 145              (setq e-sol ($setify e-sol))
 146              (setq e-sol
 147                    ($subset e-sol #'(lambda (w) ($checksolution w e-eqs nz))))
 148              (setq e-sol
 149                    ($subset e-sol #'(lambda (w) (my-freeof $%rnum_list w))))
 150              (setq e-sol ($listify e-sol))
 151
 152              (cond ((not ($emptyp e-sol))
 153                     (setq e-sol (margs e-sol))
 154                     (dolist (ei e-sol)
 155                      ;; (displa `((mequal) auxeq ,ei))
 156                       (setq sol ($checkedalgsys
 157                                  ($append ei ($substitute ei eqs)) unks nz))
 158                       (if ($listp sol) (return-from bailout sol))))
 159                    (t ($checkedalgsys ($append e-eqs eqs) unks nz))))
 160             (t (return-from bailout ($checkedalgsys eqs unks nz)))))))
 161
 162
 163
 164
 165
 166
 167
 168
 169
 170
 171
 172
 173
 174
 175
 176
 177
 178
 179
 180
 181
 182