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