share/diffequations/hstep.lisp

   1 ;;
   2 ;; Copyright (C) 2010, 2011 Mark H. Weaver <mhw@netris.org>
   3 ;;
   4 ;; hstep: Heaviside step function support for Maxima
   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 '$hstep 1 '$version)
  24
  25 (defprop $hstep %hstep verb)
  26 (defprop %hstep $hstep noun)
  27
  28 (defprop $hstep %hstep alias)
  29 (defprop %hstep $hstep reversealias)
  30
  31 (defprop %hstep simp-hstep operators)
  32 (setf (get '%hstep 'simplim%function) 'simplim%hstep)
  33
  34 (setf (get '%hstep 'real-valued) t)
  35
  36 ;; TODO: other properties which would be nice to declare about hstep:
  37 ;;   non-negative
  38 ;;   non-decreasing
  39
  40 (defprop %hstep ((x) (($delta) x)) grad)
  41 (defprop $delta ((x) ((%hstep) x)) integral)
  42
  43 (defun $hstep (z) (take '(%hstep) z))
  44
  45 ;;
  46 ;; TODO: should the following rule be included somehow?
  47 ;;
  48 ;;   hstep(-x)  -->  1 - hstep(x)
  49 ;;
  50 ;; It would also be nice to simplify products
  51 ;; containing more than one hstep.
  52 ;;
  53 (defun simp-hstep (expr z simpflag)
  54   (oneargcheck expr)
  55   (setq z (simpcheck (cadr expr) simpflag))
  56   (let ((sgn (csign z)))
  57     (cond ((eq sgn '$neg) 0)
  58           ((eq sgn '$zero) 1//2)
  59           ((eq sgn '$pos) 1)
  60           (t
  61            ;; positive * x --> x and negative * x --> -1 * x.
  62            (if (mtimesp z)
  63                (setq z (muln (mapcar #'(lambda (s)
  64                                          (let ((sgn (csign s)))
  65                                            (cond ((eq sgn '$neg) -1)
  66                                                  ((eq sgn '$pos) 1)
  67                                                  (t s))))
  68                                      (margs z))
  69                              t)))
  70            (eqtest (list '(%hstep) z) expr)))))
  71
  72 (defun simplim%hstep (e x pt)
  73   (let* ((e (limit (cadr e) x pt 'think))
  74          (sgn (mnqp e 0)))
  75     (cond ((eq t sgn) ($hstep e))     ;; limit of arg is not zero
  76           ((eq nil sgn) '$und)        ;; limit of arg is zero
  77           (t (throw 'limit nil)))))   ;; don't know