1 ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; The data in this file contains enhancements. ;;;;;
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 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
13 (macsyma-module laplac
)
15 (declare-top (special context
))
17 ;;; The properties NOUN and VERB give correct linear display
19 (defprop $laplace %laplace verb
)
20 (defprop %laplace $laplace noun
)
22 (defprop $ilt %ilt verb
)
23 (defprop %ilt $ilt noun
)
25 (defun exponentiate (pow)
26 ;;COMPUTES %E**Z WHERE Z IS AN ARBITRARY EXPRESSION TAKING SOME OF THE WORK AWAY FROM SIMPEXPT
27 (cond ((zerop1 pow
) 1)
29 (t (power '$%e pow
))))
31 (defun fixuprest (rest)
32 ;;REST IS A PRODUCT WITHOUT THE MTIMES.FIXUPREST PUTS BACK THE MTIMES
34 ((cdr rest
) (cons '(mtimes) rest
))
37 (defun isquadraticp (e x
)
38 (let ((b (sdiff e x
)))
39 (cond ((zerop1 b
) (list 0 0 e
))
40 ((freeof x b
) (list 0 b
(maxima-substitute 0 x e
)))
41 ((setq b
(islinear b x
))
42 (list (div* (car b
) 2) (cdr b
) (maxima-substitute 0 x e
))))))
44 ;;;INITIALIZES SOME GLOBAL VARIABLES THEN CALLS THE DISPATCHING FUNCTION
46 (defmfun $laplace
(fun time-var parm
)
47 (setq fun
(mratcheck fun
))
48 (cond ((or *nounsflag
* (member '%laplace
*nounl
* :test
#'eq
))
49 (setq fun
(remlaplace fun
))))
50 (laplace fun time-var parm
))
52 ;;;LAMBDA BINDS SOME SPECIAL VARIABLES TO NIL AND DISPATCHES
57 (cons (delete 'laplace
(append (car e
) nil
) :count
1 :test
#'eq
)
58 (mapcar #'remlaplace
(cdr e
)))))
60 (defun laplace (fun time-var parm
&optional
(dvar nil
))
62 ;; Handles easy cases and calls appropriate function on others.
63 (cond ((mbagp fun
) (cons (car fun
) (mapcar #'(lambda (e) (laplace e time-var parm
)) (cdr fun
))))
66 (cond ((zerop1 parm
) (simplify (list '($delta
) 0)))
68 ((alike1 fun time-var
) (power parm -
2))
69 ((or (atom fun
) (freeof time-var fun
))
70 (cond ((zerop1 parm
) (mul2 fun
(simplify (list '($delta
) 0))))
71 (t (mul2 fun
(power parm -
1)))))
73 (let ((op (caar fun
)))
74 (let ((result ; We store the result of laplace for further work.
76 (laplus fun time-var parm
))
78 (laptimes (cdr fun
) time-var parm
))
80 (lapexpt fun nil time-var parm
))
82 (lapsin fun nil nil time-var parm
))
84 (lapsin fun nil t time-var parm
))
86 (lapsinh fun nil nil time-var parm
))
88 (lapsinh fun nil t time-var parm
))
90 (laplog fun time-var parm
))
92 (lapdiff fun time-var parm
))
94 (lapint fun time-var parm dvar
))
97 (laplace (cadr fun
) time-var parm
)
102 (laperf fun time-var parm
))
103 ((and (eq op
'%ilt
)(eq (cadddr fun
) time-var
))
104 (cond ((eq parm
(caddr fun
))(cadr fun
))
105 (t (subst parm
(caddr fun
)(cadr fun
)))))
107 (lapdelta fun nil time-var parm
))
108 ((member op
'(%hstep $unit_step
))
109 (laphstep fun nil time-var parm
))
110 ((setq op
($get op
'$laplace
))
111 (mcall op fun time-var parm
))
113 (lapdefint fun time-var parm
)))))
114 (when (isinop result
'%integrate
)
115 ;; Laplace has not found a result but returns a definite
116 ;; integral. This integral can contain internal integration
117 ;; variables. Replace such a result with the noun form.
118 (setq result
(list '(%laplace
) fun time-var parm
)))
119 ;; Check if we have a result, when not call $specint.
120 (check-call-to-$specint result fun time-var parm
)))))))
122 ;;; Check if laplace has found a result, when not try $specint.
124 (defun check-call-to-$specint
(result fun time-var parm
)
126 ((or (isinop result
'%laplace
)
127 (isinop result
'%limit
) ; Try $specint for incomplete results
128 (isinop result
'%at
)) ; which contain %limit or %at too.
129 ;; laplace returns a noun form or a result which contains %limit or %at.
130 ;; We pass the function to $specint to look for more results.
132 ;; laplace assumes the parameter s to be positive and does a
133 ;; declaration before an integration is done. Therefore we declare
134 ;; the parameter of the Laplace transform to be positive before
135 ;; we call $specint too.
136 (with-new-context (context)
137 (meval `(($assume
) ,@(list (list '(mgreaterp) parm
0))))
138 (setq res
($specint
(mul fun
(power '$%e
(mul -
1 time-var parm
))) time-var
)))
139 (if (or (isinop res
'%specint
) ; Both symbols are possible, that is
140 (isinop res
'$specint
)) ; not consistent! Check it! 02/2009
141 ;; $specint has not found a result.
143 ;; $specint has found a result
147 (defun laplus (fun time-var parm
)
148 (simplus (cons '(mplus) (mapcar #'(lambda (e) (laplace e time-var parm
)) (cdr fun
))) 1 t
))
150 (defun laptimes (fun time-var parm
)
151 ;;EXPECTS A LIST (PERHAPS EMPTY) OF FUNCTIONS MULTIPLIED TOGETHER WITHOUT THE MTIMES
152 ;;SEES IF IT CAN APPLY THE FIRST AS A TRANSFORMATION ON THE REST OF THE FUNCTIONS
153 (cond ((null fun
) (list '(mexpt) parm -
1.
))
154 ((null (cdr fun
)) (laplace (car fun
) time-var parm
))
155 ((freeof time-var
(car fun
))
156 (simptimes (list '(mtimes)
158 (laptimes (cdr fun
) time-var parm
))
161 ((eq (car fun
) time-var
)
162 (simptimes (list '(mtimes) -
1 (sdiff (laptimes (cdr fun
) time-var parm
) parm
))
166 (let ((op (caaar fun
)))
167 (cond ((eq op
'mexpt
)
168 (lapexpt (car fun
) (cdr fun
) time-var parm
))
170 (laplus ($multthru
(fixuprest (cdr fun
)) (car fun
)) time-var parm
))
172 (lapsin (car fun
) (cdr fun
) nil time-var parm
))
174 (lapsin (car fun
) (cdr fun
) t time-var parm
))
176 (lapsinh (car fun
) (cdr fun
) nil time-var parm
))
178 (lapsinh (car fun
) (cdr fun
) t time-var parm
))
180 (lapdelta (car fun
) (cdr fun
) time-var parm
))
181 ((member op
'(%hstep $unit_step
))
182 (laphstep (car fun
) (cdr fun
) time-var parm
))
184 (lapshift (car fun
) (cdr fun
) time-var parm
)))))))
186 (defun lapexpt (fun rest time-var 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
))
193 (freeof time-var base-of-fun
)
197 (cond ((eq base-of-fun
'$%e
) power
)
198 (t (simptimes (list '(mtimes)
205 (cond ((equal (car ab
) 0) (go %e-case-lin
))
206 ((null rest
) (go %e-case-quad
))
208 ((and (eq base-of-fun time-var
) (freeof time-var power
))
210 ((and (alike1 '((rat) -
1 2) power
) (null rest
)
211 (setq ab
(islinear base-of-fun time-var
)))
212 (setq result
(div* (cdr ab
) (car ab
)))
221 (exponentiate (list '(mtimes) result parm
))
237 (rest (sratsimp ($at
(laptimes rest time-var parm
)
247 (afixsign (cadr ab
) nil
))
249 (return (simptimes (list '(mtimes)
250 (exponentiate (caddr ab
))
255 (setq result
(afixsign (car ab
) nil
))
264 (exponentiate (div* (list '(mexpt) parm
2)
265 (list '(mtimes) 4 result
)))
278 (and (null (equal (cadr ab
) 0))
280 (maxima-substitute (list '(mplus)
287 (return (simptimes (list '(mtimes)
288 (exponentiate (caddr ab
))
291 (cond ((or (null rest
) (freeof time-var
(fixuprest rest
)))
293 (cond ((posint power
)
294 (return (afixsign (apply '$diff
295 (list (laptimes rest time-var parm
)
300 (return (mydefint (hackit power rest time-var parm
)
301 (createname parm
(- power
))
306 (simplus (list '(mplus) 1 power
) 1 t
))
307 (or (eq (asksign ($realpart power
)) '$positive
) (go noluck
))
308 (setq result
(list (list '(%gamma
) power
)
311 (afixsign power nil
))))
312 (and rest
(setq result
(nconc result rest
)))
313 (return (simptimes (cons '(mtimes) result
) 1 nil
))
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
)
327 (t (lapshift fun rest time-var parm
))))))
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))))
344 (errset ($defint f x a
'$inf
)))))))
347 (list '(%integrate
) f x a
'$inf
))))
349 ;;;CREATES UNIQUE NAMES FOR VARIABLE OF INTEGRATION
350 (defun createname (head tail
)
351 (intern (format nil
"~S~S" head tail
)))
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 time-var parm
)
355 (cond ((equal exponent -
1)
356 (let ((parm (createname parm
1)))
357 (laptimes rest time-var parm
)))
359 (mydefint (hackit (1+ exponent
) rest time-var parm
)
360 (createname parm
(- -
1 exponent
))
361 (createname parm
(- exponent
)) parm
))))
363 (defun afixsign (funct signswitch
)
364 ;;MULTIPLIES FUNCT BY -1 IF SIGNSWITCH IS NIL
365 (cond (signswitch funct
)
366 (t (simptimes (list '(mtimes) -
1 funct
) 1 t
))))
368 (defun lapshift (fun rest time-var parm
)
370 (merror "LAPSHIFT: expected a cons, not ~M" fun
))
371 ((or (member 'laplace
(car fun
) :test
#'eq
)
373 (lapdefint (cond (rest (simptimes (cons '(mtimes)
374 (cons fun rest
)) 1 t
))
378 (laptimes (append rest
379 (ncons (cons (append (car fun
)
384 ;;;COMPUTES %E**(W*B*%I)*F(S-W*A*%I) WHERE W=-1 IF SIGN IS T ELSE W=1
385 (defun mostpart (f parm sign a b
)
386 (let ((substinfun ($at f
389 (list '(mplus) parm
(afixsign (list '(mtimes) a
'$%i
) sign
))))))
393 (exponentiate (afixsign (list '(mtimes) b
'$%i
) (null sign
)))
396 ;;;IF WHICHSIGN IS NIL THEN SIN TRANSFORM ELSE COS TRANSFORM
397 (defun compose (fun parm whichsign a b
)
398 (let ((result (list '((rat) 1 2)
400 (mostpart fun parm t a b
)
401 (afixsign (mostpart fun parm nil a b
)
403 (sratsimp (simptimes (cons '(mtimes)
409 ;;;FUN IS OF THE FORM SIN(A*T+B)*REST(T) OR COS
410 (defun lapsin (fun rest trigswitch time-var parm
)
411 (let ((ab (islinear (cadr fun
) time-var
)))
414 (compose (laptimes rest time-var parm
)
422 (cond ((zerop1 (cdr ab
))
423 (if trigswitch parm
(car ab
)))
429 (list '(%cos
) (cdr ab
)))
433 (list '(%sin
) (cdr ab
)))))
438 (list '(%sin
) (cdr ab
)))
441 (list '(%cos
) (cdr ab
))))))))
444 (list '(mexpt) parm
2)
445 (list '(mexpt) (car ab
) 2))
449 (lapshift fun rest time-var parm
)))))
451 ;;;FUN IS OF THE FORM SINH(A*T+B)*REST(T) OR IS COSH
452 (defun lapsinh (fun rest switch time-var parm
)
453 (cond ((islinear (cadr fun
) time-var
)
458 (nconc (list '(mtimes)
464 (afixsign (nconc (list '(mtimes)
476 (lapshift fun rest time-var parm
))))
478 ;;;FUN IS OF THE FORM LOG(A*T)
479 (defun laplog (fun time-var parm
)
480 (let ((ab (islinear (cadr fun
) time-var
)))
481 (cond ((and ab
(zerop1 (cdr ab
)))
482 (simptimes (list '(mtimes)
484 (subfunmake '$psi
'(0) (ncons 1))
485 (list '(%log
) (car ab
))
486 (list '(mtimes) -
1 (list '(%log
) parm
)))
487 (list '(mexpt) parm -
1))
490 (lapdefint fun time-var parm
)))))
492 (defun raiseup (fbase exponent
)
493 (if (equal exponent
1)
495 (list '(mexpt) fbase exponent
)))
497 ;;TAKES TRANSFORM OF HSTEP(A*T+B)*F(T)
498 (defun laphstep (fun rest time-var parm
)
499 (let ((ab (islinear (cadr fun
) time-var
)))
506 (if (eq (asksign offset
) '$negative
)
507 (mul (exponentiate (mul offset parm
))
508 (laplace (maxima-substitute
509 (sub time-var offset
) time-var
512 (laptimes rest time-var parm
)))
514 (if (eq (asksign offset
) '$negative
)
515 (sub (laptimes rest time-var parm
)
516 (mul (exponentiate (mul offset parm
))
517 (laplace (maxima-substitute
518 (sub time-var offset
) time-var
523 (mul fun
(laptimes rest time-var parm
)))))
524 (lapshift fun rest time-var parm
))))
526 ;;TAKES TRANSFORM OF DELTA(A*T+B)*F(T)
527 (defun lapdelta (fun rest time-var parm
)
528 (let ((ab (islinear (cadr fun
) time-var
))
532 (setq recipa
(power (car ab
) -
1) ab
(div (cdr ab
) (car ab
)))
533 (setq sign
(asksign ab
) recipa
(simplifya (list '(mabs) recipa
) nil
))
534 (simplifya (cond ((eq sign
'$positive
)
538 (maxima-substitute 0 time-var
(fixuprest rest
))
542 (maxima-substitute (neg ab
) time-var
(fixuprest rest
))
543 (list '(mexpt) '$%e
(cons '(mtimes) (cons parm
(ncons ab
))))
547 (lapshift fun rest time-var parm
)))))
549 (defun laperf (fun time-var parm
)
550 (let ((ab (islinear (cadr fun
) time-var
)))
551 (cond ((and ab
(equal (cdr ab
) 0))
552 (simptimes (list '(mtimes)
553 (div* (exponentiate (div* (list '(mexpt) parm
2)
556 (list '(mexpt) (car ab
) 2))))
562 (list '(%erf
) (div* parm
(list '(mtimes) 2 (car ab
)))))))
566 (lapdefint fun time-var parm
)))))
568 (defun lapdefint (fun time-var parm
)
570 (and ($unknown fun
)(go skip
))
571 (setq mult
(simptimes (list '(mtimes) (exponentiate
572 (list '(mtimes) -
1 time-var parm
)) fun
) 1 nil
))
573 (with-new-context (context)
574 (meval `(($assume
) ,@(list (list '(mgreaterp) parm
0))))
576 ;; We don't want any errors thrown by $defint to propagate,
577 ;; so we set errset, errcatch, $errormsg and *mdebug*. If
578 ;; $defint throws a error, then no error message is printed
579 ;; and errset simply returns nil below.
584 (errset ($defint mult time-var
0 '$inf
)))))
585 (and tryint
(not (eq (and (listp (car tryint
))
588 (return (car tryint
)))
590 (return (list '(%laplace
) fun time-var parm
))))
593 (defun lapdiff (fun time-var parm
)
594 ;;FUN IS OF THE FORM DIFF(F(T),T,N) WHERE N IS A POSITIVE INTEGER
595 (prog (difflist degree frontend resultlist newdlist order arg2
)
596 (setq newdlist
(setq difflist
(copy-tree (cddr fun
))))
597 (setq arg2
(list '(mequal) time-var
0))
599 (cond ((null difflist
)
600 (return (cons '(%derivative
)
601 (cons (list '(%laplace
)
606 ((eq (car difflist
) time-var
)
607 (setq degree
(cadr difflist
)
608 difflist
(cddr difflist
))
610 (setq difflist
(cdr (setq frontend
(cdr difflist
))))
613 (cond ((null (posint degree
))
614 (return (list '(%laplace
) fun time-var parm
))))
615 (cond (frontend (rplacd frontend difflist
))
616 (t (setq newdlist difflist
)))
617 (cond (newdlist (setq fun
(cons '(%derivative
)
620 (t (setq fun
(cadr fun
))))
625 (cons (list '(mtimes)
626 (raiseup parm degree
)
627 ($at
($diff fun time-var order
) arg2
))
630 (and (> degree
0) (go loop
))
631 (setq resultlist
(cond ((cdr resultlist
)
634 (t (car resultlist
))))
635 (return (simplus (list '(mplus)
638 (laplace fun time-var parm
))
644 ;;;FUN IS OF THE FORM INTEGRATE(F(X)*G(T)*H(T-X),X,0,T)
645 (defun lapint (fun time-var parm dvar
)
646 (prog (newfun parm-list f var-list var-parm-list
)
647 (and dvar
(go convolution
))
648 (setq dvar
(cadr (setq newfun
(cdr fun
))))
650 (zerop1 (caddr newfun
))
651 (eq (cadddr newfun
) time-var
)
652 (go convolutiontest
))
654 (setq newfun
(cdr fun
))
656 (cond ((and (freeof time-var
(caddr newfun
))
657 (freeof time-var
(cadddr newfun
)))
658 (return (list '(%integrate
)
659 (laplace (car newfun
) time-var parm dvar
)
664 (t (return (list '(%integrate
)
665 (laplace (car newfun
) time-var parm dvar
)
668 (return (list '(%laplace
) fun time-var parm
))
670 (setq newfun
($factor
(car newfun
)))
671 (cond ((eq (caar newfun
) 'mtimes
)
672 (setq f
(cadr newfun
) newfun
(cddr newfun
)))
673 (t (setq f newfun newfun nil
)))
675 (cond ((freeof dvar f
)
676 (setq parm-list
(cons f parm-list
)))
677 ((freeof time-var f
) (setq var-list
(cons f var-list
)))
679 (sratsimp (maxima-substitute (list '(mplus)
684 (setq var-parm-list
(cons f var-parm-list
)))
687 (setq f
(car newfun
) newfun
(cdr newfun
))
696 (ncons (list '(%integrate
)
703 time-var parm dvar
)))
709 (laplace ($expand
(maxima-substitute time-var
711 (fixuprest var-list
)))
714 ($expand
(maxima-substitute 0 dvar
(fixuprest var-parm-list
)))
719 (defmfun $ilt
(exp ils ilt
)
720 ;;EXP IS F(S)/G(S) WHERE F AND G ARE POLYNOMIALS IN S AND DEGR(F) < DEGR(G)
721 (let (varlist ($savefactors t
) *checkfactors
* $ratfac $keepfloat
723 ;; MAKES ILS THE MAIN VARIABLE
724 (setq varlist
(list ils
))
726 (orderpointer varlist
)
727 (setq s-var
(caadr (ratrep* ils
)))
728 (cond ((and (null (atom exp
))
729 (eq (caar exp
) 'mequal
))
731 ($ilt
(cadr exp
) ils ilt
)
732 ($ilt
(caddr exp
) ils ilt
)))
735 (list '(mtimes) exp
(list '($delta
) ilt
)))
736 (t (ilt0 exp ils ilt s-var
)))))
738 (defun maxima-rationalp (le v
)
740 ((and (null (atom (car le
))) (null (freeof v
(car le
))))
743 (maxima-rationalp (cdr le
) v
))))
745 ;;;THIS FUNCTION DOES THE PARTIAL FRACTION DECOMPOSITION
746 (defun ilt0 (exp ils ilt s-var
)
747 (prog (wholepart frpart num denom y content real factor
748 apart bpart parnumer ratarg laplac-ratform
)
750 (return (simplus (cons '(mplus)
751 (mapcar #'(lambda (f) ($ilt f ils ilt
)) (cdr exp
))) 1 t
)))
752 (and (null (atom exp
))
753 (eq (caar exp
) '%laplace
)
754 (eq (cadddr exp
) ils
)
755 (return (cond ((eq (caddr exp
) ilt
) (cadr exp
))
759 (setq ratarg
(ratrep* exp
))
760 (or (maxima-rationalp varlist ils
)
761 (return (list '(%ilt
) exp ils ilt
)))
762 (setq laplac-ratform
(car ratarg
))
763 (setq denom
(ratdenominator (cdr ratarg
)))
764 (setq frpart
(pdivide (ratnumerator (cdr ratarg
)) denom
))
765 (setq wholepart
(car frpart
))
766 (setq frpart
(ratqu (cadr frpart
) denom
))
767 (cond ((not (zerop1 (car wholepart
)))
768 (return (list '(%ilt
) exp ils ilt
)))
769 ((zerop1 (car frpart
)) (return 0)))
770 (setq num
(car frpart
) denom
(cdr frpart
))
771 (setq y
(oldcontent denom
))
772 (setq content
(car y
))
774 (setq factor
(pfactor real
))
776 (cond ((null (cddr factor
))
781 (setq apart
(pexpt (car factor
) (cadr factor
)))
782 (setq bpart
(car (ratqu real apart
)))
783 (setq y
(bprog apart bpart s-var
))
786 (cdr (ratdivide (ratti (ratnumerator num
)
789 (ratti (ratdenominator num
)
790 (ratti content apart t
)
794 (cons (ilt1 (ratqu (ratnumerator frpart
)
795 (ratti (ratdenominator frpart
)
796 (ratti (ratdenominator num
)
807 (setq factor
(cddr factor
))
809 (return (simplus (cons '(mplus) parnumer
) 1 t
))))
810 (setq num
(cdr (ratdivide (ratti num
(car y
) t
)
811 (ratti content bpart t
))))
815 (defun ilt1 (p q k laplac-ratform ils ilt s-var
)
818 (ilt3 p laplac-ratform ils ilt q s-var
))
820 (setq z
(bprog q
(pderivative q s-var
) s-var
))
821 (ilt2 p k laplac-ratform ils ilt z q s-var
)))))
824 ;;;INVERTS P(S)/Q(S)**K WHERE Q(S) IS IRREDUCIBLE
825 ;;;DOESN'T CALL ILT3 IF Q(S) IS LINEAR
826 (defun ilt2 (p k laplac-ratform ils ilt z q s-var
)
828 (and (onep1 k
) (return (ilt3 p laplac-ratform ils ilt q s-var
)))
830 (setq a
(ratti p
(car z
) t
))
831 (setq b
(ratti p
(cdr z
) t
))
834 ((or (null (equal (pdegree q s-var
) 1))
835 (> (pdegree (car p
) s-var
) 0))
841 (cdr (ratdivide (ratplus a
(ratqu (ratderivative b s-var
) k
)) y
))
849 ($multthru
(simptimes (list '(mtimes)
852 (ilt2 (cdr (ratdivide b y
))
864 (setq a
(disrep (polcoef q
1 s-var
) laplac-ratform
)
865 b
(disrep (polcoef q
0 s-var
) laplac-ratform
))
867 (simptimes (list '(mtimes)
868 (disrep p laplac-ratform
)
870 (simpexpt (list '(mexpt)
876 (list '(mexpt) a -
1)))
888 ;;(DEFUN COEF MACRO (POL) (SUBST (CADR POL) (QUOTE DEG)
889 ;; '(DISREP (RATQU (POLCOEF (CAR P) DEG) (CDR P)))))
891 ;;;INVERTS P(S)/Q(S) WHERE Q(S) IS IRREDUCIBLE
892 (defun ilt3 (p laplac-ratform ils ilt q s-var
)
894 ((coef (pol coef-ratform
)
895 (disrep (ratqu (polcoef (car p
) pol s-var
) (cdr p
)) coef-ratform
))
896 ;; FIXME: Consider replacing these 3 functions with
897 ;; make-mplus-l, make-mtimes-l, and make-expt-l.
899 (cons '(mplus) args
))
900 (lapprod (&rest args
)
901 (cons '(mtimes) args
))
903 (cons '(mexpt) args
)))
904 (prog (discrim sign a c d e b1 b0 r term1 term2 degr
)
905 (setq e
(disrep (polcoef q
0 s-var
) laplac-ratform
)
906 d
(disrep (polcoef q
1 s-var
) laplac-ratform
)
907 degr
(pdegree q s-var
))
910 (simptimes (lapprod (disrep p laplac-ratform
)
912 (expo '$%e
(lapprod -
1 ilt e
(expo d -
1))))
915 (setq c
(disrep (polcoef q
2 s-var
) laplac-ratform
))
916 (and (equal degr
2) (go quadratic
))
917 (and (equal degr
3) (zerop1 c
) (zerop1 d
)
919 (return (list '(%ilt
) (div* (disrep p laplac-ratform
)
920 (disrep q laplac-ratform
))
923 (setq a
(disrep (polcoef q
3 s-var
) laplac-ratform
)
924 r
(simpnrt (div* e a
) 3))
925 (setq d
(div* (disrep p laplac-ratform
)
926 (lapprod a
(lapsum (expo ils
3)
928 (return (ilt0 (maxima-substitute r
'%r
($partfrac d ils
)) ils ilt s-var
))
930 (setq b0
(coef 0 laplac-ratform
) b1
(coef 1 laplac-ratform
))
938 (setq sign
(cond ((free discrim
'$%i
) (asksign discrim
)) (t '$positive
))
941 (setq degr
(expo '$%e
(lapprod ilt d
(power c -
1) '((rat) -
1 2))))
942 (cond ((eq sign
'$zero
)
943 (return (simptimes (lapprod degr
946 (div* (lapsum (lapprod 2 b0 c
)
951 ((eq sign
'$negative
)
954 discrim
(simptimes (lapprod -
1 discrim
) 1 t
))))
955 (setq discrim
(simpnrt discrim
2))
965 (setq c
(power c -
1))
966 (setq discrim
(simptimes (lapprod discrim
976 (lapprod b1
(list term1 discrim
))
977 (lapprod sign
(list term2 discrim
))))