diff --git a/changelog b/changelog index 48493c1..d085c48 100644 --- a/changelog +++ b/changelog @@ -1,3 +1,6 @@ +20070920 tpd src/input/Makefile add bug101.input regression test +20070920 tpd src/input/bug101.input test laplace(log(z),z,w) bug 101 +20070920 wxh src/algebra/laplace.spad fix laplace(log(z),z,w) bug 101 20070916 tpd src/input/Makefile add bug103.input regression test 20070916 tpd src/input/bug103.input test solve(z=z,z) bug fix 20070916 tpd src/algebra/polycat.spad solve(z=z,z) bug fix diff --git a/src/algebra/laplace.spad.pamphlet b/src/algebra/laplace.spad.pamphlet index 6f81b2e..76cc837 100644 --- a/src/algebra/laplace.spad.pamphlet +++ b/src/algebra/laplace.spad.pamphlet @@ -161,41 +161,68 @@ LaplaceTransform(R, F): Exports == Implementation where lapkernel(f, t, tt, ss) == (k := retractIfCan(f)@Union(K, "failed")) case "failed" => "failed" - empty?(arg := argument(k::K)) or not empty? rest arg => "failed" + empty?(arg := argument(k::K)) => "failed" + is?(op := operator k, "%diff"::SE) => + not( #arg = 3) => "failed" + not(is?(arg.3, t)) => "failed" + fint := eval(arg.1, arg.2, tt) + s := name operator (kernels(ss).1) + ss * locallaplace(fint, t, tt, s, ss) - eval(fint, tt = 0) + not (empty?(rest arg)) => "failed" member?(t, variables(a := first(arg) / tt)) => "failed" is?(op := operator k, "Si"::SE) => atan(a / ss) / ss is?(op, "Ci"::SE) => log((ss**2 + a**2) / a**2) / (2 * ss) is?(op, "Ei"::SE) => log((ss + a) / a) / ss + -- digamma (or Gamma) needs SpecialFunctionCategory + -- which we do not have here + -- is?(op, "log"::SE) => (digamma(1) - log(a) - log(ss)) / ss "failed" + -- Below we try to apply one of the texbook rules for computing + -- Laplace transforms, either reducing problem to simpler cases + -- or using one of known base cases locallaplace(f, t, tt, s, ss) == zero? f => 0 -- one? f => inv ss (f = 1) => inv ss + + -- laplace(f(t)/t,t,s) + -- = integrate(laplace(f(t),t,v), v = s..%plusInfinity) (x := tdenom(f, tt)) case F => g := locallaplace(x::F, t, tt, vv := new()$SE, vvv := vv::F) (x := intlaplace(f, ss, g, vv, vvv)) case F => x::F oplap(f, tt, ss) + + -- Use linearity (u := mkPlus f) case List(F) => +/[locallaplace(g, t, tt, s, ss) for g in u::List(F)] (rec := splitConstant(f, t)).const ^= 1 => rec.const * locallaplace(rec.nconst, t, tt, s, ss) + + -- laplace(t^n*f(t),t,s) = (-1)^n*D(laplace(f(t),t,s), s, n)) (v := atn(f, t)) case Record(coef:F, deg:PI) => vv := v::Record(coef:F, deg:PI) is?(la := locallaplace(vv.coef, t, tt, s, ss), oplap) => oplap(f,tt,ss) (-1$Integer)**(vv.deg) * differentiate(la, s, vv.deg) + + -- Complex shift rule (w := aexp(f, t)) case Record(coef:F, coef1:F, coef0:F) => ww := w::Record(coef:F, coef1:F, coef0:F) exp(ww.coef0) * locallaplace(ww.coef,t,tt,s,ss - ww.coef1) + + -- Try base cases (x := lapkernel(f, t, tt, ss)) case F => x::F - -- last chance option: try to use the fact that - -- laplace(f(t),t,s) = s laplace(g(t),t,s) - g(0) where dg/dt = f(t) - elem?(int := lfintegrate(f, t)) and (rint := retractIfCan int) case F => - fint := rint :: F - -- to avoid infinite loops, we don't call laplace recursively - -- if the integral has no new logs and f is an algebraic function - empty?(logpart int) and algebraic?(f, t) => oplap(fint, tt, ss) - ss * locallaplace(fint, t, tt, s, ss) - eval(fint, tt = 0) + +-- -- The following does not seem to help computing transforms, but +-- -- quite frequently leads to loops, so I (wh) disabled it for now +-- -- last chance option: try to use the fact that +-- -- laplace(f(t),t,s) = s laplace(g(t),t,s) - g(0) where dg/dt = f(t) +-- elem?(int := lfintegrate(f, t)) and (rint := retractIfCan int) case F => +-- fint := rint :: F +-- -- to avoid infinite loops, we don't call laplace recursively +-- -- if the integral has no new logs and f is an algebraic function +-- empty?(logpart int) and algebraic?(f, t) => oplap(fint, tt, ss) +-- ss * locallaplace(fint, t, tt, s, ss) - eval(fint, tt = 0) oplap(f, tt, ss) setProperty(oplap,SPECIALDIFF,dvlap@((List F,SE)->F) pretend None) @@ -228,6 +255,7 @@ InverseLaplaceTransform(R, F): Exports == Implementation where ++ Laplace transform of \spad{f(s)} ++ using t as the new variable or "failed" if unable to find ++ a closed form. + ++ Handles only rational \spad{f(s)}. Implementation ==> add @@ -246,8 +274,12 @@ InverseLaplaceTransform(R, F): Exports == Implementation where ilt(expr,var,t) == expr = 0 => 0 r := univariate(expr,kernel(var)) + + -- Check that r is a rational function such that degree of + -- the numarator is lower that degree of denominator not(numer(r) quo denom(r) = 0) => "failed" not( freeOf?(numer r,var) and freeOf?(denom r,var)) => "failed" + ilt1(r,t::F) hintpac := TranscendentalHermiteIntegration(F, UP) diff --git a/src/input/Makefile.pamphlet b/src/input/Makefile.pamphlet index b286aa6..b2d5fe9 100644 --- a/src/input/Makefile.pamphlet +++ b/src/input/Makefile.pamphlet @@ -287,7 +287,8 @@ REGRES= algaggr.regress algbrbf.regress algfacob.regress alist.regress \ arrows.regress assign.regress atansqrt.regress \ asec.regress bags.regress bbtree.regress \ binary.regress bop.regress bstree.regress bouquet.regress \ - bug100.regress bug103.regress bug10069.regress \ + bug100.regress bug101.regress \ + bug103.regress bug10069.regress \ bugs.regress bug10312.regress bug6357.regress bug9057.regress \ calcprob.regress \ calculus2.regress calculus.regress cardinal.regress card.regress \ @@ -502,7 +503,8 @@ FILES= ${OUT}/algaggr.input ${OUT}/algbrbf.input ${OUT}/algfacob.input \ ${OUT}/bags.input ${OUT}/bbtree.input ${OUT}/bern.input \ ${OUT}/bernpoly.input ${OUT}/binary.input ${OUT}/bop.input \ ${OUT}/bouquet.input ${OUT}/bstree.input ${OUT}/bug6357.input \ - ${OUT}/bug9057.input ${OUT}/bug100.input ${OUT}/bug103.input \ + ${OUT}/bug9057.input ${OUT}/bug100.input ${OUT}/bug101.input \ + ${OUT}/bug103.input \ ${OUT}/bug10069.input ${OUT}/bug10312.input \ ${OUT}/calcprob.input ${OUT}/calculus.input \ ${OUT}/cardinal.input ${OUT}/card.input ${OUT}/carten.input \ @@ -678,7 +680,8 @@ DOCFILES= \ ${DOC}/bernpoly.input.dvi ${DOC}/binary.input.dvi \ ${DOC}/bop.input.dvi ${DOC}/bouquet.input.dvi \ ${DOC}/bstree.input.dvi ${DOC}/bug10069.input.dvi \ - ${DOC}/bug100.input.dvi ${DOC}/bug103.input.dvi \ + ${DOC}/bug100.input.dvi ${DOC}/bug101.input.dvi \ + ${DOC}/bug103.input.dvi \ ${DOC}/bug10312.input.dvi ${DOC}/bug6357.input.dvi \ ${DOC}/bug9057.input.dvi ${DOC}/bugs.input.dvi \ ${DOC}/c02aff.input.dvi ${DOC}/c02agf.input.dvi \ diff --git a/src/input/bug101.input.pamphlet b/src/input/bug101.input.pamphlet new file mode 100644 index 0000000..b82caf0 --- /dev/null +++ b/src/input/bug101.input.pamphlet @@ -0,0 +1,59 @@ +\documentclass{article} +\usepackage{axiom} +\begin{document} +\title{\$SPAD/src/input bug101.input} +\author{Timothy Daly} +\maketitle +\begin{abstract} +\end{abstract} +\eject +\tableofcontents +\eject +@ +The call +\begin{verbatim} +laplace(log(z),z,w) +\end{verbatim} +causes an infinite loop between lfintegrate and monomialIntegrate. + +Waldek modified the code to fail and return unevaluated. + +A more reasonable return value, not supported by Axiom, would be: +\begin{verbatim} + -log(w)-gamma + ------------- + w +\end{verbatim} + +<<*>>= +)spool bug101.output +)set message test on +)set message auto off +)clear all + +--S 1 of 2 +laplace(log(z),z,w) +--R +--R (1) laplace(log(z),z,w) +--R Type: Expression Integer +--E 1 + +--S 2 of 2 +laplace(log(z),w,z) +--R +--R log(z) +--R (2) ------ +--R z +--R Type: Expression Integer +--E 2 +@ +<<*>>= +)spool +)lisp (bye) + +@ +\eject +\begin{thebibliography}{99} +\bibitem{1} nothing +\end{thebibliography} +\end{document}