Fix AW/101

	* algebra/laplace.spad.pamphlet (lapkernel): Handle derivatives.

1 files changed, 24 insertions, 1 deletions

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)