Add some basic letsimp tests based on bug #3950
[maxima.git] / src / laplac.lisp
blob20d08e9f049725961239d48896ca0199412257d8
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 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
11 (in-package :maxima)
13 (macsyma-module laplac)
15 (declare-top (special var $savefactors
16 checkfactors $ratfac $keepfloat *nounl* *nounsflag*
17 errcatch $errormsg))
19 ;;; The properties NOUN and VERB give correct linear display
21 (defprop $laplace %laplace verb)
22 (defprop %laplace $laplace noun)
24 (defprop $ilt %ilt verb)
25 (defprop %ilt $ilt noun)
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))))
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))))
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))))))
46 ;;;INITIALIZES SOME GLOBAL VARIABLES THEN CALLS THE DISPATCHING FUNCTION
48 (defmfun $laplace (fun var parm)
49 (setq fun (mratcheck fun))
50 (cond ((or *nounsflag* (member '%laplace *nounl* :test #'eq))
51 (setq fun (remlaplace fun))))
52 (laplace fun parm))
54 ;;;LAMBDA BINDS SOME SPECIAL VARIABLES TO NIL AND DISPATCHES
56 (defun remlaplace (e)
57 (if (atom e)
59 (cons (delete 'laplace (append (car e) nil) :count 1 :test #'eq)
60 (mapcar #'remlaplace (cdr e)))))
62 (defun laplace (fun parm &optional (dvar nil))
63 (let ()
64 ;;; Handles easy cases and calls appropriate function on others.
65 (cond ((mbagp fun) (cons (car fun) (mapcar #'(lambda (e) (laplace e parm)) (cdr fun))))
66 ((equal fun 0) 0)
67 ((equal fun 1)
68 (cond ((zerop1 parm) (simplify (list '($delta) 0)))
69 (t (power parm -1))))
70 ((alike1 fun var) (power parm -2))
71 ((or (atom fun) (freeof var fun))
72 (cond ((zerop1 parm) (mul2 fun (simplify (list '($delta) 0))))
73 (t (mul2 fun (power parm -1)))))
74 (t
75 (let ((op (caar fun)))
76 (let ((result ; We store the result of laplace for further work.
77 (cond ((eq op 'mplus)
78 (laplus fun parm))
79 ((eq op 'mtimes)
80 (laptimes (cdr fun) parm))
81 ((eq op 'mexpt)
82 (lapexpt fun nil parm))
83 ((eq op '%sin)
84 (lapsin fun nil nil parm))
85 ((eq op '%cos)
86 (lapsin fun nil t parm))
87 ((eq op '%sinh)
88 (lapsinh fun nil nil parm))
89 ((eq op '%cosh)
90 (lapsinh fun nil t parm))
91 ((eq op '%log)
92 (laplog fun parm))
93 ((eq op '%derivative)
94 (lapdiff fun parm))
95 ((eq op '%integrate)
96 (lapint fun parm dvar))
97 ((eq op '%sum)
98 (list '(%sum)
99 (laplace (cadr fun) parm)
100 (caddr fun)
101 (cadddr fun)
102 (car (cddddr fun))))
103 ((eq op '%erf)
104 (laperf fun parm))
105 ((and (eq op '%ilt)(eq (cadddr fun) var))
106 (cond ((eq parm (caddr fun))(cadr fun))
107 (t (subst parm (caddr fun)(cadr fun)))))
108 ((eq op '$delta)
109 (lapdelta fun nil parm))
110 ((member op '(%hstep $unit_step))
111 (laphstep fun nil parm))
112 ((setq op ($get op '$laplace))
113 (mcall op fun var parm))
114 (t (lapdefint fun parm)))))
115 (when (isinop result '%integrate)
116 ;; Laplace has not found a result but returns a definit
117 ;; integral. This integral can contain internal integration
118 ;; variables. Replace such a result with the noun form.
119 (setq result (list '(%laplace) fun var parm)))
120 ;; Check if we have a result, when not call $specint.
121 (check-call-to-$specint result fun parm)))))))
123 ;;; Check if laplace has found a result, when not try $specint.
125 (defun check-call-to-$specint (result fun parm)
126 (cond
127 ((or (isinop result '%laplace)
128 (isinop result '%limit) ; Try $specint for incomplete results
129 (isinop result '%at)) ; which contain %limit or %at too.
130 ;; laplace returns a noun form or a result which contains %limit or %at.
131 ;; We pass the function to $specint to look for more results.
132 (let (res)
133 ;; laplace assumes the parameter s to be positive and does a
134 ;; declaration before an integration is done. Therefore we declare
135 ;; the parameter of the Laplace transform to be positive before
136 ;; we call $specint too.
137 (with-new-context (context)
138 (meval `(($assume) ,@(list (list '(mgreaterp) parm 0))))
139 (setq res ($specint (mul fun (power '$%e (mul -1 var parm))) var)))
140 (if (or (isinop res '%specint) ; Both symobls are possible, that is
141 (isinop res '$specint)) ; not consistent! Check it! 02/2009
142 ;; $specint has not found a result.
143 result
144 ;; $specint has found a result
145 res)))
146 (t result)))
148 (defun laplus (fun parm)
149 (simplus (cons '(mplus) (mapcar #'(lambda (e) (laplace e parm)) (cdr fun))) 1 t))
151 (defun laptimes (fun parm)
152 ;;;EXPECTS A LIST (PERHAPS EMPTY) OF FUNCTIONS MULTIPLIED TOGETHER WITHOUT THE MTIMES
153 ;;;SEES IF IT CAN APPLY THE FIRST AS A TRANSFORMATION ON THE REST OF THE FUNCTIONS
154 (cond ((null fun) (list '(mexpt) parm -1.))
155 ((null (cdr fun)) (laplace (car fun) parm))
156 ((freeof var (car fun))
157 (simptimes (list '(mtimes)
158 (car fun)
159 (laptimes (cdr fun) parm))
162 ((eq (car fun) var)
163 (simptimes (list '(mtimes) -1 (sdiff (laptimes (cdr fun) parm) parm))
167 (let ((op (caaar fun)))
168 (cond ((eq op 'mexpt)
169 (lapexpt (car fun) (cdr fun) parm))
170 ((eq op 'mplus)
171 (laplus ($multthru (fixuprest (cdr fun)) (car fun)) parm))
172 ((eq op '%sin)
173 (lapsin (car fun) (cdr fun) nil parm))
174 ((eq op '%cos)
175 (lapsin (car fun) (cdr fun) t parm))
176 ((eq op '%sinh)
177 (lapsinh (car fun) (cdr fun) nil parm))
178 ((eq op '%cosh)
179 (lapsinh (car fun) (cdr fun) t parm))
180 ((eq op '$delta)
181 (lapdelta (car fun) (cdr fun) parm))
182 ((member op '(%hstep $unit_step))
183 (laphstep (car fun) (cdr fun) parm))
185 (lapshift (car fun) (cdr fun) parm)))))))
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
197 (isquadraticp
198 (cond ((eq base-of-fun '$%e) power)
199 (t (simptimes (list '(mtimes)
200 power
201 (list '(%log)
202 base-of-fun))
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)
225 (list '(mtimes)
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)
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))
265 (exponentiate (div* (list '(mexpt) parm 2)
266 (list '(mtimes) 4 result)))
267 (list '(mplus)
269 (list '(mtimes)
271 (list '(%erf)
272 (div* parm
273 (list '(mtimes)
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)
285 (cadr ab)))
286 parm
287 result)))
288 (return (simptimes (list '(mtimes)
289 (exponentiate (caddr ab))
290 result) 1 nil))
291 var-case
292 (cond ((or (null rest) (freeof var (fixuprest rest)))
293 (go var-easy-case)))
294 (cond ((posint power)
295 (return (afixsign (apply '$diff
296 (list (laptimes rest parm)
297 parm
298 power))
299 (even power))))
300 ((negint power)
301 (return (mydefint (hackit power rest parm)
302 (createname parm (- power))
303 parm parm)))
304 (t (go noluck)))
305 var-easy-case
306 (setq power
307 (simplus (list '(mplus) 1 power) 1 t))
308 (or (eq (asksign ($realpart power)) '$positive) (go noluck))
309 (setq result (list (list '(%gamma) power)
310 (list '(mexpt)
311 parm
312 (afixsign power nil))))
313 (and rest (setq result (nconc result rest)))
314 (return (simptimes (cons '(mtimes) result) 1 nil))
315 noluck
316 (return
317 (cond
318 ((and (posint power)
319 (member (caar base-of-fun)
320 '(mplus %sin %cos %sinh %cosh) :test #'eq))
321 (laptimes (cons base-of-fun
322 (cons (cond ((= power 2) base-of-fun)
323 (t (list '(mexpt)
324 base-of-fun
325 (1- power))))
326 rest)) parm))
327 (t (lapshift fun rest parm))))))
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 parm)
355 (cond ((equal exponent -1)
356 (let ((parm (createname parm 1)))
357 (laptimes rest parm)))
358 (t (mydefint (hackit (1+ exponent) rest parm)
359 (createname parm (- -1 exponent))
360 (createname parm (- exponent)) parm))))
362 (defun afixsign (funct signswitch)
363 ;;;MULTIPLIES FUNCT BY -1 IF SIGNSWITCH IS NIL
364 (cond (signswitch funct)
365 (t (simptimes (list '(mtimes) -1 funct) 1 t))))
367 (defun lapshift (fun rest parm)
368 (cond ((atom fun) (merror "LAPSHIFT: expected a cons, not ~M" fun))
369 ((or (member 'laplace (car fun) :test #'eq) (null rest))
370 (lapdefint (cond (rest (simptimes (cons '(mtimes)
371 (cons fun rest)) 1 t))
372 (t fun)) parm))
373 (t (laptimes (append rest
374 (ncons (cons (append (car fun)
375 '(laplace))
376 (cdr fun)))) parm))))
378 ;;;COMPUTES %E**(W*B*%I)*F(S-W*A*%I) WHERE W=-1 IF SIGN IS T ELSE W=1
379 (defun mostpart (f parm sign a b)
380 (let ((substinfun ($at f
381 (list '(mequal)
382 parm
383 (list '(mplus) parm (afixsign (list '(mtimes) a '$%i) sign))))))
384 (if (zerop1 b)
385 substinfun
386 (list '(mtimes)
387 (exponentiate (afixsign (list '(mtimes) b '$%i) (null sign)))
388 substinfun))))
390 ;;;IF WHICHSIGN IS NIL THEN SIN TRANSFORM ELSE COS TRANSFORM
391 (defun compose (fun parm whichsign a b)
392 (let ((result (list '((rat) 1 2)
393 (list '(mplus)
394 (mostpart fun parm t a b)
395 (afixsign (mostpart fun parm nil a b)
396 whichsign)))))
397 (sratsimp (simptimes (cons '(mtimes)
398 (if whichsign
399 result
400 (cons '$%i result)))
401 1 nil))))
403 ;;;FUN IS OF THE FORM SIN(A*T+B)*REST(T) OR COS
404 (defun lapsin (fun rest trigswitch parm)
405 (let ((ab (islinear (cadr fun) var)))
406 (cond (ab
407 (cond (rest
408 (compose (laptimes rest parm)
409 parm
410 trigswitch
411 (car ab)
412 (cdr ab)))
414 (simptimes
415 (list '(mtimes)
416 (cond ((zerop1 (cdr ab))
417 (if trigswitch parm (car ab)))
419 (cond (trigswitch
420 (list '(mplus)
421 (list '(mtimes)
422 parm
423 (list '(%cos) (cdr ab)))
424 (list '(mtimes)
426 (car ab)
427 (list '(%sin) (cdr ab)))))
429 (list '(mplus)
430 (list '(mtimes)
431 parm
432 (list '(%sin) (cdr ab)))
433 (list '(mtimes)
434 (car ab)
435 (list '(%cos) (cdr ab))))))))
436 (list '(mexpt)
437 (list '(mplus)
438 (list '(mexpt) parm 2)
439 (list '(mexpt) (car ab) 2))
440 -1))
441 1 nil))))
443 (lapshift fun rest parm)))))
445 ;;;FUN IS OF THE FORM SINH(A*T+B)*REST(T) OR IS COSH
446 (defun lapsinh (fun rest switch parm)
447 (cond ((islinear (cadr fun) var)
448 (sratsimp
449 (laplus
450 (simplus
451 (list '(mplus)
452 (nconc (list '(mtimes)
453 (list '(mexpt)
454 '$%e
455 (cadr fun))
456 '((rat) 1 2))
457 rest)
458 (afixsign (nconc (list '(mtimes)
459 (list '(mexpt)
460 '$%e
461 (afixsign (cadr fun)
462 nil))
463 '((rat) 1 2))
464 rest)
465 switch))
467 nil) parm)))
468 (t (lapshift fun rest parm))))
470 ;;;FUN IS OF THE FORM LOG(A*T)
471 (defun laplog (fun parm)
472 (let ((ab (islinear (cadr fun) var)))
473 (cond ((and ab (zerop1 (cdr ab)))
474 (simptimes (list '(mtimes)
475 (list '(mplus)
476 (subfunmake '$psi '(0) (ncons 1))
477 (list '(%log) (car ab))
478 (list '(mtimes) -1 (list '(%log) parm)))
479 (list '(mexpt) parm -1))
480 1 nil))
482 (lapdefint fun parm)))))
484 (defun raiseup (fbase exponent)
485 (if (equal exponent 1)
486 fbase
487 (list '(mexpt) fbase exponent)))
489 ;;TAKES TRANSFORM OF HSTEP(A*T+B)*F(T)
490 (defun laphstep (fun rest parm)
491 (let ((ab (islinear (cadr fun) var)))
492 (if ab
493 (let* ((a (car ab))
494 (b (cdr ab))
495 (offset (div b a)))
496 (case (asksign a)
497 ($positive (if (eq (asksign offset) '$negative)
498 (mul (exponentiate (mul offset parm))
499 (laplace (maxima-substitute
500 (sub var offset) var
501 (fixuprest rest)) parm))
502 (laptimes rest parm)))
503 ($negative (if (eq (asksign offset) '$negative)
504 (sub (laptimes rest parm)
505 (mul (exponentiate (mul offset parm))
506 (laplace (maxima-substitute
507 (sub var offset) var
508 (fixuprest rest)) parm)))
510 (t (mul fun (laptimes rest parm)))))
511 (lapshift fun rest parm))))
513 ;;TAKES TRANSFORM OF DELTA(A*T+B)*F(T)
514 (defun lapdelta (fun rest parm)
515 (let ((ab (islinear (cadr fun) var))
516 (sign nil)
517 (recipa nil))
518 (cond (ab
519 (setq recipa (power (car ab) -1) ab (div (cdr ab) (car ab)))
520 (setq sign (asksign ab) recipa (simplifya (list '(mabs) recipa) nil))
521 (simplifya (cond ((eq sign '$positive)
523 ((eq sign '$zero)
524 (list '(mtimes)
525 (maxima-substitute 0 var (fixuprest rest))
526 recipa))
528 (list '(mtimes)
529 (maxima-substitute (neg ab) var (fixuprest rest))
530 (list '(mexpt) '$%e (cons '(mtimes) (cons parm (ncons ab))))
531 recipa)))
532 nil))
534 (lapshift fun rest parm)))))
536 (defun laperf (fun parm)
537 (let ((ab (islinear (cadr fun) var)))
538 (cond ((and ab (equal (cdr ab) 0))
539 (simptimes (list '(mtimes)
540 (div* (exponentiate (div* (list '(mexpt) parm 2)
541 (list '(mtimes)
543 (list '(mexpt) (car ab) 2))))
544 parm)
545 (list '(mplus)
547 (list '(mtimes)
549 (list '(%erf) (div* parm (list '(mtimes) 2 (car ab)))))))
551 nil))
553 (lapdefint fun parm)))))
555 (defun lapdefint (fun parm)
556 (prog (tryint mult)
557 (and ($unknown fun)(go skip))
558 (setq mult (simptimes (list '(mtimes) (exponentiate
559 (list '(mtimes) -1 var parm)) fun) 1 nil))
560 (with-new-context (context)
561 (meval `(($assume) ,@(list (list '(mgreaterp) parm 0))))
562 (setq tryint
563 ;; We don't want any errors thrown by $defint to propagate,
564 ;; so we set errset, errcatch, $errormsg and *mdebug*. If
565 ;; $defint throws a error, then no error message is printed
566 ;; and errset simply returns nil below.
567 (let ((errset nil)
568 (errcatch t)
569 ($errormsg nil)
570 (*mdebug* nil))
571 (errset ($defint mult var 0 '$inf)))))
572 (and tryint (not (eq (and (listp (car tryint))
573 (caaar tryint))
574 '%integrate))
575 (return (car tryint)))
576 skip (return (list '(%laplace) fun var parm))))
579 (defun lapdiff (fun parm)
580 ;;;FUN IS OF THE FORM DIFF(F(T),T,N) WHERE N IS A POSITIVE INTEGER
581 (prog (difflist degree frontend resultlist newdlist order arg2)
582 (setq newdlist (setq difflist (copy-tree (cddr fun))))
583 (setq arg2 (list '(mequal) var 0))
584 a (cond ((null difflist)
585 (return (cons '(%derivative)
586 (cons (list '(%laplace)
587 (cadr fun)
589 parm)
590 newdlist))))
591 ((eq (car difflist) var)
592 (setq degree (cadr difflist)
593 difflist (cddr difflist))
594 (go out)))
595 (setq difflist (cdr (setq frontend (cdr difflist))))
596 (go a)
597 out (cond ((null (posint degree))
598 (return (list '(%laplace) fun var parm))))
599 (cond (frontend (rplacd frontend difflist))
600 (t (setq newdlist difflist)))
601 (cond (newdlist (setq fun (cons '(%derivative)
602 (cons (cadr fun)
603 newdlist))))
604 (t (setq fun (cadr fun))))
605 (setq order 0)
606 loop (decf degree)
607 (setq resultlist
608 (cons (list '(mtimes)
609 (raiseup parm degree)
610 ($at ($diff fun var order) arg2))
611 resultlist))
612 (incf order)
613 (and (> degree 0) (go loop))
614 (setq resultlist (cond ((cdr resultlist)
615 (cons '(mplus)
616 resultlist))
617 (t (car resultlist))))
618 (return (simplus (list '(mplus)
619 (list '(mtimes)
620 (raiseup parm order)
621 (laplace fun parm))
622 (list '(mtimes)
624 resultlist))
625 1 nil))))
627 ;;;FUN IS OF THE FORM INTEGRATE(F(X)*G(T)*H(T-X),X,0,T)
628 (defun lapint (fun parm dvar)
629 (prog (newfun parm-list f var-list var-parm-list)
630 (and dvar (go convolution))
631 (setq dvar (cadr (setq newfun (cdr fun))))
632 (and (cddr newfun)
633 (zerop1 (caddr newfun))
634 (eq (cadddr newfun) var)
635 (go convolutiontest))
636 notcon
637 (setq newfun (cdr fun))
638 (cond ((cddr newfun)
639 (cond ((and (freeof var (caddr newfun))
640 (freeof var (cadddr newfun)))
641 (return (list '(%integrate)
642 (laplace (car newfun) parm dvar)
643 dvar
644 (caddr newfun)
645 (cadddr newfun))))
646 (t (go giveup))))
647 (t (return (list '(%integrate)
648 (laplace (car newfun) parm dvar)
649 dvar))))
650 giveup
651 (return (list '(%laplace) fun var parm))
652 convolutiontest
653 (setq newfun ($factor (car newfun)))
654 (cond ((eq (caar newfun) 'mtimes)
655 (setq f (cadr newfun) newfun (cddr newfun)))
656 (t (setq f newfun newfun nil)))
657 gothrulist
658 (cond ((freeof dvar f)
659 (setq parm-list (cons f parm-list)))
660 ((freeof var f) (setq var-list (cons f var-list)))
661 ((freeof dvar
662 (sratsimp (maxima-substitute (list '(mplus)
664 dvar)
666 f)))
667 (setq var-parm-list (cons f var-parm-list)))
668 (t (go notcon)))
669 (cond (newfun (setq f (car newfun) newfun (cdr newfun))
670 (go gothrulist)))
671 (and
672 parm-list
673 (return
674 (laplace
675 (cons
676 '(mtimes)
677 (nconc parm-list
678 (ncons (list '(%integrate)
679 (cons '(mtimes)
680 (append var-list
681 var-parm-list))
682 dvar
684 var)))) parm dvar)))
685 convolution
686 (return
687 (simptimes
688 (list
689 '(mtimes)
690 (laplace ($expand (maxima-substitute var
691 dvar
692 (fixuprest var-list))) parm dvar)
693 (laplace
694 ($expand (maxima-substitute 0 dvar (fixuprest var-parm-list))) parm dvar))
696 t))))
698 (declare-top (special varlist ratform ils ilt))
700 (defmfun $ilt (exp ils ilt)
701 ;;;EXP IS F(S)/G(S) WHERE F AND G ARE POLYNOMIALS IN S AND DEGR(F) < DEGR(G)
702 (let (varlist ($savefactors t) checkfactors $ratfac $keepfloat)
703 ;;; MAKES ILS THE MAIN VARIABLE
704 (setq varlist (list ils))
705 (newvar exp)
706 (orderpointer varlist)
707 (setq var (caadr (ratrep* ils)))
708 (cond ((and (null (atom exp))
709 (eq (caar exp) 'mequal))
710 (list '(mequal)
711 ($ilt (cadr exp) ils ilt)
712 ($ilt (caddr exp) ils ilt)))
713 ((zerop1 exp) 0)
714 ((freeof ils exp)
715 (list '(mtimes) exp (list '($delta) ilt)))
716 (t (ilt0 exp)))))
718 (defun maxima-rationalp (le v)
719 (cond ((null le))
720 ((and (null (atom (car le))) (null (freeof v (car le))))
721 nil)
722 (t (maxima-rationalp (cdr le) v))))
724 ;;;THIS FUNCTION DOES THE PARTIAL FRACTION DECOMPOSITION
725 (defun ilt0 (exp)
726 (prog (wholepart frpart num denom y content real factor
727 apart bpart parnumer ratarg ratform)
728 (and (mplusp exp)
729 (return (simplus (cons '(mplus)
730 (mapcar #'(lambda (f) ($ilt f ils ilt)) (cdr exp))) 1 t)))
731 (and (null (atom exp))
732 (eq (caar exp) '%laplace)
733 (eq (cadddr exp) ils)
734 (return (cond ((eq (caddr exp) ilt) (cadr exp))
735 (t (subst ilt
736 (caddr exp)
737 (cadr exp))))))
738 (setq ratarg (ratrep* exp))
739 (or (maxima-rationalp varlist ils)
740 (return (list '(%ilt) exp ils ilt)))
741 (setq ratform (car ratarg))
742 (setq denom (ratdenominator (cdr ratarg)))
743 (setq frpart (pdivide (ratnumerator (cdr ratarg)) denom))
744 (setq wholepart (car frpart))
745 (setq frpart (ratqu (cadr frpart) denom))
746 (cond ((not (zerop1 (car wholepart)))
747 (return (list '(%ilt) exp ils ilt)))
748 ((zerop1 (car frpart)) (return 0)))
749 (setq num (car frpart) denom (cdr frpart))
750 (setq y (oldcontent denom))
751 (setq content (car y))
752 (setq real (cadr y))
753 (setq factor (pfactor real))
754 loop (cond ((null (cddr factor))
755 (setq apart real
756 bpart 1
757 y '((0 . 1) 1 . 1))
758 (go skip)))
759 (setq apart (pexpt (car factor) (cadr factor)))
760 (setq bpart (car (ratqu real apart)))
761 (setq y (bprog apart bpart))
762 skip (setq frpart
763 (cdr (ratdivide (ratti (ratnumerator num)
764 (cdr y)
766 (ratti (ratdenominator num)
767 (ratti content apart t)
768 t))))
769 (setq
770 parnumer
771 (cons (ilt1 (ratqu (ratnumerator frpart)
772 (ratti (ratdenominator frpart)
773 (ratti (ratdenominator num)
774 content
777 (car factor)
778 (cadr factor))
779 parnumer))
780 (setq factor (cddr factor))
781 (cond ((null factor)
782 (return (simplus (cons '(mplus) parnumer) 1 t))))
783 (setq num (cdr (ratdivide (ratti num (car y) t)
784 (ratti content bpart t))))
785 (setq real bpart)
786 (go loop)))
788 (declare-top (special q z))
790 (defun ilt1 (p q k)
791 (let (z)
792 (cond ((onep1 k) (ilt3 p ))
793 (t (setq z (bprog q (pderivative q var)))
794 (ilt2 p k)))))
797 ;;;INVERTS P(S)/Q(S)**K WHERE Q(S) IS IRREDUCIBLE
798 ;;;DOESN'T CALL ILT3 IF Q(S) IS LINEAR
799 (defun ilt2 (p k)
800 (prog (y a b)
801 (and (onep1 k) (return (ilt3 p)))
802 (decf k)
803 (setq a (ratti p (car z) t))
804 (setq b (ratti p (cdr z) t))
805 (setq y (pexpt q k))
806 (cond
807 ((or (null (equal (pdegree q var) 1))
808 (> (pdegree (car p) var) 0))
809 (return
810 (simplus
811 (list
812 '(mplus)
813 (ilt2
814 (cdr (ratdivide (ratplus a (ratqu (ratderivative b var) k)) y))
816 ($multthru (simptimes (list '(mtimes)
818 (power k -1)
819 (ilt2 (cdr (ratdivide b y)) k))
821 t)))
823 t))))
824 (setq a (disrep (polcoef q 1))
825 b (disrep (polcoef q 0)))
826 (return
827 (simptimes (list '(mtimes)
828 (disrep p)
829 (raiseup ilt k)
830 (simpexpt (list '(mexpt)
831 '$%e
832 (list '(mtimes)
836 (list '(mexpt) a -1)))
838 nil)
839 (list '(mexpt)
841 (- -1 k))
842 (list '(mexpt)
843 (factorial k)
844 -1))
846 nil))))
848 ;;(DEFUN COEF MACRO (POL) (SUBST (CADR POL) (QUOTE DEG)
849 ;; '(DISREP (RATQU (POLCOEF (CAR P) DEG) (CDR P)))))
851 (defmacro coef (pol)
852 `(disrep (ratqu (polcoef (car p) ,pol) (cdr p))))
854 (defun lapsum (&rest args)
855 (cons '(mplus) args))
857 (defun lapprod (&rest args)
858 (cons '(mtimes) args))
860 (defun expo (&rest args)
861 (cons '(mexpt) args))
863 ;;;INVERTS P(S)/Q(S) WHERE Q(S) IS IRREDUCIBLE
864 (defun ilt3 (p)
865 (prog (discrim sign a c d e b1 b0 r term1 term2 degr)
866 (setq e (disrep (polcoef q 0))
867 d (disrep (polcoef q 1))
868 degr (pdegree q var))
869 (and (equal degr 1)
870 (return
871 (simptimes (lapprod
872 (disrep p)
873 (expo d -1)
874 (expo '$%e (lapprod -1 ilt e (expo d -1))))
876 nil)))
877 (setq c (disrep (polcoef q 2)))
878 (and (equal degr 2) (go quadratic))
879 (and (equal degr 3) (zerop1 c) (zerop1 d)
880 (go cubic))
881 (return (list '(%ilt) (div* (disrep p)(disrep q)) ils ilt))
882 cubic (setq a (disrep (polcoef q 3))
883 r (simpnrt (div* e a) 3))
884 (setq d (div* (disrep p)(lapprod a (lapsum
885 (expo ils 3)(expo '%r 3)))))
886 (return (ilt0 (maxima-substitute r '%r ($partfrac d ils))))
887 quadratic (setq b0 (coef 0) b1 (coef 1))
889 (setq discrim
890 (simplus (lapsum
891 (lapprod 4 e c)
892 (lapprod -1 d d))
894 nil))
895 (setq sign (cond ((free discrim '$%i) (asksign discrim)) (t '$positive))
896 term1 '(%cos)
897 term2 '(%sin))
898 (setq degr (expo '$%e (lapprod ilt d (power c -1) '((rat) -1 2))))
899 (cond ((eq sign '$zero)
900 (return (simptimes (lapprod degr (lapsum (div* b1 c)
901 (lapprod
902 (div* (lapsum (lapprod 2 b0 c) (lapprod -1 b1 d))
903 (lapprod 2 c c)) ilt))) 1 nil))
904 ) ((eq sign '$negative)
905 (setq term1 '(%cosh)
906 term2 '(%sinh)
907 discrim (simptimes (lapprod -1 discrim) 1 t))))
908 (setq discrim (simpnrt discrim 2))
909 (setq sign
910 (simptimes
911 (lapprod
912 (lapsum
913 (lapprod 2 b0 c)
914 (lapprod -1 b1 d))
915 (expo discrim -1))
917 nil))
918 (setq c (power c -1))
919 (setq discrim (simptimes (lapprod
920 discrim
922 '((rat) 1 2)
926 (return
927 (simptimes
928 (lapprod c degr
929 (lapsum
930 (lapprod b1 (list term1 discrim))
931 (lapprod sign (list term2 discrim))))
933 nil))))
935 (declare-top (unspecial ils ilt *nounl* q ratform var
936 varlist z))