axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / LAPLACE.spad

)abbrev package LAPLACE LaplaceTransform
++ Author: Manuel Bronstein
++ Date Created: 30 May 1990
++ Date Last Updated: 13 December 1995
++ Description: 
++ This package computes the forward Laplace Transform.

LaplaceTransform(R, F) : SIG == CODE where
  R : Join(EuclideanDomain, OrderedSet, CharacteristicZero,
           RetractableTo Integer, LinearlyExplicitRingOver Integer)
  F : Join(TranscendentalFunctionCategory, PrimitiveFunctionCategory,
           AlgebraicallyClosedFunctionSpace R)

  SE  ==> Symbol
  PI  ==> PositiveInteger
  N   ==> NonNegativeInteger
  K   ==> Kernel F
  OFE ==> OrderedCompletion F
  EQ  ==> Equation OFE
  ALGOP       ==> "%alg"
  SPECIALDIFF ==> "%specialDiff"

  SIG ==> with

    laplace : (F, SE, SE) -> F
      ++ laplace(f, t, s) returns the Laplace transform of \spad{f(t)}
      ++ using \spad{s} as the new variable.
      ++ This is \spad{integral(exp(-s*t)*f(t), t = 0..%plusInfinity)}.
      ++ Returns the formal object \spad{laplace(f, t, s)} if it cannot
      ++ compute the transform.

  CODE ==> add

    import IntegrationTools(R, F)
    import ElementaryIntegration(R, F)
    import PatternMatchIntegration(R, F)
    import PowerSeriesLimitPackage(R, F)
    import FunctionSpaceIntegration(R, F)
    import TrigonometricManipulations(R, F)

    locallaplace : (F, SE, F, SE, F) -> F
    lapkernel    : (F, SE, F, F) -> Union(F, "failed")
    intlaplace   : (F, F, F, SE, F) -> Union(F, "failed")
    isLinear     : (F, SE) -> Union(Record(const:F, nconst:F), "failed")
    mkPlus       : F -> Union(List F, "failed")
    dvlap        : (List F, SE) -> F
    tdenom       : (F, F) -> Union(F, "failed")
    atn          : (F, SE) -> Union(Record(coef:F, deg:PI), "failed")
    aexp         : (F, SE) -> Union(Record(coef:F, coef1:F, coef0:F), "failed")
    algebraic?   : (F, SE) -> Boolean

    oplap := operator("laplace"::Symbol, 3)$BasicOperator

    laplace(f,t,s) == locallaplace(complexElementary(f,t),t,t::F,s,s::F)

    -- returns true if the highest kernel of f is algebraic over something
    algebraic?(f, t) ==
      l := varselect(kernels f, t)
      m:N := reduce(max, [height k for k in l], 0)$List(N)
      for k in l repeat
         height k = m and has?(operator k, ALGOP) => return true
      false

    -- differentiate a kernel of the form  laplace(l.1,l.2,l.3) w.r.t x.
    -- note that x is not necessarily l.3
    -- if x = l.3, then there is no use recomputing the laplace transform,
    -- it will remain formal anyways
    dvlap(l, x) ==
      l1 := first l
      l2 := second l
      x = (v := retract(l3 := third l)@SE) => - oplap(l2 * l1, l2, l3)
      e := exp(- l3 * l2)
      locallaplace(differentiate(e * l1, x) / e, retract(l2)@SE, l2, v, l3)

    -- returns [b, c] iff f = c * t + b
    -- and b and c do not involve t
    isLinear(f, t) ==
      ff := univariate(f, kernel(t)@K)
      ((d := retractIfCan(denom ff)@Union(F, "failed")) case "failed")
        or (degree(numer ff) > 1) => "failed"
      freeOf?(b := coefficient(numer ff, 0) / d, t) and
        freeOf?(c := coefficient(numer ff, 1) / d, t) => [b, c]
      "failed"

    -- returns [a, n] iff f = a * t**n
    atn(f, t) ==
      if ((v := isExpt f) case Record(var:K, exponent:Integer)) then
        w := v::Record(var:K, exponent:Integer)
        (w.exponent > 0) and
          ((vv := symbolIfCan(w.var)) case SE) and (vv::SE = t) =>
            return [1, w.exponent::PI]
      (u := isTimes f) case List(F) =>
        c:F  := 1
        d:N  := 0
        for g in u::List(F) repeat
          if (rec := atn(g, t)) case Record(coef:F, deg:PI) then
            r := rec::Record(coef:F, deg:PI)
            c := c * r.coef
            d := d + r.deg
          else c := c * g
        zero? d => "failed"
        [c, d::PI]
      "failed"

    -- returns [a, c, b] iff f = a * exp(c * t + b)
    -- and b and c do not involve t
    aexp(f, t) ==
      is?(f, "exp"::SE) =>
        (v := isLinear(first argument(retract(f)@K),t)) case "failed" =>
           "failed"
        [1, v.nconst, v.const]
      (u := isTimes f) case List(F) =>
        c:F := 1
        c1 := c0 := 0$F
        for g in u::List(F) repeat
          if (r := aexp(g,t)) case Record(coef:F,coef1:F,coef0:F) then
            rec := r::Record(coef:F, coef1:F, coef0:F)
            c   := c * rec.coef
            c0  := c0 + rec.coef0
            c1  := c1 + rec.coef1
          else c := c * g
        zero? c0 and zero? c1 => "failed"
        [c, c1, c0]
      if (v := isPower f) case Record(val:F, exponent:Integer) then
        w := v::Record(val:F, exponent:Integer)
        (w.exponent ^= 1) and
          ((r := aexp(w.val, t)) case Record(coef:F,coef1:F,coef0:F)) =>
            rec := r::Record(coef:F, coef1:F, coef0:F)
            return [rec.coef ** w.exponent, w.exponent * rec.coef1,
                                            w.exponent * rec.coef0]
      "failed"

    mkPlus f ==
      (u := isPlus numer f) case "failed" => "failed"
      d := denom f
      [p / d for p in u::List(SparseMultivariatePolynomial(R, K))]

    -- returns g if f = g/t
    tdenom(f, t) ==
      (denom f exquo numer t) case "failed" => "failed"
      t * f

    intlaplace(f, ss, g, v, vv) ==
      is?(g, oplap) or ((i := integrate(g, v)) case List(F)) => "failed"
      (u:=limit(i::F,equation(vv::OFE,plusInfinity()$OFE)$EQ)) case OFE =>
        (l := limit(i::F, equation(vv::OFE, ss::OFE)$EQ)) case OFE =>
          retractIfCan(u::OFE - l::OFE)@Union(F, "failed")
        "failed"
      "failed"

    lapkernel(f, t, tt, ss) ==
      (k := retractIfCan(f)@Union(K, "failed")) case "failed" => "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
      (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
      oplap(f, tt, ss)

    setProperty(oplap,SPECIALDIFF,dvlap@((List F,SE)->F) pretend None)