Import package raddenest by Gilles Schintgen, adapted from corresponding code in...
[maxima.git] / src / laplac.lisp
blobd43c7f01098742835d2aa1070edd0e3f67d934ed
1 ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; The data in this file contains enhancements. ;;;;;
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 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
11 (in-package :maxima)
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)
28 ((equal pow 1) '$%e)
29 (t (power '$%e pow))))
31 (defun fixuprest (rest)
32 ;;REST IS A PRODUCT WITHOUT THE MTIMES.FIXUPREST PUTS BACK THE MTIMES
33 (cond ((null rest) 1)
34 ((cdr rest) (cons '(mtimes) rest))
35 (t (car 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
54 (defun remlaplace (e)
55 (if (atom e)
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))
61 (let ()
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))))
64 ((equal fun 0) 0)
65 ((equal fun 1)
66 (cond ((zerop1 parm) (simplify (list '($delta) 0)))
67 (t (power parm -1))))
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)))))
72 (t
73 (let ((op (caar fun)))
74 (let ((result ; We store the result of laplace for further work.
75 (cond ((eq op 'mplus)
76 (laplus fun time-var parm))
77 ((eq op 'mtimes)
78 (laptimes (cdr fun) time-var parm))
79 ((eq op 'mexpt)
80 (lapexpt fun nil time-var parm))
81 ((eq op '%sin)
82 (lapsin fun nil nil time-var parm))
83 ((eq op '%cos)
84 (lapsin fun nil t time-var parm))
85 ((eq op '%sinh)
86 (lapsinh fun nil nil time-var parm))
87 ((eq op '%cosh)
88 (lapsinh fun nil t time-var parm))
89 ((eq op '%log)
90 (laplog fun time-var parm))
91 ((eq op '%derivative)
92 (lapdiff fun time-var parm))
93 ((eq op '%integrate)
94 (lapint fun time-var parm dvar))
95 ((eq op '%sum)
96 (list '(%sum)
97 (laplace (cadr fun) time-var parm)
98 (caddr fun)
99 (cadddr fun)
100 (car (cddddr fun))))
101 ((eq op '%erf)
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)))))
106 ((eq op '$delta)
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 definit
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)
125 (cond
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.
131 (let (res)
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 symobls are possible, that is
140 (isinop res '$specint)) ; not consistent! Check it! 02/2009
141 ;; $specint has not found a result.
142 result
143 ;; $specint has found a result
144 res)))
145 (t 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)
157 (car fun)
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))
169 ((eq op 'mplus)
170 (laplus ($multthru (fixuprest (cdr fun)) (car fun)) time-var parm))
171 ((eq op '%sin)
172 (lapsin (car fun) (cdr fun) nil time-var parm))
173 ((eq op '%cos)
174 (lapsin (car fun) (cdr fun) t time-var parm))
175 ((eq op '%sinh)
176 (lapsinh (car fun) (cdr fun) nil time-var parm))
177 ((eq op '%cosh)
178 (lapsinh (car fun) (cdr fun) t time-var parm))
179 ((eq op '$delta)
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))
191 (cond
192 ((and
193 (freeof time-var base-of-fun)
194 (setq
196 (isquadraticp
197 (cond ((eq base-of-fun '$%e) power)
198 (t (simptimes (list '(mtimes)
199 power
200 (list '(%log)
201 base-of-fun))
203 nil)))
204 time-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 time-var) (freeof time-var power))
209 (go var-case))
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)))
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)
224 (list '(mtimes)
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 time-var 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)
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))
264 (exponentiate (div* (list '(mexpt) parm 2)
265 (list '(mtimes) 4 result)))
266 (list '(mplus)
268 (list '(mtimes)
270 (list '(%erf)
271 (div* parm
272 (list '(mtimes)
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)
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 time-var (fixuprest rest)))
292 (go var-easy-case)))
293 (cond ((posint power)
294 (return (afixsign (apply '$diff
295 (list (laptimes rest time-var parm)
296 parm
297 power))
298 (even power))))
299 ((negint power)
300 (return (mydefint (hackit power rest time-var 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 ($realpart 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))
326 time-var parm))
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))))
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))))
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 parm time-var)))
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)
369 (cond ((atom fun)
370 (merror "LAPSHIFT: expected a cons, not ~M" fun))
371 ((or (member 'laplace (car fun) :test #'eq)
372 (null rest))
373 (lapdefint (cond (rest (simptimes (cons '(mtimes)
374 (cons fun rest)) 1 t))
375 (t fun))
376 time-var parm))
378 (laptimes (append rest
379 (ncons (cons (append (car fun)
380 '(laplace))
381 (cdr fun))))
382 time-var parm))))
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
387 (list '(mequal)
388 parm
389 (list '(mplus) parm (afixsign (list '(mtimes) a '$%i) sign))))))
390 (if (zerop1 b)
391 substinfun
392 (list '(mtimes)
393 (exponentiate (afixsign (list '(mtimes) b '$%i) (null sign)))
394 substinfun))))
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)
399 (list '(mplus)
400 (mostpart fun parm t a b)
401 (afixsign (mostpart fun parm nil a b)
402 whichsign)))))
403 (sratsimp (simptimes (cons '(mtimes)
404 (if whichsign
405 result
406 (cons '$%i result)))
407 1 nil))))
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)))
412 (cond (ab
413 (cond (rest
414 (compose (laptimes rest parm time-var)
415 parm
416 trigswitch
417 (car ab)
418 (cdr ab)))
420 (simptimes
421 (list '(mtimes)
422 (cond ((zerop1 (cdr ab))
423 (if trigswitch parm (car ab)))
425 (cond (trigswitch
426 (list '(mplus)
427 (list '(mtimes)
428 parm
429 (list '(%cos) (cdr ab)))
430 (list '(mtimes)
432 (car ab)
433 (list '(%sin) (cdr ab)))))
435 (list '(mplus)
436 (list '(mtimes)
437 parm
438 (list '(%sin) (cdr ab)))
439 (list '(mtimes)
440 (car ab)
441 (list '(%cos) (cdr ab))))))))
442 (list '(mexpt)
443 (list '(mplus)
444 (list '(mexpt) parm 2)
445 (list '(mexpt) (car ab) 2))
446 -1))
447 1 nil))))
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)
454 (sratsimp
455 (laplus
456 (simplus
457 (list '(mplus)
458 (nconc (list '(mtimes)
459 (list '(mexpt)
460 '$%e
461 (cadr fun))
462 '((rat) 1 2))
463 rest)
464 (afixsign (nconc (list '(mtimes)
465 (list '(mexpt)
466 '$%e
467 (afixsign (cadr fun)
468 nil))
469 '((rat) 1 2))
470 rest)
471 switch))
473 nil)
474 time-var parm)))
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)
483 (list '(mplus)
484 (subfunmake '$psi '(0) (ncons 1))
485 (list '(%log) (car ab))
486 (list '(mtimes) -1 (list '(%log) parm)))
487 (list '(mexpt) parm -1))
488 1 nil))
490 (lapdefint fun time-var parm)))))
492 (defun raiseup (fbase exponent)
493 (if (equal exponent 1)
494 fbase
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)))
500 (if ab
501 (let* ((a (car ab))
502 (b (cdr ab))
503 (offset (div b a)))
504 (case (asksign a)
505 ($positive
506 (if (eq (asksign offset) '$negative)
507 (mul (exponentiate (mul offset parm))
508 (laplace (maxima-substitute
509 (sub time-var offset) time-var
510 (fixuprest rest))
511 time-var parm))
512 (laptimes rest time-var parm)))
513 ($negative
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
519 (fixuprest rest))
520 time-var parm)))
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))
529 (sign nil)
530 (recipa nil))
531 (cond (ab
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)
536 ((eq sign '$zero)
537 (list '(mtimes)
538 (maxima-substitute 0 time-var (fixuprest rest))
539 recipa))
541 (list '(mtimes)
542 (maxima-substitute (neg ab) time-var (fixuprest rest))
543 (list '(mexpt) '$%e (cons '(mtimes) (cons parm (ncons ab))))
544 recipa)))
545 nil))
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)
554 (list '(mtimes)
556 (list '(mexpt) (car ab) 2))))
557 parm)
558 (list '(mplus)
560 (list '(mtimes)
562 (list '(%erf) (div* parm (list '(mtimes) 2 (car ab)))))))
564 nil))
566 (lapdefint fun time-var parm)))))
568 (defun lapdefint (fun time-var parm)
569 (prog (tryint mult)
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))))
575 (setq tryint
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.
580 (let ((errset nil)
581 (errcatch t)
582 ($errormsg nil)
583 (*mdebug* nil))
584 (errset ($defint mult time-var 0 '$inf)))))
585 (and tryint (not (eq (and (listp (car tryint))
586 (caaar tryint))
587 '%integrate))
588 (return (car tryint)))
589 skip
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)
602 (cadr fun)
603 time-var
604 parm)
605 newdlist))))
606 ((eq (car difflist) time-var)
607 (setq degree (cadr difflist)
608 difflist (cddr difflist))
609 (go out)))
610 (setq difflist (cdr (setq frontend (cdr difflist))))
611 (go a)
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)
618 (cons (cadr fun)
619 newdlist))))
620 (t (setq fun (cadr fun))))
621 (setq order 0)
622 loop
623 (decf degree)
624 (setq resultlist
625 (cons (list '(mtimes)
626 (raiseup parm degree)
627 ($at ($diff fun time-var order) arg2))
628 resultlist))
629 (incf order)
630 (and (> degree 0) (go loop))
631 (setq resultlist (cond ((cdr resultlist)
632 (cons '(mplus)
633 resultlist))
634 (t (car resultlist))))
635 (return (simplus (list '(mplus)
636 (list '(mtimes)
637 (raiseup parm order)
638 (laplace fun time-var parm))
639 (list '(mtimes)
641 resultlist))
642 1 nil))))
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))))
649 (and (cddr newfun)
650 (zerop1 (caddr newfun))
651 (eq (cadddr newfun) time-var)
652 (go convolutiontest))
653 notcon
654 (setq newfun (cdr fun))
655 (cond ((cddr newfun)
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)
660 dvar
661 (caddr newfun)
662 (cadddr newfun))))
663 (t (go giveup))))
664 (t (return (list '(%integrate)
665 (laplace (car newfun) time-var parm dvar)
666 dvar))))
667 giveup
668 (return (list '(%laplace) fun time-var parm))
669 convolutiontest
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)))
674 gothrulist
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)))
678 ((freeof dvar
679 (sratsimp (maxima-substitute (list '(mplus)
680 time-var
681 dvar)
682 time-var
683 f)))
684 (setq var-parm-list (cons f var-parm-list)))
685 (t (go notcon)))
686 (cond (newfun
687 (setq f (car newfun) newfun (cdr newfun))
688 (go gothrulist)))
689 (and
690 parm-list
691 (return
692 (laplace
693 (cons
694 '(mtimes)
695 (nconc parm-list
696 (ncons (list '(%integrate)
697 (cons '(mtimes)
698 (append var-list
699 var-parm-list))
700 dvar
702 time-var))))
703 time-var parm dvar)))
704 convolution
705 (return
706 (simptimes
707 (list
708 '(mtimes)
709 (laplace ($expand (maxima-substitute time-var
710 dvar
711 (fixuprest var-list)))
712 time-var parm dvar)
713 (laplace
714 ($expand (maxima-substitute 0 dvar (fixuprest var-parm-list)))
715 time-var parm dvar))
717 t))))
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
722 s-var)
723 ;; MAKES ILS THE MAIN VARIABLE
724 (setq varlist (list ils))
725 (newvar exp)
726 (orderpointer varlist)
727 (setq s-var (caadr (ratrep* ils)))
728 (cond ((and (null (atom exp))
729 (eq (caar exp) 'mequal))
730 (list '(mequal)
731 ($ilt (cadr exp) ils ilt)
732 ($ilt (caddr exp) ils ilt)))
733 ((zerop1 exp) 0)
734 ((freeof ils exp)
735 (list '(mtimes) exp (list '($delta) ilt)))
736 (t (ilt0 exp ils ilt s-var)))))
738 (defun maxima-rationalp (le v)
739 (cond ((null le))
740 ((and (null (atom (car le))) (null (freeof v (car le))))
741 nil)
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)
749 (and (mplusp exp)
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))
756 (t (subst ilt
757 (caddr exp)
758 (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))
773 (setq real (cadr y))
774 (setq factor (pfactor real))
775 loop
776 (cond ((null (cddr factor))
777 (setq apart real
778 bpart 1
779 y '((0 . 1) 1 . 1))
780 (go skip)))
781 (setq apart (pexpt (car factor) (cadr factor)))
782 (setq bpart (car (ratqu real apart)))
783 (setq y (bprog apart bpart s-var))
784 skip
785 (setq frpart
786 (cdr (ratdivide (ratti (ratnumerator num)
787 (cdr y)
789 (ratti (ratdenominator num)
790 (ratti content apart t)
791 t))))
792 (setq
793 parnumer
794 (cons (ilt1 (ratqu (ratnumerator frpart)
795 (ratti (ratdenominator frpart)
796 (ratti (ratdenominator num)
797 content
800 (car factor)
801 (cadr factor)
802 laplac-ratform
805 s-var)
806 parnumer))
807 (setq factor (cddr factor))
808 (cond ((null factor)
809 (return (simplus (cons '(mplus) parnumer) 1 t))))
810 (setq num (cdr (ratdivide (ratti num (car y) t)
811 (ratti content bpart t))))
812 (setq real bpart)
813 (go loop)))
815 (defun ilt1 (p q k laplac-ratform ils ilt s-var)
816 (let (z)
817 (cond ((onep1 k)
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)
827 (prog (y a b)
828 (and (onep1 k) (return (ilt3 p laplac-ratform ils ilt q s-var)))
829 (decf k)
830 (setq a (ratti p (car z) t))
831 (setq b (ratti p (cdr z) t))
832 (setq y (pexpt q k))
833 (cond
834 ((or (null (equal (pdegree q s-var) 1))
835 (> (pdegree (car p) s-var) 0))
836 (return
837 (simplus
838 (list
839 '(mplus)
840 (ilt2
841 (cdr (ratdivide (ratplus a (ratqu (ratderivative b s-var) k)) y))
843 laplac-ratform
848 s-var)
849 ($multthru (simptimes (list '(mtimes)
851 (power k -1)
852 (ilt2 (cdr (ratdivide b y))
854 laplac-ratform
859 s-var))
861 t)))
863 t))))
864 (setq a (disrep (polcoef q 1 s-var) laplac-ratform)
865 b (disrep (polcoef q 0 s-var) laplac-ratform))
866 (return
867 (simptimes (list '(mtimes)
868 (disrep p laplac-ratform)
869 (raiseup ilt k)
870 (simpexpt (list '(mexpt)
871 '$%e
872 (list '(mtimes)
876 (list '(mexpt) a -1)))
878 nil)
879 (list '(mexpt)
881 (- -1 k))
882 (list '(mexpt)
883 (factorial k)
884 -1))
886 nil))))
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)
893 (labels
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.
898 (lapsum (&rest args)
899 (cons '(mplus) args))
900 (lapprod (&rest args)
901 (cons '(mtimes) args))
902 (expo (&rest 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))
908 (and (equal degr 1)
909 (return
910 (simptimes (lapprod (disrep p laplac-ratform)
911 (expo d -1)
912 (expo '$%e (lapprod -1 ilt e (expo d -1))))
914 nil)))
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)
918 (go cubic))
919 (return (list '(%ilt) (div* (disrep p laplac-ratform)
920 (disrep q laplac-ratform))
921 ils ilt))
922 cubic
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)
927 (expo '%r 3)))))
928 (return (ilt0 (maxima-substitute r '%r ($partfrac d ils)) ils ilt s-var))
929 quadratic
930 (setq b0 (coef 0 laplac-ratform) b1 (coef 1 laplac-ratform))
932 (setq discrim
933 (simplus (lapsum
934 (lapprod 4 e c)
935 (lapprod -1 d d))
937 nil))
938 (setq sign (cond ((free discrim '$%i) (asksign discrim)) (t '$positive))
939 term1 '(%cos)
940 term2 '(%sin))
941 (setq degr (expo '$%e (lapprod ilt d (power c -1) '((rat) -1 2))))
942 (cond ((eq sign '$zero)
943 (return (simptimes (lapprod degr
944 (lapsum (div* b1 c)
945 (lapprod
946 (div* (lapsum (lapprod 2 b0 c)
947 (lapprod -1 b1 d))
948 (lapprod 2 c c))
949 ilt)))
950 1 nil)))
951 ((eq sign '$negative)
952 (setq term1 '(%cosh)
953 term2 '(%sinh)
954 discrim (simptimes (lapprod -1 discrim) 1 t))))
955 (setq discrim (simpnrt discrim 2))
956 (setq sign
957 (simptimes
958 (lapprod
959 (lapsum
960 (lapprod 2 b0 c)
961 (lapprod -1 b1 d))
962 (expo discrim -1))
964 nil))
965 (setq c (power c -1))
966 (setq discrim (simptimes (lapprod discrim
968 '((rat) 1 2)
972 (return
973 (simptimes
974 (lapprod c degr
975 (lapsum
976 (lapprod b1 (list term1 discrim))
977 (lapprod sign (list term2 discrim))))
979 nil)))))