From e193c16d6eec5d174a24987a8590247c4c4227d1 Mon Sep 17 00:00:00 2001 From: dos-reis Date: Thu, 16 Oct 2008 04:41:31 +0000 Subject: Fix AW/101 * algebra/laplace.spad.pamphlet (lapkernel): Handle derivatives. --- src/ChangeLog | 5 +++++ src/algebra/laplace.spad.pamphlet | 25 ++++++++++++++++++++++++- src/testsuite/interpreter/aw-101.input | 2 ++ 3 files changed, 31 insertions(+), 1 deletion(-) create mode 100644 src/testsuite/interpreter/aw-101.input diff --git a/src/ChangeLog b/src/ChangeLog index e3cbf6c4..ea652e11 100644 --- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,8 @@ +2008-10-15 Waldek Hebisch + + Fix AW/101 + * algebra/laplace.spad.pamphlet (lapkernel): Handle derivatives. + 2008-10-15 Gabriel Dos Reis * interp/compiler.boot (backendCompile1): Move to c-util.boot. diff --git a/src/algebra/laplace.spad.pamphlet b/src/algebra/laplace.spad.pamphlet index b3a6767c..74921607 100644 --- a/src/algebra/laplace.spad.pamphlet +++ b/src/algebra/laplace.spad.pamphlet @@ -161,32 +161,55 @@ 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 + if F has SpecialFunctionCategory then + 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) diff --git a/src/testsuite/interpreter/aw-101.input b/src/testsuite/interpreter/aw-101.input new file mode 100644 index 00000000..294be6d4 --- /dev/null +++ b/src/testsuite/interpreter/aw-101.input @@ -0,0 +1,2 @@ +-- Used to crash the system. +laplace(log z, z, w) -- cgit v1.2.3