share/diffequations/pwilt.lisp

   1 ;;
   2 ;; Copyright (C) 2010, 2011 Mark H. Weaver <mhw@netris.org>
   3 ;;
   4 ;; pwilt: inverse laplace transforms for piecewise functions
   5 ;;
   6 ;; This program is free software; you can redistribute it and/or
   7 ;; modify it under the terms of the GNU General Public License
   8 ;; as published by the Free Software Foundation; either version 2
   9 ;; of the License, or (at your option) any later version.
  10 ;;
  11 ;; This program is distributed in the hope that it will be useful,
  12 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 ;; GNU General Public License for more details.
  15 ;;
  16 ;; You should have received a copy of the GNU General Public License
  17 ;; along with this program; if not, write to the Free Software
  18 ;; Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  19 ;;
  20
  21 (in-package :maxima)
  22
  23 ($put '$pwilt 1 '$version)
  24
  25 (if (not ($get '$hstep '$version)) ($load '$hstep))
  26
  27 (defmvar $pwilt_hstep_all nil
  28   "If true, pwilt includes hstep in all terms of the result.")
  29
  30 (defun pwilt-add-parts (a b)
  31   (flet ((merge (ap bp) (cons (car ap) (add (cdr ap) (cdr bp)))))
  32     (let ((r '()))
  33       (loop (cond ((endp a) (return (nreconc r b)))
  34                   ((endp b) (return (nreconc r a)))
  35                   ((alike1 (caar a) (caar b))
  36                    (push (merge (pop a) (pop b)) r))
  37                   ((great (caar a) (caar b))
  38                    (push (pop a) r))
  39                   (t
  40                    (push (pop b) r)))))))
  41
  42 (defun pwilt-mult-parts (a b)
  43   (reduce #'pwilt-add-parts
  44           (mapcar #'(lambda (ap)
  45                       (mapcar #'(lambda (bp)
  46                                   (cons (add (car ap) (car bp))
  47                                         (mul (cdr ap) (cdr bp))))
  48                               b))
  49                   a)))
  50
  51 (defun pwilt-parts (e s)
  52   (setq e (specrepcheck e))
  53   (cond ((freeof s e) (list (cons 0 e)))
  54         ((mplusp e) (reduce #'pwilt-add-parts
  55                             (mapcar #'(lambda (e) (pwilt-parts e s))
  56                                     (margs e))))
  57         ((mtimesp e) (reduce #'pwilt-mult-parts
  58                              (mapcar #'(lambda (e) (pwilt-parts e s))
  59                                      (margs e))))
  60         ((and (mexptp e) (eq (mexpt-base e) '$%e))
  61          (let ((l (islinear (mexpt-expt e) s)))
  62            (if l (list (cons (car l) (power '$%e (cdr l))))
  63              (list (cons 0 e)))))
  64         ((and (mexptp e) (eql (mexpt-expt e) -1))
  65          (let ((dp (pwilt-parts (mexpt-base e) s)))
  66            (cond ((= 1 (length dp))
  67                   (list (cons (neg (caar dp))
  68                               (inv (cdar dp)))))
  69                  ((and (= 2 (length dp))
  70                        (freeof '%pwilt_periodic (car (first dp)))
  71                        (freeof '%pwilt_periodic (car (second dp)))
  72                        (zerop1 (add (cdr (first dp))
  73                                     (cdr (second dp)))))
  74                   (let* ((offset (car (first dp)))
  75                          (period (sub offset (car (second dp))))
  76                          (scale (cdr (first dp))))
  77                     (case (asksign period)
  78                       ($positive)
  79                       ($negative (setq offset (sub offset period)
  80                                        period (neg period)
  81                                        scale (neg scale)))
  82                       (t (return-from pwilt-parts (list (cons 0 e)))))
  83                     (list (cons (add offset (take '(%pwilt_periodic) period))
  84                                 scale))))
  85                  (t (list (cons 0 e))))))
  86         (t (list (cons 0 e)))))
  87
  88 (defun pwilt-one-part (p s x xtemp)
  89   (let ((xnew (add x (car p)))
  90         (lexp (cdr p)))
  91     (cond ((freeof s lexp) (mul lexp (take '($delta) xnew)))
  92           ((or $pwilt_hstep_all (not (zerop1 (car p))))
  93            (mul ($hstep xnew)
  94                 (maxima-substitute xnew xtemp
  95                                    (pwilt1 lexp s xtemp))))
  96           (t (pwilt1 lexp s x)))))
  97
  98 (defun pwilt-subst-periodic (e)
  99   (let ((counter 0))
 100     (labels
 101         ((recurse (e) (cond ((and (not (atom e))
 102                                   (eq (caar e) '%pwilt_periodic))
 103                              (incf counter)
 104                              (mul -1 '$%k (cadr e)))
 105                             ((mplusp e)
 106                              (addn (mapcar #'recurse (cdr e)) t))
 107                             (t e))))
 108       (let ((result (recurse e)))
 109         (and (< counter 2)
 110              (freeof '%pwilt_periodic result)
 111              result)))))
 112
 113 ;;
 114 ;; Returns the degree of E in S, or NIL if E is not a polynomial in S
 115 ;;
 116 ;; XXX Note that this does not expand or simplify, so it may
 117 ;;     overestimate the degree if the highest-order term gets
 118 ;;     cancelled out!
 119 ;;
 120 (defun pwilt-polynomial-degree (e s)
 121   (labels
 122       ((recurse (e)
 123          (setq e (specrepcheck e))
 124          (cond ((alike1 s e) 1)
 125                ((freeof s e) 0)
 126                ((mplusp e) (reduce #'max (mapcar #'recurse (cdr e))))
 127                ((mtimesp e) (reduce #'+ (mapcar #'recurse (cdr e))))
 128                ((and (mexptp e)
 129                      (integerp (mexpt-expt e))
 130                      (>= (mexpt-expt e) 0))
 131                 (* (mexpt-expt e) (recurse (mexpt-base e))))
 132                (t (return-from pwilt-polynomial-degree nil)))))
 133     (recurse e)))
 134
 135 (defun pwilt1 (e s x)
 136   (setq e (specrepcheck e))
 137   (cond ((mbagp e) (cons (car e)
 138                          (mapcar #'(lambda (ee) (pwilt1 ee s x))
 139                                  (cdr e))))
 140         ((freeof s e) (mul e (take '($delta) x)))
 141         ((mplusp e) (addn (mapcar #'(lambda (ee) (pwilt1 ee s x))
 142                                   (cdr e))
 143                           t))
 144         ((and (null (atom e))
 145               (eq (caar e) '%laplace)
 146               (eq (cadddr e) s))
 147          (if (eq x (caddr e)) (cadr e)
 148            (subst x (caddr e) (cadr e))))
 149         (t (pwilt2 e s x))))
 150
 151 ;;
 152 ;; Calculate the inverse laplace transform by the residue method.
 153 ;;
 154 ;; This works only when e is of the form P(s)/Q(s), where P(s) and
 155 ;; Q(s) are polynomials in s, and where the degree of P(s) is less
 156 ;; than the degree of Q(s).
 157 ;;
 158 (defun pwilt2 (e s x)
 159   (flet ((fail () (return-from pwilt2 (list '(%ilt simp) e s x))))
 160     (let* ((n ($num e))
 161            (d ($denom e))
 162            (n-deg (pwilt-polynomial-degree n s))
 163            (d-deg (pwilt-polynomial-degree d s)))
 164       (unless (and n-deg d-deg (< n-deg d-deg))
 165         (let ((e2 ($partfrac e s)))
 166           (if (alike1 e e2)
 167               (fail)
 168             (return-from pwilt2 (pwilt1 e2 s x)))))
 169       (let* (($multiplicities nil)
 170              (solns (let (($breakup nil) ($programmode t)
 171                           ($globalsolve nil) ($solvedecomposes t)
 172                           ($solveexplicit t) ($solvefactors t))
 173                       ($solve d s)))
 174              ($demoivre t))
 175         (unless (and ($listp solns)
 176                      ($listp $multiplicities)
 177                      (= (length solns)
 178                         (length $multiplicities))
 179                      (= d-deg (reduce #'+ (cdr $multiplicities)
 180                                       :initial-value 0)))
 181           (fail))
 182         (addn (mapcar #'(lambda (soln m)
 183                           (unless (and (mequalp soln)
 184                                        (eq s (mequal-lhs soln)))
 185                             (fail))
 186                           (maxima-substitute
 187                            (mequal-rhs soln) s
 188                            ($diff (mul (power '$%e (mul s x))
 189                                        (div n ($diff d s m))
 190                                        m)
 191                                   s (- m 1))))
 192                       (cdr solns) (cdr $multiplicities))
 193               t)))))
 194
 195 (defun $pwilt (e s x)
 196   (setq e (specrepcheck e))
 197   (if (mbagp e)
 198       (cons (delsimp (car e))
 199             (mapcar #'(lambda (ee) ($pwilt ee s x))
 200                     (cdr e)))
 201     (let* ((xtemp (gensym))
 202            (periodic nil)
 203            (result
 204             (addn (mapcar #'(lambda (p)
 205                               (if (freeof '%pwilt_periodic (car p))
 206                                   (pwilt-one-part p s x xtemp)
 207                                 (let* ((offset
 208                                         (or (pwilt-subst-periodic (car p))
 209                                             (return-from $pwilt
 210                                               (list '(%pwilt) e s x))))
 211                                        (p (cons offset (cdr p))))
 212                                   (setq periodic t)
 213                                   (pwilt-one-part p s x xtemp))))
 214                           (pwilt-parts e s))
 215                   t)))
 216       (if periodic
 217           (list '(%sum simp) result '$%k 0 '$inf)
 218         result))))