module laplace; % Package for Laplace and inverse Laplace transforms.
% Authors: C. Kazasov, M. Spiridonova, V. Tomov.
% Date: 24 October 1988.
% Revisions:
% 5 Nov 1993 H. Melenk: adapt code for REDUCE 3.5:
% - safe restoration of environment.
% - moved *mcd/*exp:=nil after initial
% simp/reval call for safer pattern
% match
% - enable invlap(1/x^n,x,t) (wrong termination
% condition)
% - repair fctrf call in invlap (incomplete input
% conversion and incomplete result test)
% - repair of pattern matching for rules
% with 2 argument laplace and invlap
% expressions as used in the xmpl file
%
%
% 2 Dec 1988. Commented out rule for sqrt(-x), since it interferes
% with integrator.
% 20 Nov 1988. Converted to lower case and tabs removed.
%
%*******************************************************************
%* *
%* L A P L A C E 2.0 *
%* *
%* AN EXPERIMENTAL PACKAGE FOR PERFORMING IN REDUCE 3 *
%* DIRECT AND INVERSE LAPLACE TRANSFORMATIONS *
%* *
%* SOFIA UNIVERSITY - B U L G A R I A *
%* *
%*******************************************************************
create!-package('(laplace),'(contrib misc));
fluid '(!*exp !*limitedfactors !*mcd !*precise !*rounded depl!* kord!*
subfg!* transvar!* varstack*);
global '(lpsm!* lpcm!* lpshm!* lpchm!* lpse!* lpce!* lpshe!* lpche!*
lpexpt!* ile1!* ile2!* ile3!* ile4!* ile5!* lpvar!* ilvar!*
lpshift!* !*lmsg !*lmon !*ltrig !*lhyp !*ldone !*lione );
switch lhyp,lmon,ltrig;
% Default value:
!*lmsg:= t;
% put('intl,'simpfn,'simpiden);
% put('one, 'simpfn,'simpiden);
% put('delta,'simpfn,'simpiden);
% put('gamma,'simpfn,'simpiden);
if not (gettype 'intl = 'operator) then algebraic operator intl;
if not (gettype 'one = 'operator) then algebraic operator one;
if not (gettype 'delta = 'operator) then algebraic operator delta;
if not (gettype 'gamma = 'operator) then algebraic operator gamma;
%*******************************************************************
%* *
%* Save and restore environment *
%* *
%*******************************************************************
symbolic procedure lap!-save!-environment();
begin scalar u;
u:={ !*exp,!*mcd,kord!*,depl!*,
get('expt,'opmtch),
get('sin,'opmtch),
get('cos,'opmtch),
get('sinh,'opmtch),
get('cosh,'opmtch),
get('gamma,'simpfn),
get('one,'simpfn),
get('delta,'simpfn),
get('intl,'simpfn),
get('laplace,'simpfn),
get('invlap,'simpfn)
};
% copy lists such that rplac* don't touch the environment
kord!* := append(kord!*,nil);
depl!*:=for each d in depl!* collect append(d,nil);
return u;
end;
symbolic procedure lap!-restore!-environment(u);
begin
!*exp := car u; u := cdr u;
!*mcd := car u; u := cdr u;
kord!*:= car u; u := cdr u;
depl!*:= car u; u := cdr u;
put('expt,'opmtch, car u); u:=cdr u;
put('sin,'opmtch, car u); u:=cdr u;
put('cos,'opmtch, car u); u:=cdr u;
put('sinh,'opmtch, car u); u:=cdr u;
put('cosh,'opmtch, car u); u:=cdr u;
put('gamma,'simpfn, car u); u:=cdr u;
put('one,'simpfn, car u); u:=cdr u;
put('delta,'simpfn, car u); u:=cdr u;
put('intl,'simpfn, car u); u:=cdr u;
put('laplace,'simpfn, car u); u:=cdr u;
put('invlap,'simpfn, car u); u:=cdr u;
end;
%*******************************************************************
%* *
%* DIRECT LAPLACE TRANSFORMATION *
%* *
%*******************************************************************
put('laplace, 'simpfn, 'simplaplace);
lpsm!*:='( ((minus !=x))
(nil depends (reval (quote !=x)) lpvar!* )
(minus (times (one (minus !=x)) (sin !=x)) ) nil );
lpcm!*:='( (( minus !=x ))
(nil depends (reval (quote !=x)) lpvar!* )
(times (one (minus !=x)) (cos !=x)) nil );
lpshm!*:='( ((minus !=x))
(nil depends (reval (quote !=x)) lpvar!* )
(minus (times (one (minus !=x)) (sinh !=x)) ) nil );
lpchm!*:='( (( minus !=x ))
(nil depends (reval (quote !=x)) lpvar!* )
(times (one (minus !=x)) (cosh !=x)) nil );
lpse!*:= '( (!=x) (nil depends (reval(quote !=x)) lpvar!* )
(times (one !=x) (quotient (difference (expt e (times i !=x))
(expt e (minus (times i !=x))) )
(times 2 i) ) ) nil ) ;
lpce!*:= '( (!=x) (nil depends (reval(quote !=x)) lpvar!* )
(times (one !=x) (quotient (plus (expt e (times i !=x))
(expt e (minus (times i !=x))) )
2 ) ) nil ) ;
lpshe!*:= '( (!=x) (nil depends (reval(quote !=x)) lpvar!* )
(times (one !=x) (quotient (difference (expt e !=x)
(expt e (minus !=x)) )
2 ) ) nil );
lpche!*:= '( (!=x) (nil depends (reval(quote !=x)) lpvar!* )
(times (one !=x) (quotient (plus (expt e !=x)
(expt e (minus !=x)) )
2 ) ) nil );
lpexpt!*:= '( (e (plus !=x !=y)) (nil . t)
(times (expt e !=x) (expt e !=y) (one (plus !=x !=y)) ) nil );
symbolic procedure simplaplace u;
begin scalar e,r;
e:=lap!-save!-environment();
r:=errorset({'simplaplace!*,mkquote u},nil,nil);
lap!-restore!-environment(e);
if errorp r then typerr('laplace.u,"Laplace form")
else return laplace_fixup car r
end;
symbolic procedure laplace_fixup u;
% For some reason, results do not always come out in the most
% natural form. This is an attempt to fix this.
<<put('laplace,'simpfn,'simpiden);
u := simp aeval!* prepsq u;
put('laplace,'simpfn,'simplaplace);
u>> where varstack!* = nil;
symbolic procedure simplaplace!* u;
% Main procedure for Laplace transformation.
% U is in prefix form: (<expr> <lp.var> <il.var>), where
% <expr> is the object function,
% <lp.var> is the var. of the object function (intern. lp!&),
% <il.var> is the var. of the laplace transform(intern. il!&),
% and can be omitted - then il!& is assumed.
% Returns a standard quotient of Laplace transform.
begin scalar !*exp,!*mcd,v,w,transvar!*,!*precise;
% We need to make this run with precise on.
if null subfg!* then return mksq('laplace . u, 1);
if cddr u and null idp(w:=caddr u) or null idp(v:=cadr u)
then go to err;
v:= caaaar simp v;
transvar!* := w; % Needed for returning a Laplace form.
% Should the following be an error?
if null transvar!* then transvar!* := 'il!&;
if null idp v then go to err;
u:= car u ;
% Make environment for Laplace transform.
!*mcd := !*exp := t;
kord!*:= 'lp!& . 'il!& . kord!* ;
put('one,'simpfn,'lpsimp1);
put('gamma,'simpfn,'lpsimpg);
if !*ldone then put('expt,'opmtch,lpexpt!*.get('expt,'opmtch));
if !*lmon then
<< put('sin,'opmtch, lpse!* . get('sin,'opmtch));
put('cos,'opmtch, lpce!* . get('cos,'opmtch));
put('sinh,'opmtch, lpshe!* . get('sinh,'opmtch));
put('cosh,'opmtch, lpche!* . get('cosh,'opmtch)) >>
else
<< put('sin,'opmtch, lpsm!* . get('sin,'opmtch));
put('cos,'opmtch, lpcm!* . get('cos,'opmtch));
put('sinh,'opmtch, lpshm!* . get('sinh,'opmtch));
put('cosh,'opmtch, lpchm!* . get('cosh,'opmtch)) >>;
lpvar!*:= v; lpshift!*:=t;
if v neq 'lp!& then kord!*:=v . kord!*;
for each x in depl!* do if v memq cdr x then rplacd(x,'lp!& . cdr x);
% HM: resimplify u for rules before mcd goes off.
% ACH: However, it gives wrong results e.g. for laplace(sin(-x),x,p)
% rmsubs(); u := reval u;
off mcd;
u:= laplace1 list(u,v);
if w then u:=subf(numr u, list('il!& . w));
% Restore old env.
for each x in depl!* do
if 'lp!& memq cdr x then rplacd(x,delete('lp!&,cdr x));
put('one,'simpfn,'simpiden);
put('gamma,'simpfn,'simpiden);
kord!*:= cddr kord!*;
put('sin,'opmtch, cdr get('sin,'opmtch) );
put('cos,'opmtch, cdr get('cos,'opmtch) );
put('sinh,'opmtch, cdr get('sinh,'opmtch) );
put('cosh,'opmtch, cdr get('cosh,'opmtch) );
if !*ldone then put('expt,'opmtch,cdr get('expt,'opmtch) );
if erfg!* then erfg!*:=nil;
return u;
err: msgpri("Laplace operator incorrect",nil,nil,nil,t)
end where !*exp = !*exp, !*mcd = !*mcd;
put('sin,'lpfn,'(quotient k (plus (expt il!& 2) (expt k 2) )) );
put('cos,'lpfn,'(quotient il!& (plus (expt il!& 2) (expt k 2) )) );
put('sinh,'lpfn,'(quotient k (plus (expt il!& 2)
(minus (expt k 2)) )) );
put('cosh,'lpfn,'(quotient il!& (plus (expt il!& 2)
(minus (expt k 2)) )) );
put('one,'lpfn,'(quotient 1 il!&) );
put('expt,'lpfn,'(quotient (times (expt k d) (gamma (plus d 1)) )
(expt il!& (plus d 1)) ) );
put('delta,'lpfn, 1 );
symbolic procedure laplace1 u;
% Car u is in pref. form, cadr u is the var of the object function.
% Returns standard quotient of Laplace transform.
begin scalar v,w,z;
v := cadr u;
u := car u;
z:= simp!* u;
if denr z neq 1 then z := simp prepsq z; % *SQ must have occurred.
if denr z neq 1 then rederr list(u,"has non-trivial denominator");
z := numr z;
if v neq 'lp!& then << kord!*:=cdr kord!*;
z:=subla(list(v.'lp!&),z); z:=reorder z >>;
if erfg!* then return !*kk2q list
('laplace, subla(list('lp!& . lpvar!*), u), lpvar!*,transvar!*);
w:= nil ./ 1; u:=z; !*exp:=nil;
while u do
if domainp u
then << w:=addsq(w, lpdom u); u:=nil >>
else << w:=addsq(w, if (z:=lptermx lt u) then z
else !*kk2q list('laplace, subla
(list('lp!&.lpvar!*),prepsq !*t2q lt u),lpvar!*,transvar!*));
u:= red u >>;
return w;
end;
symbolic procedure lptermx u ;
% U is standard term, which may contain integer power of lp!&.
% Returns standard quot or nil, if Laplace transform is impossible.
begin scalar w ; integer n ;
if tvar u neq 'lp!& then return lpterm u
else if fixp cdar u
then if (n:=cdar u)>0 then nil else return lpunknown u
else return lpterm
( (list('expt,'lp!&,prepsq(cdar u ./ 1)) to 1) .* cdr u );
if (w:=lpform cdr u) then nil else return nil ;
a: % We use here the rule:
% laplace(x*fun(x),x)=-df(laplace(fun(x),x),il!&) ,or
% laplace(x**n*fun(x),x)=(-1)**n*df(laplace(fun(x),x),il!&,n);
if n=0 then return w;
w:=negsq diffsq(w,'il!&);
n:=n-1; go to a;
end;
symbolic procedure lpdom u ;
% We use here the rule: laplace(const,lp!&)=const/lp!&.
% U is domain. Returns standard quotient.
!*t2q (('il!& to -1) .* u) ;
symbolic procedure lpform u ;
% U is standard form, not containing integer powers of lp!&.
% Returns standard quot or nil, if Laplace transform is impossible.
begin scalar y,z ;
if domainp u
then return lpdom u
else if red u
then return
( if (y:=lpterm lt u) and (z:=lpform red u)
then addsq(y,z) else nil )
else return lpterm lt u ;
end ;
symbolic procedure lpterm u ;
% U is standard term, not containing integer powers of lp!&.
% Returns standard quot or nil, if Laplace transform is impossible.
begin scalar v,w,w1,y,z ;
v:=car u; % l.pow. - the first factor.
w:=cdr u; % l.coeff. - i.e. st.f.
if atom (y:=car v) or atom car y % I.e. atom or Lisp func.
then if not depends(y,'lp!&)
then return if (z:=lpform w)
then multpq(v,z) else nil
else if atom y then return lpunknown u
else if car y = 'expt
then return lpexpt(v,nil,w) else nil % Go next.
else return if not depends(prepsq(y./1),'lp!&)
then if (z:=lpform w)
then multpq(v,z) else nil
else lpunknown u;
% We can't handle v now, because nothing is known for w for now.
if domainp w then return lpfunc(v,w);
% If we have sum, and off exp.
if cdr w then return if (y:=lpterm list(v,car w)) and
(z:=lpterm(v . cdr w))then addsq(y,z) else nil;
w1:=cdar w; % l.coeff - i.e. st.f.
w :=caar w; % l.pow. - the second factor.
if not depends(if domainp(y:=car w) then y else prepsq(y./1),'lp!&)
then return if (z:=lpterm(v.w1)) then multpq(w,z) else nil
else if car y = 'expt then return lpexpt(w,v,w1);
% Now we have multiply of two functions.
if caar v = 'one and caar w = 'one
then return lpmult1(v,w,w1)
else return lpunknown u;
end ;
symbolic procedure lpunknown u ;
% Try to apply any previously given let rules for Laplace operator.
% U is standard term.
% Returns standard quotient or nil if matching not successful.
begin scalar d,z,w;
if domainp (d:=cdr u) and not !:onep d
then (u:= !*p2q car u) else (u:= !*t2q u);
u:= list('laplace, prepsq u, 'lp!&,transvar!*);
w:= list('laplace, cadr u,'lp!&); % HM: short rule form
if get('laplace,'opmtch) and
( (z:=opmtch u) or (z:=opmtch w))
then << !*exp:=t;
put('laplace,'simpfn,'laplace1);
z:=simp z; !*exp:=nil;
put('laplace,'simpfn,'simplaplace) >>;
if null z then return if !*lmsg
then msgpri("Laplace for", subla(list('lp!& . lpvar!*), cadr u),
if !*lmon or atom cadr u then "not known"
else "not known - try ON LMON",nil,nil)
else nil;
z:=subla(list('lp!&.lpvar!*), z);
return if domainp d and not !:onep d then multsq(z,d./1) else z;
end ;
symbolic procedure lpsimp1 u ;
% Simplify the one-function. % U is in prefix form.
% Returns standard quotient or nil ./ 1 if an error occurs.
begin scalar v,l,r ;
v:=subla(list(lpvar!* . 'lp!&),u);
if not depends(car v,'lp!&) then return 1 ./ 1;
v:= car simpcar v; % Standard form.
if mvar v neq 'lp!& then << !*mcd:=t; v:=subf(v,nil); !*mcd:=nil;
v:=multf(car v, recipf!* cdr v) >>;
if not(mvar v eq 'lp!& and !:onep ldeg v)
then go to err;
l:=lc v; r:=red v; % Standard form.
if null r then if minusf l then go to err else return 1 ./ 1;
v:=if minusf l then multsq(negf r ./ 1, 1 ./ negf l)
else multsq(r ./ 1, 1 ./ l);
if not minusf numr v then return 1 ./ 1;
if null lpshift!* then go to err
else return mksq(list('one,prepsq addsq(!*k2q 'lp!&, v)), 1);
err: if !*lmsg then msgpri("Laplace induces", 'one.u,
" which is not allowed", nil, 'hold);
return nil ./ 1;
end ;
symbolic procedure lpsimpg u ;
% Simplifies gamma(k), if k is rational and semiinteger.
% U is in prefix form. Returns standard quotient.
begin scalar n,v ;
u:= simpcar cdr u; % Gamma is now flagged "full".
if denr u neq 1
% Maybe we can do better than this.
then return mksq(list('gamma,prepsq u),1);
u := car u;
if domainp u and eqcar(u,'!:rn!:) and (cddr u = 2) % Semiint.
then return if (n:=cadr u) = 1
then mksq(list('sqrt,'pi),1)
else if n > 0 then
<< v:='!:rn!: . difference(n,2) . 2 ;
resimp !*t2q ( (list('gamma,rnprep!: v) to 1) .* v ) >>
else % N negative.
resimp !*t2q ( (list('gamma,rnprep!:('!:rn!:.plus(n,2) . 2)) to 1)
.* ('!:rn!:.(-2).(-n)) )
else return mksq(list('gamma,prepsq(u./1)),1);
end ;
symbolic procedure lpmult1 (u,v,w) ;
% Perform: one(l1*lp!&-r1)*one(l2*lp!&-r2) = one(l*lp!&-r),
% where l,r are those for the rightmost shifted one-function.
% U and v are standard powers for one-func., w is leading coeff.
% Returns standard quotient if all coeff. are domains, otherwise nil.
begin scalar u1,v1,l1,r1,l2,r2 ;
u1:= car simp cadar u;
v1:= car simp cadar v;
l1:=lc u1; l2:=lc v1;
r1:=red u1; r2:=red v1;
if domainp l1 and domainp l2 and domainp r1 and domainp r2
then if !:minusp adddm(multdm(r1,l2), !:minus multdm(r2,l1))
then return lpterm(u . w)
else return lpterm(v . w)
else return lpunknown list(u, v.w);
end ;
symbolic procedure lpexpt (u,v,w) ;
% Perform the rule: laplace(e**(l*lp!&)*fun(lp!&), lp!&) =
% sub(il!&=il!&-l, laplace(fun(lp!&),lp!&)),
% or call lpfunc for gamma-function.
% U is lpow for expt-func, v is other lpow or nil. W is lcoeff.
% Returns standard quotient or nil.
begin scalar p,q,r,z,l,la ;
r:=cdr u; % Degree for expt-func.
p:=cadar u; % First arg for expt.
q:=caddar u; % Second arg for expt.
if depends(p,'lp!&) then go to gamma;
!*exp:=t; q:=car simp q;
if mvar q neq 'lp!& then << !*mcd:=t; q:=subf(q,nil); !*mcd:=nil;
q:=multf(car q, recipf!* cdr q) >>;
if not !:onep r then q:=multf(q,r); !*exp:=nil;
if not(mvar q eq 'lp!& and !:onep ldeg q)
then return if null v then lpunknown(u . w)
else lpunknown list(u, v . w);
if (r:=red q) then
<< if !*ldone then << !*exp:=t;
w:=multf(w, car lpsimp1 list prepsq(q./1)); !*exp:=nil >>;
q:=list(lt q); r:=!*p2q(list('expt,p,prepsq(r./1)) to 1) >>;
if p neq 'e then q:=multf(q, !*kk2f list('log,p) );
z:= if null v then lpform w else lpterm(v.w);
if null z then return nil;
l:= prepsq !*f2q lc q;
la:=list('il!& . list('difference,'il!&,l) );
% Provide for those forms that contain the true transform variable.
if not(transvar!* eq 'il!&)
then z := subsq(z,list(transvar!* . 'il!&));
z:=subf(numr z,la);
return if r then multsq(r,z) else z;
gamma: % Check and call lpfunc for gamma-func.
return if null v
then if domainp w
then lpfunc(u,w)
else % if off exp
% if red w then if (z:=lpexpt(u,v,list(car w)) ) and
% (l:=lpexpt(u,v,cdr w)) then addsq(z,l) else nil else
if not depends((l:=mvar w),'lp!&)
then if (z:=lpexpt(u,nil,lc w))
then multpq(lpow w,z) else nil
else if not atom l and car l = 'expt
then lpexpt(lpow w,u,lc w)
else lpunknown(u . w)
else lpunknown list(u, v . w);
end ;
symbolic procedure lpfunc (u,v) ;
% Perform Laplace transform for intl-operator and simple functions:
% expt(arg,const), sin,cos,sinh,cosh,one,
% with args: k*lp!&-tau, where k>0, tau>=0 are const.
% U is standard power, v a domain element.
% Returns standard quotient or nil.
begin scalar ld,fn,w,var,ex,k,tau,c ;
ld:=cdr u; % Degree of func.
w:=car u; % Func in prefix form.
fn:=car w; % Name of func.
lintl: if fn neq 'intl then go to lexpt;
% Perform Laplace(intl(<expr>,<var>,0,lp!&), lp!&).
if not ( !:onep ld and cadddr w =0 and
car cddddr w = 'lp!& and idp(var:=caddr w) )
then return if !*lmsg then msgpri("Laplace integral",
subla(list('lp!& . lpvar!*), prepsq !*p2q u),
"not allowed", nil, nil) else nil;
ex:= subla(list(var . 'lp!&), cadr w);
lpshift!*:=nil; w:= laplace1 list(ex,'lp!&); lpshift!*:=t;
return if w then multsq(multd(v,!*p2f('il!& to -1))./1, w) else nil;
lexpt: if fn neq 'expt then go to lfunc;
% Perform Laplace(expt,(k*lp!&-tau),d), for d - not int. const.
ld:= multf(ld, car simp caddr w);
if minusf(addd(1,ld)) or depends(prepsq(ld./1), 'lp!&)
then return lpunknown(u.v);
ld:= prepsq !*f2q ld;
lfunc: % Perform Laplace transform for simple and one-function.
if fn = 'expt or (fn = 'one) or !:onep ld
then nil else return lpunknown(u.v);
!*exp:=t; ex:= car simp cadr w; !*exp:=nil;
if not( mvar ex = 'lp!& and !:onep ldeg ex )
then return lpunknown(u.v);
k:=lc ex; tau:=red ex;
if minusf k or (null lpshift!* and tau) then return
if !*lmsg then msgpri("Laplace for",
subla(list('lp!&.lpvar!*), w),"not allowed",nil,nil) else nil;
if tau and not minusf tau then return lpunknown(u.v);
c:= prepsq !*f2q k;
% Ind. lpfn gives Laplace transform for func(k*lp!&).
if (w:= get(fn,'lpfn))
then w:=car simp subla(list('k.c, 'd.ld), w);
return if null w
then lpunknown(u.v)
else if null tau
then multd(v, w) ./ 1
else multd(v, multf( w,!*kk2f list
('expt,'e,prepsq multsq(!*k2q 'il!&, quotsq(tau./1, k./1)) )
) ) ./ 1 ;
end ;
% Tables for Explicit Transforms for Delta Function. Note explicit
% construction for difference of arguments to reflect parser.
algebraic;
for all x,y,z let laplace(z*delta x,x,y) = sub(x=0,z);
for all k,x,y,z let laplace(z*delta(x+(-k)),x,y) = e**(y*-k)*sub(x=k,z);
for all x,y let laplace(df(delta x,x),x,y) = y;
for all n,x,y let laplace(df(delta x,x,n),x,y) = y**n;
for all k,x,y let laplace(df(delta(x+(-k)),x),x,y) = y*e**(-k*y);
for all k,n,x,y let laplace(df(delta(x+(-k)),x,n),x,y) = y**n*e**(-k*y);
symbolic;
%*******************************************************************
%* *
%* INVERSE LAPLACE TRANSFORMATION *
%* *
%*******************************************************************
put('invlap, 'simpfn, 'simpinvlap);
ile1!*:='( (e (times i !=x))
(nil depends(reval (quote !=x)) lpvar!*)
(plus (cos !=x) (times i (sin !=x))) nil );
ile2!*:='( (e (minus (times i !=x)))
(nil depends(reval (quote !=x)) lpvar!*)
(difference (cos !=x) (times i (sin !=x))) nil );
ile3!*:='( (e !=x )
(nil depends(reval (quote !=x)) lpvar!*)
(plus (cosh !=x) (sinh !=x)) nil );
ile4!*:='( (e (minus !=x))
(nil depends(reval (quote !=x)) lpvar!*)
(difference (cosh !=x) (sinh !=x)) nil );
ile5!*:='( (e (plus !=x !=y))
(nil and (not(depends(reval(quote !=x)) (quote i)))
(depends(reval(quote !=y)) (quote i)) )
(times (expt e !=x) (expt e !=y)) nil );
symbolic procedure simpinvlap u;
begin scalar r,e;
e:=lap!-save!-environment();
r:=errorset({'simpinvlap!*,mkquote u},nil,nil);
lap!-restore!-environment e;
if errorp r then typerr('invlap.u,"Laplace form")
else return invlap_fixup car r
end;
symbolic procedure invlap_fixup u;
% For some reason, results do not always come out in the most
% natural form. This is an attempt to fix this.
<<put('invlap,'simpfn,'simpiden);
u := simp aeval!* prepsq u;
put('invlap,'simpfn,'simpinvlap);
u>> where varstack!* = nil;
symbolic procedure simpinvlap!* u ;
% Main procedure for inverse Laplace transformation.
% U is in prefix form: (<expr> <il.var> <lp.var>) ,where
% <expr> is the laplace transform,
% <il.var> is the var. of the Laplace transform (intern. il!&),
% <lp.var> is the var. of the object function (intern. lp!&),
% and can be omitted - then lp!& is assumed.
% Returns a standard quotient of inverse Laplace transform.
begin scalar !*exp,!*mcd,v,w,!*precise;
% We need to make this run with precise on.
if null subfg!* then return mksq('invlap . u, 1);
if cddr u and null idp(w:=caddr u) then go to err;
v:= caaaar simp cadr u;
transvar!* := w;
if null idp v then go to err;
u:= car u ;
% Make environment for invlap transform.
!*exp := !*mcd := nil;
kord!*:= 'il!& . 'lp!& . kord!* ;
put('gamma,'simpfn,'lpsimpg);
put('one,'simpfn,'ilsimp1);
ilvar!*:=v; if v neq 'il!& then kord!*:=v.kord!*;
for each x in depl!* do if v memq cdr x then rplacd(x,'il!& . cdr x);
u:= invlap1 list(u,v);
put('invlap,'simpfn,'simpiden);
if w then << lpvar!*:=w; u:=subla(list('lp!& . w), u) >>
else lpvar!*:='lp!& ;
if !*ltrig or !*lhyp then << !*exp:=t;
if !*lhyp then put('expt,'opmtch,ile3!*.ile4!*.get('expt,'opmtch));
if !*ltrig then put('expt,'opmtch,ile1!*.ile2!*.get('expt,'opmtch));
put('expt,'opmtch, ile5!*.get('expt,'opmtch));
u:= simp prepsq u;
if !*ltrig and !*lhyp
then put('expt,'opmtch, cdr cddddr get('expt,'opmtch))
else put('expt,'opmtch, cdddr get('expt,'opmtch)) >>
else u:= resimp u;
% Restore old env.
for each x in depl!* do
if 'il!& memq cdr x then rplacd(x,delete('il!&,cdr x));
put('gamma,'simpfn,'simpiden);
put('one,'simpfn,'simpiden);
kord!*:= cddr kord!*;
return u;
err: msgpri("Invlap operator incorrect",nil,nil,nil,t);
end where !*exp = !*exp, !*mcd = !*mcd;
symbolic procedure invlap1 u;
% Car U is in prefix form, cadr u is the var of the Laplace transform.
% Returns standard quotient of inverse Laplace transform.
begin scalar v,w,z;
v := cadr u;
u := car u;
z:= simp!* u;
if denr z neq 1 then z := simp prepsq z; % *SQ must have occurred.
if denr z neq 1 then rederr list(u,"has non-trivial denominator");
z := numr z;
u := z;
if v neq 'il!& then << kord!*:=cdr kord!*;
u:=subla(list(v.'il!&),u); u:=reorder u >>;
w:= nil ./ 1;
while u do
if domainp u
then << w:=addsq(w, !*t2q((list('delta,'lp!&) to 1) .* u) );
u:= nil >>
else << w:=addsq(w, if (z:=ilterm (lt u,1,1,nil)) then z
else !*kk2q list('invlap, subla
(list('il!&.ilvar!*),prepsq !*t2q lt u), ilvar!*,transvar!*));
u:= red u >>;
return w;
end;
symbolic procedure ilterm (u, numf, denf, rootl) ;
% U is standard term, numf is standard form, with one term, and
% contains only powers from numerator of expression, depends on il!&,
% but not exponent. Denf is standard form, with one term, and
% contains only powers from denominator of expression, depends on il!&
% but not exponent. Rootl is assoc. list of: (<root> . <multiplity>).
% Returns standard quotient, or nil if inverse Laplace transform is
% impossible.
begin scalar v,v1,v2,w,y,z,p,p1 ;
v:=car u; w:=cdr u; v1:=car v; v2:=cdr v;
if not depends(if domainp v1 then v1 else prepsq(v1./1), 'il!&)
then return if (z:=ilform(w,numf,denf,rootl))
then multpq (v,z) else nil;
% V depends on il!&.
if atom v1
% the following clause "if n1 neq il& then" introduced by HM
then (if not(v1 = 'il!&) then return ilunknown(u,numf,denf))
else if atom car v1 % I.e. Lisp func.
then return
if car v1 = 'expt
then ilexpt(v,nil,w,numf,denf,rootl)
else if domainp w
then ilexptfn(v,w,numf,denf)
else if cdr w
then if(y:=ilterm(list(v,lt w),numf,denf,rootl))
and (z:=ilterm(v.cdr w,numf,denf,rootl))
then addsq(y,z)
else nil
else ilterm(list(lpow w,v.(lc w)),numf,denf,rootl);
% May be infinite recursion above, if mult. of two unknown func.
% Mvar is atom 'il!& or standard form, since exp off.
if numberp v2 and fixp v2
then
if v2 > 0
then
if atom v1
then return ilform(w, multf(!*p2f v,numf), denf, rootl)
else nil
else return ilroot(v, w, numf, denf, rootl)
else return
ilexpt(list('expt,
if domainp v1 then v1 else prepsq(v1./1),
prepsq(v2./1)) to 1, nil, w, numf,
denf, rootl);
% Now v1 remains as a standard form and v2>0.
v:= if !:onep v2 then v1 else !*p2f v;
if red v1 then
<< !*exp:=t; y:=numr subf(v,nil); z:=y;
while z do if domainp z then z:=nil
else if ldeg z < 0 then if depends
(if domainp(p1:=mvar z) then p1 else prepsq(p1 ./1), 'il!&)
then << p:=t; z:=nil >> else z:=addf(lc z, red z)
else z:=addf(lc z,red z);
if p then w:=multf(y, w) else numf:=multf(v,numf);
!*exp:=nil >> else numf:=multf(v,numf);
return ilform(w,numf,denf,rootl);
end;
symbolic procedure ilform (u, numf, denf, rootl) ;
% U is a standard form. Numf, denf, rootl are the same as in ILTERM.
% Returns standard quotient or nil if invlap is impossible.
begin scalar y,z ;
return if domainp u
then if (z:=ilresid(numf,denf,rootl))
then multsq(u ./ 1, z) else nil
else if null red u
then ilterm(lt u,numf,denf,rootl)
else if (y:=ilterm(lt u,numf,denf,rootl)) and
(z:=ilform(red u,numf,denf,rootl))
then addsq(y,z) else nil;
end ;
symbolic procedure ilunknown (u, numf, denf) ;
% We try here to apply any previously given let rules for Laplace
% operator. U is standard term, numf, denf are the same.
% Returns standard quotient or nil if matching not successful.
begin scalar d,z,w;
if domainp (d:=cdr u) then if !:onep d
then u:=!*t2q u else u:=!*p2q car u
else u:=!*t2q u;
if numf neq 1 then u:=multsq(u, numf./1);
if denf neq 1 then u:=multsq(u,1 ./denf);
u:= list('invlap, prepsq u,'il!&,transvar!*);
% HM: alternative shorter form for rule match
w:= list('invlap, cadr u, 'il!&);
if get('invlap,'opmtch) and
((z:=opmtch u) or (z:=opmtch w))
then << !*exp:=t;
put('invlap,'simpfn,'invlap1);
z:=simp z; !*exp:=nil;
put('invlap,'simpfn,'simpinvlap) >>;
if null z and !*lmsg then msgpri("Invlap for",
subla(list('il!& . ilvar!*), cadr u),
"not known", nil, nil);
return if null z then nil
else if domainp d and not !:onep d
then multsq(z, d ./ 1) else z;
end ;
symbolic procedure ilsimp1 u ;
% Simplify the one-function. U is in prefix form.
% Returns standard quotient.
if atom car u then 1 ./ 1 else mksq('one . u, 1);
symbolic procedure ilexpt (u, v, w, numf, denf, rootl) ;
% Perform the rule: invlap(e**(-l*il!&)*fun(il!&), il!&) =
% sub(lp!&=lp!&-l, invlap(fun(il!&),il!&)), for l > 0,
% or call ilfunc for gamma-function.
% U is lpow for expt-function, v is other lpow or nil,
% W is lcoeff (standard form), numf, denf, rootl are the same.
% Returns standard quotient or nil.
begin scalar p,q,r,z,l ;
r:=cdr u; % Degree for expt-func.
p:=cadar u; % First arg for expt.
q:=caddar u; % Second arg for expt.
if depends(p,'il!&)then go to gamma;
!*exp:=t; q:=car simp q;
if mvar q neq 'il!& then << !*mcd:=t; q:=subf(q,nil); !*mcd:=nil;
q:=multf(car q, recipf!* cdr q) >>;
if not !:onep r then q:=multf(q,r); !*exp:=nil;
if not((mvar q = 'il!&) and !:onep ldeg q and minusf lc q)
then return if null v then ilunknown(u.w,numf,denf)
else ilunknown(list(u,v.w),numf,denf);
if (r:=red q) then<< q:=list(lt q);
r:=!*p2q(list('expt,p,prepsq(r./1)) to 1) >>;
if p neq 'e then q:=multf(q, !*kk2f list('log,p) );
z:= if null v then ilform(w,numf,denf,rootl)
else ilterm(v.w,numf,denf,rootl);
if null z then return nil;
l:= list('plus, 'lp!&, prepsq((lc q)./1));
z:= subf(numr z, list('lp!& . l) ) ; % Standard quotient.
% If you want shifted one-func. to remain always in obj. func.
if !*lione then z:=multsq(z, !*kk2q list('one,l) );
return if r then multsq(r,z) else z ;
gamma: % Check and call ilfunc if gamma-func. case.
return if null v
then if domainp w
then ilexptfn(u,w,numf,denf)
else if red w
then if (z:=ilexpt(u,nil,list(car w),numf,denf,rootl)) and
(l:=ilexpt(u,nil,cdr w,numf,denf,rootl))
then addsq(z,l) else nil
else if not depends(if domainp(l:=mvar w) then l
else prepsq(l./1), 'il!&)
then if (z:=ilexpt(u,nil,lc w,numf,denf,rootl))
then multpq(lpow w,z) else nil
else if not atom l and (car l = 'expt)
then ilexpt(lpow w,u,lc w,numf,denf,rootl)
else if atom l or not atom car l
then ilterm(list(lpow w,u.(lc w)),numf,denf,rootl)
else ilunknown(u.w,numf,denf)
else ilunknown(list(u,v.w),numf,denf) ;
end ;
symbolic procedure ilexptfn (u, v, numf, denf) ;
% Perform invlap for expt function - i.e., gamma-function case.
% U is standard power for expt, v is domain, numf, denf the same.
% Returns standard quotient or nil.
begin scalar ex,dg,fn,k,a,b,y,d ;
ex:=car u; dg:=cdr u; fn:=car ex;
if fn neq 'expt then go to unk;
d:=caddr ex; if atom(ex:=cadr ex) then k:=t;
!*exp:=t; ex:=car simp ex;
dg:=multd(dg,car simp d); a:=lc ex;
if not(domainp a and !:onep a) then
<< ex:=multf(ex, recipf!* a);
a:=!*kk2f list('expt,prepsq(a./1),prepsq(dg./1)) >>;
b:=red ex; !*exp:=nil;
if (mvar ex neq 'il!&) or (ldeg ex neq 1) or
depends(prepsq(b./1),'il!&) then go to unk;
if (numf=1) and (denf=1) then go to ret;
% We must have identical monomials in numf, denf and in expt-func.
y:= multf(multf(numf, !*kk2f list('expt,
prepsq(ex./1),prepsq(dg./1)) ), recipf!* denf);
if cdr y or (lc lc y neq 1) or (car mvar lc y neq 'expt)
or (not k and (mvar y neq ex))
or (k and (mvar y neq mvar ex)) then go to unk;
dg:=addd(ldeg y,dg);
ret: if minusf dg then d:=prepsq(negf dg ./ 1) else go to unk;
if (y:=get(fn,'ilfn))
then y:=car simp subla(list('d.d), y) else go to unk;
if b then y:=multd(v, multf(y, !*kk2f list
('expt,'e,prepsq(multf(!*k2f 'lp!&,negf b) ./1)) ))
else y:=multd(v, y);
return if domainp a and !:onep a then y./1 else multf(a,y)./1;
unk: return ilunknown(u.v, numf, denf);
end ;
put('expt,'ilfn,'(quotient (expt lp!& (plus d (minus 1))) (gamma d)));
symbolic procedure addrootl (root,mltpl,rootl) ;
% Add roots with multiplity at head of rootl - an assoc. list.
begin scalar parr ;
parr:=assoc(root,rootl);
if parr then << mltpl:= mltpl + cdr parr;
rootl:= delete(parr,rootl) >>;
return (root . mltpl) . rootl ;
end ;
symbolic procedure recipf!* u ;
% U is standard form. Returns st.f. for u to (-1), by off mcd.
begin scalar d;
if domainp u then if !:onep u then return 1
else if !:onep negf u then return -1
else if fieldp u then nil
else if (d:=get(dmode!*,'i2d))
then u:=apply1(d,u) else u:=mkratnum u
else return if cdr u then !*p2f(u to (-1))
else multf(!*p2f(mvar u to (-ldeg u)), recipf!* lc u);
return dcombine(1,u,'quotient);
end ;
symbolic procedure ilroot (u,v,numf,denf,rootl);
% Find the roots of polynomial of first and second degree.
% U is standard power - the polynomial, v is the remaining st.f.
% Numf, denf, rootl are the same. Returns standard quot or nil.
begin scalar dg,ex,a,b,c,z,x1,x2 ;
dg:=-cdr u; ex:=car u; % dg>0;
if atom ex then return ilform(v,numf,
multf(!*p2f('il!& to dg),denf), addrootl(nil,dg,rootl) );
if atom car ex then return ilunknown(u.v,numf,denf);
!*exp:=t; ex:=subf(ex,nil); !*exp:=nil;
if not depends(prepsq ex, 'il!&) then return
if (z:=ilform(v,numf,denf,rootl)) then multpq(u,z) else nil;
ex:=car ex;
if ldeg ex > 2 then return il3pol(u,v,numf,denf,rootl);
a:=lc ex;
if depends(prepsq(a./1),'il!&)
then return ilunknown(u.v,numf,denf);
if not(domainp a and !:onep a) then
<< !*exp:=t; a:=recipf!* a; ex:=multf(ex,a);
if dg>1 then a:=exptf(a,dg); !*exp:=nil >>;
if ldeg ex = 2 then go to lbin;
lmon: if (b:=red ex)
then << rootl:=addrootl(negf b, dg, rootl);
denf:= if !:onep dg then multf(ex, denf)
else multpf(ex to dg, denf) >>
else << rootl:=addrootl(nil, dg, rootl);
denf:= multpf('il!& to dg, denf) >>;
go to ret;
lbin: if (b:=red ex)
then if domainp b then << c:=b; b:=nil >>
else if mvar b = 'il!& then << c:=red b; b:=lc b >>
else << c:=b; b:=nil >>
else c:=nil ;
if depends(prepsq(b./1),'il!&) or depends(prepsq(c./1),'il!&)
then return ilunknown(u.v,numf,denf);
if null b and null c
then << rootl:=addrootl(nil, 2*dg, rootl);
denf:=multpf('il!& to (2*dg), denf) >>
else << !*exp:=t; b:=multd('!:rn!: . ((-1) . 2), b);
c := simp list('sqrt,prepsq(addf(multf(b,b),negf c)./1));
if fixp denr c
then c := multd(('!:rn!: . 1 . denr c),numr c)
else rederr {"invalid laplace denominator",denr c};
x1:=addf(b,c); x2:=addf(b,negf c); !*exp:=nil;
if x1 = x2 then << rootl:=addrootl(x1,2*dg,rootl);
x1:=(('il!& to 1).*1) .+ negf x1;
denf:=multpf(x1 to (2*dg),denf) >>
else << rootl:=addrootl(x2,dg,addrootl(x1,dg,rootl));
x1:=(('il!& to 1).*1) .+ negf x1;
x2:=(('il!& to 1).*1) .+ negf x2;
if not !:onep dg then
<< x1:=!*p2f(x1 to dg); x2:=!*p2f(x2 to dg) >>;
denf:=multf(x2,multf(x1,denf)) >> >>;
ret: z:=ilform(v,numf,denf,rootl);
return if (domainp a and !:onep a) then z
else if null z then nil else multsq(a./1, z);
end;
symbolic procedure il3pol (u, v, numf, denf, rootl) ;
% Find the roots of polynomial of third and higher degree.
% U is standard power - the polynomial, v is the remaining st.f.
% Numf, denf, rootl are the same. Returns standard quot or nil.
(begin scalar a,d,p,y,z,w;
if !*rounded then go to unk;
d:=-cdr u; p:=car u;
!*exp:=t; !*mcd:=t;
% We must now convert rationals, if any, to standard quotients.
% Since MCD was previously off, we must use limitedfactors here,
% since the regular factorization turns EZGCD on.
!*limitedfactors := t;
y:=p; p:=nil./1;
while y do if domainp y then << p:=addsq(p,!*d2q y); y:=nil >>
else << a:=1; z:=list car y; % S.F. with 1 term only.
while not domainp z do
<< w:=lpow z;
% distinguish between mvar=kernel/form
w:=if kernlp car w then !*p2f w else
exptf(car w,cdr w);
a:=multf(a,w);
z:=lc z
>>;
p:=addsq(p,multsq(a./1,!*d2q z)); y:=red y >>;
if ((a:=cdr p) neq 1) and (d neq 1) then a:=exptf(a,d);
z := fctrf car p;
!*exp:=nil; !*mcd:=nil;
% if length z = 2 then go to unk; % No factors.
% corrected (HM):
if length z=2 and cdr cadr z=1 then go to unk;
if car z neq 1 then errach list(car z,"found in IL3POL");
z:=cdr z; y:=v;
while z do << p:= caar z;
if cdar z neq 1 then p := exptf(p,cdar z);
if d neq 1 then p:=exptf(p,d);
y:=multf(y,recipf!* p); z:=cdr z >>;
y:=ilform(y,numf,denf,rootl);
if null y then go to unk
else return if a = 1 then y else multsq(a./1, y);
unk: return ilunknown(u.v,numf,denf);
end) where !*limitedfactors := !*limitedfactors;
symbolic procedure ilresid (numf, denf, rootl) ;
% Apply the residue theorem at last.
% Numf, denf, rootl are standard forms. Returns standard quot or nil.
begin scalar n,d,ndeg,ddeg,m,x,y,z,w ;
!*exp:=t; n:=numr subf(numf,nil); !*exp:=nil;
z:=nil ./ 1; w:=nil ./ 1; x:=n; % Result accumulated in w.
while x and not domainp x do
<< y:=car x; x:=cdr x;
if depends(prepsq(cdr y./1),'il!&) or (caar y neq 'il!&
and depends(caar y,'il!&) )
then if (z:=ilterm(y,1,denf,rootl))
then << w:=addsq(w,z); n:=delete(y,n) >>
else x:=nil >> ;
if null z then return ;
% Now n is polynomial of il!& with constant coeff.
ndeg:=if not domainp n and mvar n = 'il!& then ldeg n else 0;
!*exp:=t; d:=numr subf(denf,nil); !*exp:=nil;
ddeg:=if not domainp d and mvar d = 'il!& then ldeg d else 0;
if ndeg < ddeg then go to resid;
!*exp:=t; y:=qremf(n,d); !*exp:=nil;
n:=cdr y; x:=car y; % N is remainder polynomial.
while x do if domainp x
then << w:=addsq(w, !*t2q(('(delta lp!&) to 1).* x)); x:=nil >>
else if mvar x neq 'il!&
then << w:=addsq(w, multsq(!*kk2q '(delta lp!&), !*t2q lt x) );
x:=red x >>
else << w:=addsq(w, multsq(!*kk2q list('df,list('delta,'lp!&),
'lp!&,ldeg x), lc x ./ 1) ); x:=red x >> ;
resid: if null rootl then return w ;
x:=caar rootl; m:=cdar rootl;
if null x then y:=!*p2f('il!& to m)
else << y:=(('il!& to 1) .* 1) .+ negf x;
if m neq 1 then y:=!*p2f(y to m) >>;
!*exp:=t; y:=numr subf(y,nil);
y:=car qremf(d,y); !*exp:=nil; % D is quotient - remainder = 0.
z:=multpf('(expt e (times il!& lp!&)) to 1, n); % Numerator.
y:=recipf!* y;
z:=multf(z,y) ./ 1;
while (m:=m-1) > 0 do z:=diffsq(z, 'il!&);
x:= if null x then 0 else prepsq(x./1); % Root in prefix form.
!*exp:=t; z:=subf(numr z, list('il!&.x)); % One residue as st.q.
if not depends(prepsq z, 'lp!&)
then z:=multsq(z, !*kk2q '(one lp!&));
if (m:=cdar rootl) > 2 then while (m:=m-1) > 1 do
z:=multf(car z,'!:rn!: . 1 . m)./1;
w:=addsq(w,z); !*exp:=nil;
rootl:=cdr rootl; go to resid;
end;
endmodule;
end;