src/laplac.lisp

   1 ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
   2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   3 ;;;     The data in this file contains enhancments.                    ;;;;;
   4 ;;;                                                                    ;;;;;
   5 ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
   6 ;;;     All rights reserved                                            ;;;;;
   7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   8 ;;;     (c) Copyright 1981 Massachusetts Institute of Technology         ;;
   9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10
  11 (in-package :maxima)
  12
  13 (macsyma-module laplac)
  14
  15 (declare-top (special var $savefactors
  16                       checkfactors $ratfac $keepfloat *nounl* *nounsflag*
  17                       errcatch $errormsg))
  18
  19 ;;; The properties NOUN and VERB give correct linear display
  20
  21 (defprop $laplace %laplace verb)
  22 (defprop %laplace $laplace noun)
  23
  24 (defprop $ilt %ilt verb)
  25 (defprop %ilt $ilt noun)
  26
  27 (defun exponentiate (pow)
  28        ;;;COMPUTES %E**Z WHERE Z IS AN ARBITRARY EXPRESSION TAKING SOME OF THE WORK AWAY FROM SIMPEXPT
  29   (cond ((zerop1 pow) 1)
  30         ((equal pow 1) '$%e)
  31         (t (power '$%e pow))))
  32
  33 (defun fixuprest (rest)
  34        ;;;REST IS A PRODUCT WITHOUT THE MTIMES.FIXUPREST PUTS BACK THE MTIMES
  35   (cond ((null rest) 1)
  36         ((cdr rest) (cons '(mtimes) rest))
  37         (t (car rest))))
  38
  39 (defun isquadraticp (e x)
  40   (let ((b (sdiff e x)))
  41     (cond ((zerop1 b) (list 0 0 e))
  42           ((freeof x b) (list 0 b (maxima-substitute 0 x e)))
  43           ((setq b (islinear b x))
  44            (list (div* (car b) 2) (cdr b) (maxima-substitute 0 x e))))))
  45
  46 ;;;INITIALIZES SOME GLOBAL VARIABLES THEN CALLS THE DISPATCHING FUNCTION
  47
  48 (defmfun $laplace (fun var parm)
  49   (setq fun (mratcheck fun))
  50   (cond ((or *nounsflag* (member '%laplace *nounl* :test #'eq))
  51          (setq fun (remlaplace fun))))
  52    (laplace fun parm))
  53
  54 ;;;LAMBDA BINDS SOME SPECIAL VARIABLES TO NIL AND DISPATCHES
  55
  56 (defun remlaplace (e)
  57   (if (atom e)
  58       e
  59       (cons (delete 'laplace (append (car e) nil) :count 1 :test #'eq)
  60             (mapcar #'remlaplace (cdr e)))))
  61
  62 (defun laplace (fun parm &optional (dvar nil))
  63   (let ()
  64 ;;; Handles easy cases and calls appropriate function on others.
  65     (cond ((mbagp fun) (cons (car fun) (mapcar #'(lambda (e) (laplace e parm)) (cdr fun))))
  66           ((equal fun 0) 0)
  67           ((equal fun 1)
  68            (cond ((zerop1 parm) (simplify (list '($delta) 0)))
  69                  (t (power parm -1))))
  70           ((alike1 fun var) (power parm -2))
  71           ((or (atom fun) (freeof var fun))
  72            (cond ((zerop1 parm) (mul2 fun (simplify (list '($delta) 0))))
  73                  (t (mul2 fun (power parm -1)))))
  74           (t
  75            (let ((op (caar fun)))
  76              (let ((result ; We store the result of laplace for further work.
  77                (cond ((eq op 'mplus)
  78                       (laplus fun parm))
  79                      ((eq op 'mtimes)
  80                       (laptimes (cdr fun) parm))
  81                      ((eq op 'mexpt)
  82                       (lapexpt fun nil parm))
  83                      ((eq op '%sin)
  84                       (lapsin fun nil nil parm))
  85                      ((eq op '%cos)
  86                       (lapsin fun nil t parm))
  87                      ((eq op '%sinh)
  88                       (lapsinh fun nil nil parm))
  89                      ((eq op '%cosh)
  90                       (lapsinh fun nil t parm))
  91                      ((eq op '%log)
  92                       (laplog fun parm))
  93                      ((eq op '%derivative)
  94                       (lapdiff fun parm))
  95                      ((eq op '%integrate)
  96                       (lapint fun parm dvar))
  97                      ((eq op '%sum)
  98                       (list '(%sum)
  99                             (laplace (cadr fun) parm)
 100                             (caddr fun)
 101                             (cadddr fun)
 102                             (car (cddddr fun))))
 103                      ((eq op '%erf)
 104                       (laperf fun parm))
 105                      ((and (eq op '%ilt)(eq (cadddr fun) var))
 106                       (cond ((eq parm (caddr fun))(cadr fun))
 107                             (t (subst parm (caddr fun)(cadr fun)))))
 108                      ((eq op '$delta)
 109                       (lapdelta fun nil parm))
 110                      ((member op '(%hstep $unit_step))
 111                       (laphstep fun nil parm))
 112                      ((setq op ($get op '$laplace))
 113                       (mcall op fun var parm))
 114                      (t (lapdefint fun parm)))))
 115               (when (isinop result '%integrate)
 116                 ;; Laplace has not found a result but returns a definit
 117                 ;; integral. This integral can contain internal integration
 118                 ;; variables. Replace such a result with the noun form.
 119                 (setq result (list '(%laplace) fun var parm)))
 120               ;; Check if we have a result, when not call $specint.
 121               (check-call-to-$specint result fun parm)))))))
 122
 123 ;;; Check if laplace has found a result, when not try $specint.
 124
 125 (defun check-call-to-$specint (result fun parm)
 126   (cond
 127     ((or (isinop result '%laplace)
 128          (isinop result '%limit)   ; Try $specint for incomplete results
 129          (isinop result '%at))     ; which contain %limit or %at too.
 130      ;; laplace returns a noun form or a result which contains %limit or %at.
 131      ;; We pass the function to $specint to look for more results.
 132      (let (res)
 133        ;; laplace assumes the parameter s to be positive and does a
 134        ;; declaration before an integration is done. Therefore we declare
 135        ;; the parameter of the Laplace transform to be positive before
 136        ;; we call $specint too.
 137        (with-new-context (context)
 138          (meval `(($assume) ,@(list (list '(mgreaterp) parm 0))))
 139          (setq res ($specint (mul fun (power '$%e (mul -1 var parm))) var)))
 140        (if (or (isinop res '%specint)  ; Both symobls are possible, that is
 141                (isinop res '$specint)) ; not consistent! Check it! 02/2009
 142            ;; $specint has not found a result.
 143            result
 144            ;; $specint has found a result
 145            res)))
 146        (t result)))
 147
 148 (defun laplus (fun parm)
 149   (simplus (cons '(mplus) (mapcar #'(lambda (e) (laplace e parm)) (cdr fun))) 1 t))
 150
 151 (defun laptimes (fun parm)
 152        ;;;EXPECTS A LIST (PERHAPS EMPTY) OF FUNCTIONS MULTIPLIED TOGETHER WITHOUT THE MTIMES
 153        ;;;SEES IF IT CAN APPLY THE FIRST AS A TRANSFORMATION ON THE REST OF THE FUNCTIONS
 154   (cond ((null fun) (list '(mexpt) parm -1.))
 155         ((null (cdr fun)) (laplace (car fun) parm))
 156         ((freeof var (car fun))
 157          (simptimes (list '(mtimes)
 158                           (car fun)
 159                           (laptimes (cdr fun) parm))
 160                     1
 161                     t))
 162         ((eq (car fun) var)
 163          (simptimes (list '(mtimes) -1 (sdiff (laptimes (cdr fun) parm) parm))
 164                     1
 165                     t))
 166         (t
 167          (let ((op (caaar fun)))
 168            (cond ((eq op 'mexpt)
 169                   (lapexpt (car fun) (cdr fun) parm))
 170                  ((eq op 'mplus)
 171                   (laplus ($multthru (fixuprest (cdr fun)) (car fun)) parm))
 172                  ((eq op '%sin)
 173                   (lapsin (car fun) (cdr fun) nil parm))
 174                  ((eq op '%cos)
 175                   (lapsin (car fun) (cdr fun) t parm))
 176                  ((eq op '%sinh)
 177                   (lapsinh (car fun) (cdr fun) nil parm))
 178                  ((eq op '%cosh)
 179                   (lapsinh (car fun) (cdr fun) t parm))
 180                  ((eq op '$delta)
 181                   (lapdelta (car fun) (cdr fun) parm))
 182                  ((member op '(%hstep $unit_step))
 183                   (laphstep (car fun) (cdr fun) parm))
 184                  (t
 185                   (lapshift (car fun) (cdr fun) parm)))))))
 186
 187 (defun lapexpt (fun rest parm)
 188        ;;;HANDLES %E**(A*T+B)*REST(T), %E**(A*T**2+B*T+C),
 189        ;;; 1/SQRT(A*T+B), OR T**K*REST(T)
 190   (prog (ab base-of-fun power result)
 191      (setq base-of-fun (cadr fun) power (caddr fun))
 192      (cond
 193        ((and
 194          (freeof var base-of-fun)
 195          (setq
 196           ab
 197           (isquadraticp
 198            (cond ((eq base-of-fun '$%e) power)
 199                  (t (simptimes (list '(mtimes)
 200                                      power
 201                                      (list '(%log)
 202                                            base-of-fun))
 203                                1
 204                                nil)))
 205            var)))
 206         (cond ((equal (car ab) 0) (go %e-case-lin))
 207               ((null rest) (go %e-case-quad))
 208               (t (go noluck))))
 209        ((and (eq base-of-fun var) (freeof var power))
 210         (go var-case))
 211        ((and (alike1 '((rat) -1 2) power) (null rest)
 212              (setq ab (islinear base-of-fun var)))
 213         (setq result (div* (cdr ab) (car ab)))
 214         (return (simptimes
 215                  (list '(mtimes)
 216                        (list '(mexpt)
 217                              (div* '$%pi
 218                                    (list '(mtimes)
 219                                          (car ab)
 220                                          parm))
 221                              '((rat) 1 2))
 222                        (exponentiate (list '(mtimes) result parm))
 223                        (list '(mplus)
 224                              1
 225                              (list '(mtimes)
 226                                    -1
 227                                    (list '(%erf)
 228                                          (list '(mexpt)
 229                                                (list '(mtimes)
 230                                                      result
 231                                                      parm)
 232                                                '((rat) 1 2)))
 233                                    ))) 1 nil)))
 234        (t (go noluck)))
 235      %e-case-lin
 236      (setq result
 237       (cond
 238         (rest (sratsimp ($at (laptimes rest parm)
 239                              (list '(mequal)
 240                                    parm
 241                                    (list '(mplus)
 242                                          parm
 243                                          (afixsign (cadr ab)
 244                                                    nil))))))
 245         (t (list '(mexpt)
 246                  (list '(mplus)
 247                        parm
 248                        (afixsign (cadr ab) nil))
 249                  -1))))
 250      (return (simptimes (list '(mtimes)
 251                               (exponentiate (caddr ab))
 252                               result)
 253                         1
 254                         nil))
 255      %e-case-quad
 256      (setq result (afixsign (car ab) nil))
 257      (setq
 258       result
 259       (list
 260        '(mtimes)
 261        (div* (list '(mexpt)
 262                    (div* '$%pi result)
 263                    '((rat) 1 2))
 264              2)
 265        (exponentiate (div* (list '(mexpt) parm 2)
 266                            (list '(mtimes) 4 result)))
 267        (list '(mplus)
 268              1
 269              (list '(mtimes)
 270                    -1
 271                    (list '(%erf)
 272                          (div* parm
 273                                (list '(mtimes)
 274                                      2
 275                                      (list '(mexpt)
 276                                            result
 277                                            '((rat) 1 2)))))
 278                    ))))
 279      (and (null (equal (cadr ab) 0))
 280           (setq result
 281                 (maxima-substitute (list '(mplus)
 282                                          parm
 283                                          (list '(mtimes)
 284                                                -1
 285                                                (cadr ab)))
 286                                    parm
 287                                    result)))
 288      (return (simptimes  (list '(mtimes)
 289                                (exponentiate (caddr ab))
 290                                result) 1 nil))
 291      var-case
 292      (cond ((or (null rest) (freeof var (fixuprest rest)))
 293             (go var-easy-case)))
 294      (cond ((posint power)
 295             (return (afixsign (apply '$diff
 296                                      (list (laptimes rest parm)
 297                                            parm
 298                                            power))
 299                               (even power))))
 300            ((negint power)
 301             (return (mydefint (hackit power rest parm)
 302                               (createname parm (- power))
 303                               parm parm)))
 304            (t (go noluck)))
 305      var-easy-case
 306      (setq power
 307            (simplus (list '(mplus) 1 power) 1 t))
 308      (or (eq (asksign ($realpart power)) '$positive) (go noluck))
 309      (setq result (list (list '(%gamma) power)
 310                         (list '(mexpt)
 311                               parm
 312                               (afixsign power nil))))
 313      (and rest (setq result (nconc result rest)))
 314      (return (simptimes (cons '(mtimes) result) 1 nil))
 315      noluck
 316      (return
 317        (cond
 318          ((and (posint power)
 319                (member (caar base-of-fun)
 320                      '(mplus %sin %cos %sinh %cosh) :test #'eq))
 321           (laptimes (cons base-of-fun
 322                           (cons (cond ((= power 2) base-of-fun)
 323                                       (t (list '(mexpt)
 324                                                base-of-fun
 325                                                (1- power))))
 326                                 rest)) parm))
 327          (t (lapshift fun rest parm))))))
 328
 329 ;;;INTEGRAL FROM A TO INFINITY OF F(X)
 330 (defun mydefint (f x a parm)
 331   (let ((tryint (and (not ($unknown f))
 332                      ;; We don't want any errors thrown by $defint to propagate,
 333                      ;; so we set errset, errcatch, $errormsg and *mdebug*.  If
 334                      ;; $defint throws a error, then no error message is printed
 335                      ;; and errset simply returns nil below.
 336                      (with-new-context (context)
 337                        (meval `(($assume) ,@(list (list '(mgreaterp) parm 0))))
 338                        (meval `(($assume) ,@(list (list '(mgreaterp) x 0))))
 339                        (meval `(($assume) ,@(list (list '(mgreaterp) a 0))))
 340                        (let ((errset nil)
 341                              (errcatch t)
 342                              ($errormsg nil)
 343                              (*mdebug* nil))
 344                          (errset ($defint f x a '$inf)))))))
 345     (if tryint
 346         (car tryint)
 347         (list '(%integrate) f x a '$inf))))
 348
 349  ;;;CREATES UNIQUE NAMES FOR VARIABLE OF INTEGRATION
 350 (defun createname (head tail)
 351   (intern (format nil "~S~S" head tail)))
 352
 353 ;;;REDUCES LAPLACE(F(T)/T**N,T,S) CASE TO LAPLACE(F(T)/T**(N-1),T,S) CASE
 354 (defun hackit (exponent rest parm)
 355   (cond ((equal exponent -1)
 356          (let ((parm (createname parm 1)))
 357            (laptimes rest parm)))
 358         (t (mydefint (hackit (1+ exponent) rest parm)
 359                      (createname parm (- -1 exponent))
 360                      (createname parm (- exponent)) parm))))
 361
 362 (defun afixsign (funct signswitch)
 363        ;;;MULTIPLIES FUNCT BY -1 IF SIGNSWITCH IS NIL
 364   (cond (signswitch funct)
 365         (t (simptimes (list '(mtimes) -1 funct) 1 t))))
 366
 367 (defun lapshift (fun rest parm)
 368   (cond ((atom fun) (merror "LAPSHIFT: expected a cons, not ~M" fun))
 369         ((or (member 'laplace (car fun) :test #'eq) (null rest))
 370          (lapdefint (cond (rest (simptimes (cons '(mtimes)
 371                                                  (cons fun rest)) 1 t))
 372                           (t fun)) parm))
 373         (t (laptimes (append rest
 374                              (ncons (cons (append (car fun)
 375                                                   '(laplace))
 376                                           (cdr fun)))) parm))))
 377
 378 ;;;COMPUTES %E**(W*B*%I)*F(S-W*A*%I) WHERE W=-1 IF SIGN IS T ELSE W=1
 379 (defun mostpart (f parm sign a b)
 380   (let ((substinfun ($at f
 381                          (list '(mequal)
 382                                parm
 383                                (list '(mplus) parm (afixsign (list '(mtimes) a '$%i) sign))))))
 384     (if (zerop1 b)
 385         substinfun
 386         (list '(mtimes)
 387               (exponentiate (afixsign (list '(mtimes) b '$%i) (null sign)))
 388               substinfun))))
 389
 390  ;;;IF WHICHSIGN IS NIL THEN SIN TRANSFORM ELSE COS TRANSFORM
 391 (defun compose (fun parm whichsign a b)
 392   (let ((result (list '((rat) 1 2)
 393                       (list '(mplus)
 394                             (mostpart fun parm t a b)
 395                             (afixsign (mostpart fun parm nil a b)
 396                                       whichsign)))))
 397     (sratsimp (simptimes (cons '(mtimes)
 398                                (if whichsign
 399                                    result
 400                                    (cons '$%i result)))
 401                          1 nil))))
 402
 403  ;;;FUN IS OF THE FORM SIN(A*T+B)*REST(T) OR COS
 404 (defun lapsin (fun rest trigswitch parm)
 405   (let ((ab (islinear (cadr fun) var)))
 406     (cond (ab
 407            (cond (rest
 408                   (compose (laptimes rest parm)
 409                            parm
 410                            trigswitch
 411                            (car ab)
 412                            (cdr ab)))
 413                  (t
 414                   (simptimes
 415                    (list '(mtimes)
 416                     (cond ((zerop1 (cdr ab))
 417                            (if trigswitch parm (car ab)))
 418                           (t
 419                            (cond (trigswitch
 420                                   (list '(mplus)
 421                                         (list '(mtimes)
 422                                               parm
 423                                               (list '(%cos) (cdr ab)))
 424                                         (list '(mtimes)
 425                                               -1
 426                                               (car ab)
 427                                               (list '(%sin) (cdr ab)))))
 428                                  (t
 429                                   (list '(mplus)
 430                                         (list '(mtimes)
 431                                               parm
 432                                               (list '(%sin) (cdr ab)))
 433                                         (list '(mtimes)
 434                                               (car ab)
 435                                               (list '(%cos) (cdr ab))))))))
 436                     (list '(mexpt)
 437                           (list '(mplus)
 438                                 (list '(mexpt) parm 2)
 439                                 (list '(mexpt) (car ab) 2))
 440                           -1))
 441                    1 nil))))
 442           (t
 443            (lapshift fun rest parm)))))
 444
 445  ;;;FUN IS OF THE FORM SINH(A*T+B)*REST(T) OR IS COSH
 446 (defun lapsinh (fun rest switch parm)
 447   (cond ((islinear (cadr fun) var)
 448          (sratsimp
 449           (laplus
 450            (simplus
 451             (list '(mplus)
 452                   (nconc (list '(mtimes)
 453                                (list '(mexpt)
 454                                      '$%e
 455                                      (cadr fun))
 456                                '((rat) 1 2))
 457                          rest)
 458                   (afixsign (nconc (list '(mtimes)
 459                                          (list '(mexpt)
 460                                                '$%e
 461                                                (afixsign (cadr fun)
 462                                                          nil))
 463                                          '((rat) 1 2))
 464                                    rest)
 465                             switch))
 466             1
 467             nil) parm)))
 468         (t (lapshift fun rest parm))))
 469
 470  ;;;FUN IS OF THE FORM LOG(A*T)
 471 (defun laplog (fun parm)
 472   (let ((ab (islinear (cadr fun) var)))
 473     (cond ((and ab (zerop1 (cdr ab)))
 474            (simptimes (list '(mtimes)
 475                             (list '(mplus)
 476                                   (subfunmake '$psi '(0) (ncons 1))
 477                                   (list '(%log) (car ab))
 478                                   (list '(mtimes) -1 (list '(%log) parm)))
 479                             (list '(mexpt) parm -1))
 480                       1 nil))
 481           (t
 482            (lapdefint fun parm)))))
 483
 484 (defun raiseup (fbase exponent)
 485   (if (equal exponent 1)
 486       fbase
 487       (list '(mexpt) fbase exponent)))
 488
 489 ;;TAKES TRANSFORM OF HSTEP(A*T+B)*F(T)
 490 (defun laphstep (fun rest parm)
 491   (let ((ab (islinear (cadr fun) var)))
 492     (if ab
 493         (let* ((a (car ab))
 494                (b (cdr ab))
 495                (offset (div b a)))
 496           (case (asksign a)
 497             ($positive (if (eq (asksign offset) '$negative)
 498                            (mul (exponentiate (mul offset parm))
 499                                 (laplace (maxima-substitute
 500                                           (sub var offset) var
 501                                           (fixuprest rest)) parm))
 502                          (laptimes rest parm)))
 503             ($negative (if (eq (asksign offset) '$negative)
 504                            (sub (laptimes rest parm)
 505                                 (mul (exponentiate (mul offset parm))
 506                                      (laplace (maxima-substitute
 507                                                (sub var offset) var
 508                                                (fixuprest rest)) parm)))
 509                          0))
 510             (t (mul fun (laptimes rest parm)))))
 511       (lapshift fun rest parm))))
 512
 513 ;;TAKES TRANSFORM OF DELTA(A*T+B)*F(T)
 514 (defun lapdelta (fun rest parm)
 515   (let ((ab (islinear (cadr fun) var))
 516         (sign nil)
 517         (recipa nil))
 518     (cond (ab
 519            (setq recipa (power (car ab) -1) ab (div (cdr ab) (car ab)))
 520            (setq sign (asksign ab) recipa (simplifya (list '(mabs) recipa) nil))
 521            (simplifya (cond ((eq sign '$positive)
 522                              0)
 523                             ((eq sign '$zero)
 524                              (list '(mtimes)
 525                                    (maxima-substitute 0 var (fixuprest rest))
 526                                    recipa))
 527                             (t
 528                              (list '(mtimes)
 529                                    (maxima-substitute (neg ab) var (fixuprest rest))
 530                                    (list '(mexpt) '$%e (cons '(mtimes) (cons parm (ncons ab))))
 531                                    recipa)))
 532                       nil))
 533           (t
 534            (lapshift fun rest parm)))))
 535
 536 (defun laperf (fun parm)
 537   (let ((ab (islinear (cadr fun) var)))
 538     (cond ((and ab (equal (cdr ab) 0))
 539            (simptimes (list '(mtimes)
 540                             (div* (exponentiate (div* (list '(mexpt) parm 2)
 541                                                       (list '(mtimes)
 542                                                             4
 543                                                             (list '(mexpt) (car ab) 2))))
 544                                   parm)
 545                             (list '(mplus)
 546                                   1
 547                                   (list '(mtimes)
 548                                         -1
 549                                         (list '(%erf) (div* parm (list '(mtimes) 2 (car ab)))))))
 550                       1
 551                       nil))
 552           (t
 553            (lapdefint fun parm)))))
 554
 555 (defun lapdefint (fun parm)
 556   (prog (tryint mult)
 557      (and ($unknown fun)(go skip))
 558      (setq mult (simptimes (list '(mtimes) (exponentiate
 559                                             (list '(mtimes) -1 var parm)) fun) 1 nil))
 560      (with-new-context (context)
 561        (meval `(($assume) ,@(list (list '(mgreaterp) parm 0))))
 562        (setq tryint
 563              ;; We don't want any errors thrown by $defint to propagate,
 564              ;; so we set errset, errcatch, $errormsg and *mdebug*.  If
 565              ;; $defint throws a error, then no error message is printed
 566              ;; and errset simply returns nil below.
 567              (let ((errset nil)
 568                    (errcatch t)
 569                    ($errormsg nil)
 570                    (*mdebug* nil))
 571                (errset ($defint mult var 0 '$inf)))))
 572      (and tryint (not (eq (and (listp (car tryint))
 573                                (caaar tryint))
 574                           '%integrate))
 575           (return (car tryint)))
 576      skip (return (list '(%laplace) fun var parm))))
 577
 578
 579 (defun lapdiff (fun parm)
 580 ;;;FUN IS OF THE FORM DIFF(F(T),T,N) WHERE N IS A POSITIVE INTEGER
 581   (prog (difflist degree frontend resultlist newdlist order arg2)
 582      (setq newdlist (setq difflist (copy-tree (cddr fun))))
 583      (setq arg2 (list '(mequal) var 0))
 584      a    (cond ((null difflist)
 585                  (return (cons '(%derivative)
 586                                (cons (list '(%laplace)
 587                                            (cadr fun)
 588                                            var
 589                                            parm)
 590                                      newdlist))))
 591                 ((eq (car difflist) var)
 592                  (setq degree (cadr difflist)
 593                        difflist (cddr difflist))
 594                  (go out)))
 595      (setq difflist (cdr (setq frontend (cdr difflist))))
 596      (go a)
 597      out  (cond ((null (posint degree))
 598                  (return (list '(%laplace) fun var parm))))
 599      (cond (frontend (rplacd frontend difflist))
 600            (t (setq newdlist difflist)))
 601      (cond (newdlist (setq fun (cons '(%derivative)
 602                                      (cons (cadr fun)
 603                                            newdlist))))
 604            (t (setq fun (cadr fun))))
 605      (setq order 0)
 606      loop (decf degree)
 607      (setq resultlist
 608            (cons (list '(mtimes)
 609                        (raiseup parm degree)
 610                        ($at ($diff fun var order) arg2))
 611                  resultlist))
 612      (incf order)
 613      (and (> degree 0) (go loop))
 614      (setq resultlist (cond ((cdr resultlist)
 615                              (cons '(mplus)
 616                                    resultlist))
 617                             (t (car resultlist))))
 618      (return (simplus (list '(mplus)
 619                             (list '(mtimes)
 620                                   (raiseup parm order)
 621                                   (laplace fun parm))
 622                             (list '(mtimes)
 623                                   -1
 624                                   resultlist))
 625                       1 nil))))
 626
 627  ;;;FUN IS OF THE FORM INTEGRATE(F(X)*G(T)*H(T-X),X,0,T)
 628 (defun lapint (fun parm dvar)
 629   (prog (newfun parm-list f var-list var-parm-list)
 630      (and dvar (go convolution))
 631      (setq dvar (cadr (setq newfun (cdr fun))))
 632      (and (cddr newfun)
 633           (zerop1 (caddr newfun))
 634           (eq (cadddr newfun) var)
 635           (go convolutiontest))
 636      notcon
 637      (setq newfun (cdr fun))
 638      (cond ((cddr newfun)
 639             (cond ((and (freeof var (caddr newfun))
 640                         (freeof var (cadddr newfun)))
 641                    (return (list '(%integrate)
 642                                  (laplace (car newfun) parm dvar)
 643                                  dvar
 644                                  (caddr newfun)
 645                                  (cadddr newfun))))
 646                   (t (go giveup))))
 647            (t (return (list '(%integrate)
 648                             (laplace (car newfun) parm dvar)
 649                             dvar))))
 650      giveup
 651      (return (list '(%laplace) fun var parm))
 652      convolutiontest
 653      (setq newfun ($factor (car newfun)))
 654      (cond ((eq (caar newfun) 'mtimes)
 655             (setq f (cadr newfun) newfun (cddr newfun)))
 656            (t (setq f newfun newfun nil)))
 657      gothrulist
 658      (cond ((freeof dvar f)
 659             (setq parm-list (cons f parm-list)))
 660            ((freeof var f) (setq var-list (cons f var-list)))
 661            ((freeof dvar
 662                     (sratsimp (maxima-substitute (list '(mplus)
 663                                                        var
 664                                                        dvar)
 665                                                  var
 666                                                  f)))
 667             (setq var-parm-list (cons f var-parm-list)))
 668            (t (go notcon)))
 669      (cond (newfun (setq f (car newfun) newfun (cdr newfun))
 670                    (go gothrulist)))
 671      (and
 672       parm-list
 673       (return
 674         (laplace
 675          (cons
 676           '(mtimes)
 677           (nconc parm-list
 678                  (ncons (list '(%integrate)
 679                               (cons '(mtimes)
 680                                     (append var-list
 681                                             var-parm-list))
 682                               dvar
 683                               0
 684                               var)))) parm dvar)))
 685      convolution
 686      (return
 687        (simptimes
 688         (list
 689          '(mtimes)
 690          (laplace ($expand (maxima-substitute var
 691                                               dvar
 692                                               (fixuprest var-list))) parm dvar)
 693          (laplace
 694           ($expand (maxima-substitute 0 dvar (fixuprest var-parm-list))) parm dvar))
 695         1
 696         t))))
 697
 698 (declare-top (special varlist ratform ils ilt))
 699
 700 (defmfun $ilt (exp ils ilt)
 701  ;;;EXP IS F(S)/G(S) WHERE F AND G ARE POLYNOMIALS IN S AND DEGR(F) < DEGR(G)
 702   (let (varlist ($savefactors t) checkfactors $ratfac $keepfloat)
 703                 ;;; MAKES ILS THE MAIN VARIABLE
 704     (setq varlist (list ils))
 705     (newvar exp)
 706     (orderpointer varlist)
 707     (setq var (caadr (ratrep* ils)))
 708     (cond ((and (null (atom exp))
 709                 (eq (caar exp) 'mequal))
 710            (list '(mequal)
 711                  ($ilt (cadr exp) ils ilt)
 712                  ($ilt (caddr exp) ils ilt)))
 713           ((zerop1 exp) 0)
 714           ((freeof ils exp)
 715            (list '(mtimes) exp (list '($delta) ilt)))
 716           (t (ilt0 exp)))))
 717
 718 (defun maxima-rationalp (le v)
 719   (cond ((null le))
 720         ((and (null (atom (car le))) (null (freeof v (car le))))
 721          nil)
 722         (t (maxima-rationalp (cdr le) v))))
 723
 724  ;;;THIS FUNCTION DOES THE PARTIAL FRACTION DECOMPOSITION
 725 (defun ilt0 (exp)
 726   (prog (wholepart frpart num denom y content real factor
 727          apart bpart parnumer ratarg ratform)
 728      (and (mplusp exp)
 729           (return (simplus  (cons '(mplus)
 730                                   (mapcar #'(lambda (f) ($ilt f ils ilt)) (cdr exp))) 1 t)))
 731      (and (null (atom exp))
 732           (eq (caar exp) '%laplace)
 733           (eq (cadddr exp) ils)
 734           (return (cond ((eq (caddr exp) ilt) (cadr exp))
 735                         (t (subst ilt
 736                                   (caddr exp)
 737                                   (cadr exp))))))
 738      (setq ratarg (ratrep* exp))
 739      (or (maxima-rationalp varlist ils)
 740          (return (list '(%ilt) exp ils ilt)))
 741      (setq ratform (car ratarg))
 742      (setq denom (ratdenominator (cdr ratarg)))
 743      (setq frpart (pdivide (ratnumerator (cdr ratarg)) denom))
 744      (setq wholepart (car frpart))
 745      (setq frpart (ratqu (cadr frpart) denom))
 746      (cond ((not (zerop1 (car wholepart)))
 747             (return (list '(%ilt) exp ils ilt)))
 748            ((zerop1 (car frpart)) (return 0)))
 749      (setq num (car frpart) denom (cdr frpart))
 750      (setq y (oldcontent denom))
 751      (setq content (car y))
 752      (setq real (cadr y))
 753      (setq factor (pfactor real))
 754      loop (cond ((null (cddr factor))
 755                  (setq apart real
 756                        bpart 1
 757                        y '((0 . 1) 1 . 1))
 758                  (go skip)))
 759      (setq apart (pexpt (car factor) (cadr factor)))
 760      (setq bpart (car (ratqu real apart)))
 761      (setq y (bprog apart bpart))
 762      skip (setq frpart
 763                 (cdr (ratdivide (ratti (ratnumerator num)
 764                                        (cdr y)
 765                                        t)
 766                                 (ratti (ratdenominator num)
 767                                        (ratti content apart t)
 768                                        t))))
 769      (setq
 770       parnumer
 771       (cons (ilt1 (ratqu (ratnumerator frpart)
 772                          (ratti (ratdenominator frpart)
 773                                 (ratti (ratdenominator num)
 774                                        content
 775                                        t)
 776                                 t))
 777                   (car factor)
 778                   (cadr factor))
 779             parnumer))
 780      (setq factor (cddr factor))
 781      (cond ((null factor)
 782             (return (simplus (cons '(mplus) parnumer) 1 t))))
 783      (setq num (cdr (ratdivide (ratti num (car y) t)
 784                                (ratti content bpart t))))
 785      (setq real bpart)
 786      (go loop)))
 787
 788 (declare-top (special q z))
 789
 790 (defun ilt1 (p q k)
 791   (let (z)
 792     (cond ((onep1 k) (ilt3 p ))
 793           (t (setq z (bprog q (pderivative q var)))
 794              (ilt2 p k)))))
 795
 796
 797  ;;;INVERTS P(S)/Q(S)**K WHERE Q(S)  IS IRREDUCIBLE
 798  ;;;DOESN'T CALL ILT3 IF Q(S) IS LINEAR
 799 (defun ilt2 (p k)
 800   (prog (y a b)
 801      (and (onep1 k) (return (ilt3 p)))
 802      (decf k)
 803      (setq a (ratti p (car z) t))
 804      (setq b (ratti p (cdr z) t))
 805      (setq y (pexpt q k))
 806      (cond
 807        ((or (null (equal (pdegree q var) 1))
 808             (> (pdegree (car p) var) 0))
 809         (return
 810           (simplus
 811            (list
 812             '(mplus)
 813             (ilt2
 814              (cdr (ratdivide (ratplus a (ratqu (ratderivative b var) k)) y))
 815              k)
 816             ($multthru (simptimes (list '(mtimes)
 817                                         ilt
 818                                         (power k -1)
 819                                         (ilt2 (cdr (ratdivide b y)) k))
 820                                   1
 821                                   t)))
 822            1
 823            t))))
 824      (setq a (disrep (polcoef q 1))
 825            b (disrep (polcoef q 0)))
 826      (return
 827        (simptimes (list '(mtimes)
 828                         (disrep p)
 829                         (raiseup ilt k)
 830                         (simpexpt (list '(mexpt)
 831                                         '$%e
 832                                         (list '(mtimes)
 833                                               -1
 834                                               ilt
 835                                               b
 836                                               (list '(mexpt) a -1)))
 837                                   1
 838                                   nil)
 839                         (list '(mexpt)
 840                               a
 841                               (- -1 k))
 842                         (list '(mexpt)
 843                               (factorial k)
 844                               -1))
 845                   1
 846                   nil))))
 847
 848 ;;(DEFUN COEF MACRO (POL) (SUBST (CADR POL) (QUOTE DEG)
 849 ;;  '(DISREP (RATQU (POLCOEF (CAR P) DEG) (CDR P)))))
 850
 851 (defmacro coef (pol)
 852   `(disrep (ratqu (polcoef (car p) ,pol) (cdr p))))
 853
 854 (defun lapsum (&rest args)
 855   (cons '(mplus) args))
 856
 857 (defun lapprod (&rest args)
 858   (cons '(mtimes) args))
 859
 860 (defun expo (&rest args)
 861   (cons '(mexpt) args))
 862
 863 ;;;INVERTS P(S)/Q(S) WHERE Q(S) IS IRREDUCIBLE
 864 (defun ilt3 (p)
 865   (prog (discrim sign a c d e b1 b0 r term1 term2 degr)
 866      (setq e (disrep (polcoef q 0))
 867            d (disrep (polcoef q 1))
 868            degr (pdegree q var))
 869      (and (equal degr 1)
 870           (return
 871             (simptimes (lapprod
 872                         (disrep p)
 873                         (expo d -1)
 874                         (expo '$%e (lapprod -1 ilt e (expo d -1))))
 875                        1
 876                        nil)))
 877      (setq c (disrep (polcoef q 2)))
 878      (and (equal degr 2) (go quadratic))
 879      (and (equal degr 3) (zerop1 c) (zerop1 d)
 880           (go cubic))
 881      (return (list '(%ilt) (div* (disrep p)(disrep q)) ils ilt))
 882      cubic (setq  a (disrep (polcoef q 3))
 883                   r (simpnrt (div* e a) 3))
 884      (setq d (div* (disrep p)(lapprod a (lapsum
 885                                          (expo ils 3)(expo '%r 3)))))
 886      (return (ilt0 (maxima-substitute r '%r ($partfrac d ils))))
 887      quadratic (setq b0 (coef 0) b1 (coef 1))
 888
 889      (setq discrim
 890            (simplus (lapsum
 891                      (lapprod 4 e c)
 892                      (lapprod -1 d d))
 893                     1
 894                     nil))
 895      (setq sign (cond ((free discrim '$%i) (asksign discrim)) (t '$positive))
 896            term1 '(%cos)
 897            term2 '(%sin))
 898      (setq degr (expo '$%e (lapprod ilt d (power c -1) '((rat) -1 2))))
 899      (cond ((eq sign '$zero)
 900             (return (simptimes (lapprod degr (lapsum (div* b1 c)
 901                                                      (lapprod
 902                                                       (div* (lapsum (lapprod 2 b0 c) (lapprod -1 b1 d))
 903                                                             (lapprod 2 c c)) ilt))) 1 nil))
 904             )              ((eq sign '$negative)
 905             (setq term1 '(%cosh)
 906                   term2 '(%sinh)
 907                   discrim (simptimes (lapprod -1 discrim) 1 t))))
 908      (setq discrim (simpnrt discrim 2))
 909      (setq sign
 910       (simptimes
 911        (lapprod
 912         (lapsum
 913          (lapprod 2 b0 c)
 914          (lapprod -1 b1 d))
 915         (expo discrim -1))
 916        1
 917        nil))
 918      (setq c (power c -1))
 919      (setq discrim (simptimes (lapprod
 920                                discrim
 921                                ilt
 922                                '((rat) 1 2)
 923                                c)
 924                               1
 925                               t))
 926      (return
 927        (simptimes
 928         (lapprod c degr
 929          (lapsum
 930           (lapprod b1 (list term1 discrim))
 931           (lapprod sign (list term2 discrim))))
 932         1
 933         nil))))
 934
 935 (declare-top (unspecial ils ilt *nounl* q ratform var
 936                         varlist z))