LaplaceTransform.ml

(* Copyright INRIA and Microsoft Corporation, 2008-2013. *)
(* DDMF is distributed under CeCILL-B license. *)

INCLUDE "preamble.ml"

(* Exception private to this module. *)
exception Quit

let title _ = <:text<Laplace Transform>>

let truncation marshaled_sf y point order =
  (* The last two arguments (notation, parametrized) are only relevant *)
  (* for the descr, but here we need the obj. *)
  let linear_comb =
    LocalExpansion.obj ((marshaled_sf, point, y, "", true), ()) in
  if <:bool< $(linear_comb) = FAIL >> then raise Quit;
  let trunc = TruncatedSeries.obj ((linear_comb, ""), order) in
  << DDMF:-dynamow_to_maple($(trunc)) >>


(* BS May 2010. Preliminary version by Shiv Shankar rewritten from scratch. *)
let_service LaplaceTransform
  (marshaled_sf : string)
  (y : name maple)
  (def : diffeq maple) :
  DC.sec_entities * bool with { title = title } =

  (* The calculation may fail at many places, jumping to the end by raising *)
  (* the exception Quit. *)
  try

    let sf = DB.unmarshal_sf marshaled_sf in
    let x = SF.var_of_t sf in

    <:unit<
      deq := gfun:-formatdiffeq([$(def), $(y)($(x))]);
      sing := {solve(deq[-1], $(x))};
    >> ;

    (* Check that no singularity lies on (0,infinity). *)
    (* The Laplace transform might still exist, but we do not treat this yet. *)
    if <:bool< select(type,sing,positive) <> {} >> then raise Quit ;

    (* This expansion order is not justified, but seems ok for all functions *)
    (* currently in the dictionary. *)
    let sufficient_order = 4 in
    let tr point = truncation marshaled_sf y point sufficient_order in
    let expansion_at_0, expansion_at_infinity = tr << 0 >>, tr << infinity >> in

    (* Check that the behaviour at 0 is ok. *)
    (* Deal only with absolutely convergent integrals in this version. *)
    if <:bool< member(0,sing) and
         evalb(limit($(x) * $(expansion_at_0), $(x)=0, right) <> 0) >>
    then raise Quit ;

    (* Compute half-plane of convergence and check that it is nonempty. *)
    <:unit<
      bpoint := limit(log( eval($(expansion_at_infinity),
                                [sin=1,cos=1]) ) / $(x),
                      $(x)=infinity)
    >> ;
    if <:bool< bpoint = infinity >> then raise Quit ;

    (* Compute differential equation satisfied by the Laplace transform. *)
    <:unit<
      try
        deq_laplace := gfun:-Laplace($(def),$(y)($(x)),'diffeq') ;
      catch "no valid initial conditions":
      end try
    >> ;
    if <:bool< not assigned(deq_laplace) >> then raise Quit ;

    (* Solve it. *)
    <:unit<
      if has(deq_laplace, diff) then # this is a true differential equation
        if type(deq_laplace,set) then
          deq_laplace:=op(remove(type,deq_laplace,`=`))
        fi;
        laplacetransf:=subs(dsolve(deq_laplace,$(y)($(x))),$(y)($(x)));
        inds:=indets(laplacetransf,name);
        inds:=inds intersect {(seq(_C || i,i=1..nops(inds)))};
        dim:=nops(inds);
        behaviour_at_0:=convert($(expansion_at_0),polynom);
        if has(behaviour_at_0,ln($(x))) then # Do not try to solve (yet).
          laplacetransf:=FAIL
        else
          ser1:=MultiSeries:-multiseries(
            laplacetransf -
            add(coeff(behaviour_at_0,$(x),i)*i!*$(x)^i,i=0..nops(inds)) +
              O($(x)^(nops(inds)+1)),$(x),nops(inds)+2);
          # extract all coefficients
          ser1:={eval(ser1,SERIES=proc(u,v) op(v) end)};
          # do not remove simplify: it is necessary for arccot(x)
          sys:=map(simplify,ser1); # set of coefficients
          # Functions of x in the coefficients, cancel their coefficients.
          if has(sys,$(x)) then
            sys:=expand(sys); # otherwise coeffs will fail
            sys:=map(coeffs,sys,select(has,indets(sys),$(x)) minus {_var[$(x)]})
          fi;
          laplacetransf:=subs(solve(sys,inds),laplacetransf);
        fi
      else
        laplacetransf:=solve(deq_laplace,$(y)($(x)));
      fi
    >> ;
    (* FIXME: Reimplement this hack so as not to depend on the implementation *)
    (* of DynaMoW. *)
    (* Cut here to get out of try...catch by dynamow, *)
    (* so that the tests below are performed in all cases. *)
    <:unit<
      if has(laplacetransf,{_C,_C1}) then laplacetransf:=FAIL;
      else laplacetransf := normal( 1/s*eval(laplacetransf,$(x)=1/s) ) ;
      fi;
      if laplacetransf = FAIL and assigned(deq_laplace) then
        if has(deq_laplace,_C) then deq_laplace:=FAIL;
        else
          deq_laplace :=
            gfun:-algebraicsubs(deq_laplace,$(x)*$(y)-1,$(y)($(x)));
          deq_laplace := subs($(x)=s,
            collect(numer(eval(deq_laplace, $(y)($(x))=$(y)($(x))/$(x))),
                                          [diff,$(y)($(x))],factor))
        fi
      fi
    >> ;

    if <:bool< laplacetransf = FAIL and deq_laplace = FAIL >> then
      raise Quit ;

    let body =
      (* Domain of convergence. *)
      (if <:bool< bpoint<>-infinity >> then
        let bp = << bpoint >> in
        <:par<For <:imath< \Re(s) > $(symb:bp) >>, >>
      else <:par<For any complex number <:imath< s >>, >>) @:@

      (* Result. *)
      if <:bool< laplacetransf <> FAIL >> then
        let var = SF.var_of_t sf in
        let rep = << subs($(var) = t, $(sf.SF.rep)) >> in
        <:par<<:dmath<
          \int_0^\infty e^{-ts} $(symb:rep) \,dt = << laplacetransf >> .
        >>>>
      else
        <:par<
          the $(t_ent:Glossary.g "Laplace transform") satisfies the following
          linear differential equation: <:dmath< << deq_laplace >> = 0 . >>
        >> in
    (DC.section (title ()) body, true)

  with Quit -> (DC.section (title ()) DC.ent_null, false)

Generated by GNU Enscript 1.6.5.90.