Artifact 3b41771db86115549c29cd8ec56e1b8d3c238d9a35f0e4d7e9070e9fddb2739b:
- File
r35/lib/laplace.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 41758) [annotate] [blame] [check-ins using] [more...]
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 laplace)); fluid '(!*exp !*limitedfactors !*mcd !*rounded depl!* kord!* transvar!*); 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 car r; end; 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 v,w,svfg,transvar!*; 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. svfg:=!*exp . !*mcd ; % !*mcd:=nil; HM: too early here (pattern match fails) !*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 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; go to ret >>; put('laplace,'simpfn,'simpiden); !*exp:=car svfg; svfg:=cdr svfg; if !*exp then u:=resimp u; !*mcd:=svfg; if !*mcd then u:=resimp u; put('laplace,'simpfn,'simplaplace); % ret: return reorder numr u ./ reorder denr u ; ret: return u ; err: msgpri("Laplace operator incorrect",nil,nil,nil,t); end ; 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:= car simpcar 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 car r; end; 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 v,w,svfg ; 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. svfg:=!*exp . !*mcd ; !*mcd:=nil; !*exp:=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)); !*exp:=car svfg; svfg:=cdr svfg; if !*exp and not(!*ltrig or !*lhyp) then u:=resimp u; !*mcd:=svfg; if !*mcd then u:=resimp u; put('invlap,'simpfn,'simpinvlap); put('one,'simpfn,'simpiden); put('gamma,'simpfn,'simpiden); kord!*:= cddr kord!*; % return reorder numr u ./ reorder denr u ; return u; err: msgpri("Invlap operator incorrect",nil,nil,nil,t); end; 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:=numr simp list('sqrt,prepsq(addf(multf(b,b),negf c)./1)); 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;