axiom-developer
[Top][All Lists]

## [Axiom-developer] 20070921.01.tpd.patch

 From: daly Subject: [Axiom-developer] 20070921.01.tpd.patch Date: Fri, 21 Sep 2007 11:00:27 -0500

The call

laplace(log(z),z,w)

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:

-log(w)-gamma
-------------
w

======================================================================
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
index 6f81b2e..76cc837 100644
@@ -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
++ using t as the new variable or "failed" if unable to find
++ a closed form.
+      ++ Handles only rational \spad{f(s)}.

@@ -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}