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