share/contrib/diffequations/sqfr.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 '$sqfr 1 '$version)
  16
  17 (eval-when
  18     (:load-toplevel :compile-toplevel :execute)
  19   ($load "odeutils"))
  20
  21 ;; If x is a symbol for a subvarp, return its general representation.
  22 ;; Otherwise signal an error---the argument f is the string name of
  23 ;; a function that requires the symbol.
  24
  25 (defun require-symbol (x f &optional a)
  26   (declare (ignore a))
  27   (setq x ($ratdisrep x))
  28   (if (or (symbolp x) ($subvarp x)) x
  29     (merror "Function ~:M requires a symbol; instead found ~:M" f x))
  30   x)
  31
  32 (defun $strictmysqfr (p x)
  33   (setq p ($expand ($ratdisrep p)))
  34   (setq x (require-symbol x "$mysqfr"))
  35   (let ((i 1) (lc 1) (r) (ai) (w) (bad) (acc nil) (q) (s)
  36         ($gcd '$spmod) ($algebraic t) ($resultant '$red) ($ratfac nil)
  37         ($ratprint nil))
  38
  39     (setq lc ($coeff p x ($hipow p x)))
  40     (setq p ($expand (div p lc)))
  41     (while (not ($freeof x p))
  42       (setq r ($gcd p ($diff p x) x))
  43       (setq q ($first ($divide p r)))
  44       (setq ai ($first ($divide q ($gcd q r x) x)))
  45       (setq ai (div ai ($coeff ai x ($hipow ai x))))
  46       (push ai acc)
  47       (setq p r)
  48       (incf i))
  49
  50     (setq acc (reverse acc))
  51     (setq r lc)
  52
  53     (setq i 0)
  54     (while acc
  55       (setq ai (pop acc))
  56       (setq r (mul r (power ai (incf i))))
  57       (setq s ($resultant ai ($diff ai x) x))
  58       (if (not ($constantp s)) (push ($factor s) bad))
  59       (dolist (wi w)
  60         (setq s ($resultant ai wi x))
  61         (if (not ($constantp s)) (push ($factor s) bad)))
  62       (push ai w))
  63
  64     (setq bad `(($set) ,@bad))
  65     (setq bad (mbag-map #'(lambda (w) `((mnotequal) ,w 0)) bad))
  66     (if (not ($emptyp bad)) (mtell "Proviso: assuming ~:M~%" bad))
  67     (values r bad)))
  68
  69 (defun $mysqfr (p x)
  70   (setq p ($expand ($ratdisrep p)))
  71   (setq x (require-symbol x "$mysqfr"))
  72   (let ((i 0) (r) (ai) (acc) (q) ($gcd '$spmod) ($algebraic t))
  73
  74     (cond ((like 0 p) 0)
  75           (t
  76            (setq p ($expand p))
  77            (setq acc ($coeff p x ($hipow p x)))
  78            (setq p (div p acc))
  79            (while (not ($freeof x p))
  80              (setq r ($gcd p ($diff p x) x))
  81              (setq q ($first ($divide p r x)))
  82              (setq ai ($first ($divide q ($gcd q r x) x)))
  83              (setq ai (div ai ($coeff ai x ($hipow ai x))))
  84              (setq acc (mul acc (power ai (incf i))))
  85              (setq p r))
  86            acc))))
  87