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 (progn 140 (meval `(($assume) ,@(list (list '(mgreaterp) parm 0)))) 141 (setq res ($specint (mul fun (power '$%e (mul -1 var parm))) var)))) 142 (if (or (isinop res '%specint) ; Both symobls are possible, that is 143 (isinop res '$specint)) ; not consistent! Check it! 02/2009 144 ;; $specint has not found a result. 145 result 146 ;; $specint has found a result 147 res))) 148 (t result))) 149 150(defun laplus (fun parm) 151 (simplus (cons '(mplus) (mapcar #'(lambda (e) (laplace e parm)) (cdr fun))) 1 t)) 152 153(defun laptimes (fun parm) 154 ;;;EXPECTS A LIST (PERHAPS EMPTY) OF FUNCTIONS MULTIPLIED TOGETHER WITHOUT THE MTIMES 155 ;;;SEES IF IT CAN APPLY THE FIRST AS A TRANSFORMATION ON THE REST OF THE FUNCTIONS 156 (cond ((null fun) (list '(mexpt) parm -1.)) 157 ((null (cdr fun)) (laplace (car fun) parm)) 158 ((freeof var (car fun)) 159 (simptimes (list '(mtimes) 160 (car fun) 161 (laptimes (cdr fun) parm)) 162 1 163 t)) 164 ((eq (car fun) var) 165 (simptimes (list '(mtimes) -1 (sdiff (laptimes (cdr fun) parm) parm)) 166 1 167 t)) 168 (t 169 (let ((op (caaar fun))) 170 (cond ((eq op 'mexpt) 171 (lapexpt (car fun) (cdr fun) parm)) 172 ((eq op 'mplus) 173 (laplus ($multthru (fixuprest (cdr fun)) (car fun)) parm)) 174 ((eq op '%sin) 175 (lapsin (car fun) (cdr fun) nil parm)) 176 ((eq op '%cos) 177 (lapsin (car fun) (cdr fun) t parm)) 178 ((eq op '%sinh) 179 (lapsinh (car fun) (cdr fun) nil parm)) 180 ((eq op '%cosh) 181 (lapsinh (car fun) (cdr fun) t parm)) 182 ((eq op '$delta) 183 (lapdelta (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 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 ;; $defint should not throw a Maxima error, 333 ;; therefore we set the flags errcatch and $errormsg. 334 ;; errset catches the error and returns nil 335 (with-new-context (context) 336 (progn 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 ((errcatch t) ($errormsg nil)) 341 (errset ($defint f x a '$inf)))))))) 342 (if tryint 343 (car tryint) 344 (list '(%integrate) f x a '$inf)))) 345 346 ;;;CREATES UNIQUE NAMES FOR VARIABLE OF INTEGRATION 347(defun createname (head tail) 348 (intern (format nil "~S~S" head tail))) 349 350;;;REDUCES LAPLACE(F(T)/T**N,T,S) CASE TO LAPLACE(F(T)/T**(N-1),T,S) CASE 351(defun hackit (exponent rest parm) 352 (cond ((equal exponent -1) 353 (let ((parm (createname parm 1))) 354 (laptimes rest parm))) 355 (t (mydefint (hackit (1+ exponent) rest parm) 356 (createname parm (- -1 exponent)) 357 (createname parm (- exponent)) parm)))) 358 359(defun afixsign (funct signswitch) 360 ;;;MULTIPLIES FUNCT BY -1 IF SIGNSWITCH IS NIL 361 (cond (signswitch funct) 362 (t (simptimes (list '(mtimes) -1 funct) 1 t)))) 363 364(defun lapshift (fun rest parm) 365 (cond ((atom fun) (merror "LAPSHIFT: expected a cons, not ~M" fun)) 366 ((or (member 'laplace (car fun) :test #'eq) (null rest)) 367 (lapdefint (cond (rest (simptimes (cons '(mtimes) 368 (cons fun rest)) 1 t)) 369 (t fun)) parm)) 370 (t (laptimes (append rest 371 (ncons (cons (append (car fun) 372 '(laplace)) 373 (cdr fun)))) parm)))) 374 375;;;COMPUTES %E**(W*B*%I)*F(S-W*A*%I) WHERE W=-1 IF SIGN IS T ELSE W=1 376(defun mostpart (f parm sign a b) 377 (let ((substinfun ($at f 378 (list '(mequal) 379 parm 380 (list '(mplus) parm (afixsign (list '(mtimes) a '$%i) sign)))))) 381 (if (zerop1 b) 382 substinfun 383 (list '(mtimes) 384 (exponentiate (afixsign (list '(mtimes) b '$%i) (null sign))) 385 substinfun)))) 386 387 ;;;IF WHICHSIGN IS NIL THEN SIN TRANSFORM ELSE COS TRANSFORM 388(defun compose (fun parm whichsign a b) 389 (let ((result (list '((rat) 1 2) 390 (list '(mplus) 391 (mostpart fun parm t a b) 392 (afixsign (mostpart fun parm nil a b) 393 whichsign))))) 394 (sratsimp (simptimes (cons '(mtimes) 395 (if whichsign 396 result 397 (cons '$%i result))) 398 1 nil)))) 399 400 ;;;FUN IS OF THE FORM SIN(A*T+B)*REST(T) OR COS 401(defun lapsin (fun rest trigswitch parm) 402 (let ((ab (islinear (cadr fun) var))) 403 (cond (ab 404 (cond (rest 405 (compose (laptimes rest parm) 406 parm 407 trigswitch 408 (car ab) 409 (cdr ab))) 410 (t 411 (simptimes 412 (list '(mtimes) 413 (cond ((zerop1 (cdr ab)) 414 (if trigswitch parm (car ab))) 415 (t 416 (cond (trigswitch 417 (list '(mplus) 418 (list '(mtimes) 419 parm 420 (list '(%cos) (cdr ab))) 421 (list '(mtimes) 422 -1 423 (car ab) 424 (list '(%sin) (cdr ab))))) 425 (t 426 (list '(mplus) 427 (list '(mtimes) 428 parm 429 (list '(%sin) (cdr ab))) 430 (list '(mtimes) 431 (car ab) 432 (list '(%cos) (cdr ab)))))))) 433 (list '(mexpt) 434 (list '(mplus) 435 (list '(mexpt) parm 2) 436 (list '(mexpt) (car ab) 2)) 437 -1)) 438 1 nil)))) 439 (t 440 (lapshift fun rest parm))))) 441 442 ;;;FUN IS OF THE FORM SINH(A*T+B)*REST(T) OR IS COSH 443(defun lapsinh (fun rest switch parm) 444 (cond ((islinear (cadr fun) var) 445 (sratsimp 446 (laplus 447 (simplus 448 (list '(mplus) 449 (nconc (list '(mtimes) 450 (list '(mexpt) 451 '$%e 452 (cadr fun)) 453 '((rat) 1 2)) 454 rest) 455 (afixsign (nconc (list '(mtimes) 456 (list '(mexpt) 457 '$%e 458 (afixsign (cadr fun) 459 nil)) 460 '((rat) 1 2)) 461 rest) 462 switch)) 463 1 464 nil) parm))) 465 (t (lapshift fun rest parm)))) 466 467 ;;;FUN IS OF THE FORM LOG(A*T) 468(defun laplog (fun parm) 469 (let ((ab (islinear (cadr fun) var))) 470 (cond ((and ab (zerop1 (cdr ab))) 471 (simptimes (list '(mtimes) 472 (list '(mplus) 473 (subfunmake '$psi '(0) (ncons 1)) 474 (list '(%log) (car ab)) 475 (list '(mtimes) -1 (list '(%log) parm))) 476 (list '(mexpt) parm -1)) 477 1 nil)) 478 (t 479 (lapdefint fun parm))))) 480 481(defun raiseup (fbase exponent) 482 (if (equal exponent 1) 483 fbase 484 (list '(mexpt) fbase exponent))) 485 486;;TAKES TRANSFORM OF DELTA(A*T+B)*F(T) 487(defun lapdelta (fun rest parm) 488 (let ((ab (islinear (cadr fun) var)) 489 (sign nil) 490 (recipa nil)) 491 (cond (ab 492 (setq recipa (power (car ab) -1) ab (div (cdr ab) (car ab))) 493 (setq sign (asksign ab) recipa (simplifya (list '(mabs) recipa) nil)) 494 (simplifya (cond ((eq sign '$positive) 495 0) 496 ((eq sign '$zero) 497 (list '(mtimes) 498 (maxima-substitute 0 var (fixuprest rest)) 499 recipa)) 500 (t 501 (list '(mtimes) 502 (maxima-substitute (neg ab) var (fixuprest rest)) 503 (list '(mexpt) '$%e (cons '(mtimes) (cons parm (ncons ab)))) 504 recipa))) 505 nil)) 506 (t 507 (lapshift fun rest parm))))) 508 509(defun laperf (fun parm) 510 (let ((ab (islinear (cadr fun) var))) 511 (cond ((and ab (equal (cdr ab) 0)) 512 (simptimes (list '(mtimes) 513 (div* (exponentiate (div* (list '(mexpt) parm 2) 514 (list '(mtimes) 515 4 516 (list '(mexpt) (car ab) 2)))) 517 parm) 518 (list '(mplus) 519 1 520 (list '(mtimes) 521 -1 522 (list '(%erf) (div* parm (list '(mtimes) 2 (car ab))))))) 523 1 524 nil)) 525 (t 526 (lapdefint fun parm))))) 527 528(defun lapdefint (fun parm) 529 (prog (tryint mult) 530 (and ($unknown fun)(go skip)) 531 (setq mult (simptimes (list '(mtimes) (exponentiate 532 (list '(mtimes) -1 var parm)) fun) 1 nil)) 533 (with-new-context (context) 534 (progn 535 (meval `(($assume) ,@(list (list '(mgreaterp) parm 0)))) 536 (setq tryint 537 ;; $defint should not throw a Maxima error. 538 ;; therefore we set the flags errcatch and errormsg. 539 ;; errset catches an error and returns nil. 540 (let ((errcatch t) ($errormsg nil)) 541 (errset ($defint mult var 0 '$inf)))))) 542 (and tryint (not (eq (and (listp (car tryint)) 543 (caaar tryint)) 544 '%integrate)) 545 (return (car tryint))) 546 skip (return (list '(%laplace) fun var parm)))) 547 548 549(defun lapdiff (fun parm) 550;;;FUN IS OF THE FORM DIFF(F(T),T,N) WHERE N IS A POSITIVE INTEGER 551 (prog (difflist degree frontend resultlist newdlist order arg2) 552 (setq newdlist (setq difflist (copy-tree (cddr fun)))) 553 (setq arg2 (list '(mequal) var 0)) 554 a (cond ((null difflist) 555 (return (cons '(%derivative) 556 (cons (list '(%laplace) 557 (cadr fun) 558 var 559 parm) 560 newdlist)))) 561 ((eq (car difflist) var) 562 (setq degree (cadr difflist) 563 difflist (cddr difflist)) 564 (go out))) 565 (setq difflist (cdr (setq frontend (cdr difflist)))) 566 (go a) 567 out (cond ((null (posint degree)) 568 (return (list '(%laplace) fun var parm)))) 569 (cond (frontend (rplacd frontend difflist)) 570 (t (setq newdlist difflist))) 571 (cond (newdlist (setq fun (cons '(%derivative) 572 (cons (cadr fun) 573 newdlist)))) 574 (t (setq fun (cadr fun)))) 575 (setq order 0) 576 loop (decf degree) 577 (setq resultlist 578 (cons (list '(mtimes) 579 (raiseup parm degree) 580 ($at ($diff fun var order) arg2)) 581 resultlist)) 582 (incf order) 583 (and (> degree 0) (go loop)) 584 (setq resultlist (cond ((cdr resultlist) 585 (cons '(mplus) 586 resultlist)) 587 (t (car resultlist)))) 588 (return (simplus (list '(mplus) 589 (list '(mtimes) 590 (raiseup parm order) 591 (laplace fun parm)) 592 (list '(mtimes) 593 -1 594 resultlist)) 595 1 nil)))) 596 597 ;;;FUN IS OF THE FORM INTEGRATE(F(X)*G(T)*H(T-X),X,0,T) 598(defun lapint (fun parm dvar) 599 (prog (newfun parm-list f var-list var-parm-list) 600 (and dvar (go convolution)) 601 (setq dvar (cadr (setq newfun (cdr fun)))) 602 (and (cddr newfun) 603 (zerop1 (caddr newfun)) 604 (eq (cadddr newfun) var) 605 (go convolutiontest)) 606 notcon 607 (setq newfun (cdr fun)) 608 (cond ((cddr newfun) 609 (cond ((and (freeof var (caddr newfun)) 610 (freeof var (cadddr newfun))) 611 (return (list '(%integrate) 612 (laplace (car newfun) parm dvar) 613 dvar 614 (caddr newfun) 615 (cadddr newfun)))) 616 (t (go giveup)))) 617 (t (return (list '(%integrate) 618 (laplace (car newfun) parm dvar) 619 dvar)))) 620 giveup 621 (return (list '(%laplace) fun var parm)) 622 convolutiontest 623 (setq newfun ($factor (car newfun))) 624 (cond ((eq (caar newfun) 'mtimes) 625 (setq f (cadr newfun) newfun (cddr newfun))) 626 (t (setq f newfun newfun nil))) 627 gothrulist 628 (cond ((freeof dvar f) 629 (setq parm-list (cons f parm-list))) 630 ((freeof var f) (setq var-list (cons f var-list))) 631 ((freeof dvar 632 (sratsimp (maxima-substitute (list '(mplus) 633 var 634 dvar) 635 var 636 f))) 637 (setq var-parm-list (cons f var-parm-list))) 638 (t (go notcon))) 639 (cond (newfun (setq f (car newfun) newfun (cdr newfun)) 640 (go gothrulist))) 641 (and 642 parm-list 643 (return 644 (laplace 645 (cons 646 '(mtimes) 647 (nconc parm-list 648 (ncons (list '(%integrate) 649 (cons '(mtimes) 650 (append var-list 651 var-parm-list)) 652 dvar 653 0 654 var)))) parm dvar))) 655 convolution 656 (return 657 (simptimes 658 (list 659 '(mtimes) 660 (laplace ($expand (maxima-substitute var 661 dvar 662 (fixuprest var-list))) parm dvar) 663 (laplace 664 ($expand (maxima-substitute 0 dvar (fixuprest var-parm-list))) parm dvar)) 665 1 666 t)))) 667 668(declare-top (special varlist ratform ils ilt)) 669 670(defmfun $ilt (exp ils ilt) 671 ;;;EXP IS F(S)/G(S) WHERE F AND G ARE POLYNOMIALS IN S AND DEGR(F) < DEGR(G) 672 (let (varlist ($savefactors t) checkfactors $ratfac $keepfloat) 673 ;;; MAKES ILS THE MAIN VARIABLE 674 (setq varlist (list ils)) 675 (newvar exp) 676 (orderpointer varlist) 677 (setq var (caadr (ratrep* ils))) 678 (cond ((and (null (atom exp)) 679 (eq (caar exp) 'mequal)) 680 (list '(mequal) 681 ($ilt (cadr exp) ils ilt) 682 ($ilt (caddr exp) ils ilt))) 683 ((zerop1 exp) 0) 684 ((freeof ils exp) 685 (list '(%ilt) exp ils ilt)) 686 (t (ilt0 exp))))) 687 688(defun maxima-rationalp (le v) 689 (cond ((null le)) 690 ((and (null (atom (car le))) (null (freeof v (car le)))) 691 nil) 692 (t (maxima-rationalp (cdr le) v)))) 693 694 ;;;THIS FUNCTION DOES THE PARTIAL FRACTION DECOMPOSITION 695(defun ilt0 (exp) 696 (prog (wholepart frpart num denom y content real factor 697 apart bpart parnumer ratarg ratform) 698 (and (mplusp exp) 699 (return (simplus (cons '(mplus) 700 (mapcar #'(lambda (f) ($ilt f ils ilt)) (cdr exp))) 1 t))) 701 (and (null (atom exp)) 702 (eq (caar exp) '%laplace) 703 (eq (cadddr exp) ils) 704 (return (cond ((eq (caddr exp) ilt) (cadr exp)) 705 (t (subst ilt 706 (caddr exp) 707 (cadr exp)))))) 708 (setq ratarg (ratrep* exp)) 709 (or (maxima-rationalp varlist ils) 710 (return (list '(%ilt) exp ils ilt))) 711 (setq ratform (car ratarg)) 712 (setq denom (ratdenominator (cdr ratarg))) 713 (setq frpart (pdivide (ratnumerator (cdr ratarg)) denom)) 714 (setq wholepart (car frpart)) 715 (setq frpart (ratqu (cadr frpart) denom)) 716 (cond ((not (zerop1 (car wholepart))) 717 (return (list '(%ilt) exp ils ilt))) 718 ((zerop1 (car frpart)) (return 0))) 719 (setq num (car frpart) denom (cdr frpart)) 720 (setq y (oldcontent denom)) 721 (setq content (car y)) 722 (setq real (cadr y)) 723 (setq factor (pfactor real)) 724 loop (cond ((null (cddr factor)) 725 (setq apart real 726 bpart 1 727 y '((0 . 1) 1 . 1)) 728 (go skip))) 729 (setq apart (pexpt (car factor) (cadr factor))) 730 (setq bpart (car (ratqu real apart))) 731 (setq y (bprog apart bpart)) 732 skip (setq frpart 733 (cdr (ratdivide (ratti (ratnumerator num) 734 (cdr y) 735 t) 736 (ratti (ratdenominator num) 737 (ratti content apart t) 738 t)))) 739 (setq 740 parnumer 741 (cons (ilt1 (ratqu (ratnumerator frpart) 742 (ratti (ratdenominator frpart) 743 (ratti (ratdenominator num) 744 content 745 t) 746 t)) 747 (car factor) 748 (cadr factor)) 749 parnumer)) 750 (setq factor (cddr factor)) 751 (cond ((null factor) 752 (return (simplus (cons '(mplus) parnumer) 1 t)))) 753 (setq num (cdr (ratdivide (ratti num (car y) t) 754 (ratti content bpart t)))) 755 (setq real bpart) 756 (go loop))) 757 758(declare-top (special q z)) 759 760(defun ilt1 (p q k) 761 (let (z) 762 (cond ((onep1 k) (ilt3 p )) 763 (t (setq z (bprog q (pderivative q var))) 764 (ilt2 p k))))) 765 766 767 ;;;INVERTS P(S)/Q(S)**K WHERE Q(S) IS IRREDUCIBLE 768 ;;;DOESN'T CALL ILT3 IF Q(S) IS LINEAR 769(defun ilt2 (p k) 770 (prog (y a b) 771 (and (onep1 k) (return (ilt3 p))) 772 (decf k) 773 (setq a (ratti p (car z) t)) 774 (setq b (ratti p (cdr z) t)) 775 (setq y (pexpt q k)) 776 (cond 777 ((or (null (equal (pdegree q var) 1)) 778 (> (pdegree (car p) var) 0)) 779 (return 780 (simplus 781 (list 782 '(mplus) 783 (ilt2 784 (cdr (ratdivide (ratplus a (ratqu (ratderivative b var) k)) y)) 785 k) 786 ($multthru (simptimes (list '(mtimes) 787 ilt 788 (power k -1) 789 (ilt2 (cdr (ratdivide b y)) k)) 790 1 791 t))) 792 1 793 t)))) 794 (setq a (disrep (polcoef q 1)) 795 b (disrep (polcoef q 0))) 796 (return 797 (simptimes (list '(mtimes) 798 (disrep p) 799 (raiseup ilt k) 800 (simpexpt (list '(mexpt) 801 '$%e 802 (list '(mtimes) 803 -1 804 ilt 805 b 806 (list '(mexpt) a -1))) 807 1 808 nil) 809 (list '(mexpt) 810 a 811 (- -1 k)) 812 (list '(mexpt) 813 (factorial k) 814 -1)) 815 1 816 nil)))) 817 818(declare-top(notype k)) 819 820;;(DEFUN COEF MACRO (POL) (SUBST (CADR POL) (QUOTE DEG) 821;; '(DISREP (RATQU (POLCOEF (CAR P) DEG) (CDR P))))) 822 823(defmacro coef (pol) 824 `(disrep (ratqu (polcoef (car p) ,pol) (cdr p)))) 825 826(defun lapsum (&rest args) 827 (cons '(mplus) args)) 828 829(defun lapprod (&rest args) 830 (cons '(mtimes) args)) 831 832(defun expo (&rest args) 833 (cons '(mexpt) args)) 834 835;;;INVERTS P(S)/Q(S) WHERE Q(S) IS IRREDUCIBLE 836(defun ilt3 (p) 837 (prog (discrim sign a c d e b1 b0 r term1 term2 degr) 838 (setq e (disrep (polcoef q 0)) 839 d (disrep (polcoef q 1)) 840 degr (pdegree q var)) 841 (and (equal degr 1) 842 (return 843 (simptimes (lapprod 844 (disrep p) 845 (expo d -1) 846 (expo '$%e (lapprod -1 ilt e (expo d -1)))) 847 1 848 nil))) 849 (setq c (disrep (polcoef q 2))) 850 (and (equal degr 2) (go quadratic)) 851 (and (equal degr 3) (zerop1 c) (zerop1 d) 852 (go cubic)) 853 (return (list '(%ilt) (div* (disrep p)(disrep q)) ils ilt)) 854 cubic (setq a (disrep (polcoef q 3)) 855 r (simpnrt (div* e a) 3)) 856 (setq d (div* (disrep p)(lapprod a (lapsum 857 (expo ils 3)(expo '%r 3))))) 858 (return (ilt0 (maxima-substitute r '%r ($partfrac d ils)))) 859 quadratic (setq b0 (coef 0) b1 (coef 1)) 860 861 (setq discrim 862 (simplus (lapsum 863 (lapprod 4 e c) 864 (lapprod -1 d d)) 865 1 866 nil)) 867 (setq sign (cond ((free discrim '$%i) (asksign discrim)) (t '$positive)) 868 term1 '(%cos) 869 term2 '(%sin)) 870 (setq degr (expo '$%e (lapprod ilt d (power c -1) '((rat) -1 2)))) 871 (cond ((eq sign '$zero) 872 (return (simptimes (lapprod degr (lapsum (div* b1 c) 873 (lapprod 874 (div* (lapsum (lapprod 2 b0 c) (lapprod -1 b1 d)) 875 (lapprod 2 c c)) ilt))) 1 nil)) 876 ) ((eq sign '$negative) 877 (setq term1 '(%cosh) 878 term2 '(%sinh) 879 discrim (simptimes (lapprod -1 discrim) 1 t)))) 880 (setq discrim (simpnrt discrim 2)) 881 (setq sign 882 (simptimes 883 (lapprod 884 (lapsum 885 (lapprod 2 b0 c) 886 (lapprod -1 b1 d)) 887 (expo discrim -1)) 888 1 889 nil)) 890 (setq c (power c -1)) 891 (setq discrim (simptimes (lapprod 892 discrim 893 ilt 894 '((rat) 1 2) 895 c) 896 1 897 t)) 898 (return 899 (simptimes 900 (lapprod c degr 901 (lapsum 902 (lapprod b1 (list term1 discrim)) 903 (lapprod sign (list term2 discrim)))) 904 1 905 nil)))) 906 907(declare-top (unspecial ils ilt *nounl* q ratform var 908 varlist z)) 909