[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[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
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}
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Axiom-developer] 20070921.01.tpd.patch,
daly <=