LocalSolutionsOrigin.ml

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

(* The service LocalSolutionsOrigin computes all transseries solutions *)
(* of a differential equation at 0. A transseries solution consists of *)
(* a Puiseux series of log-power type, multiplied by some exponential  *)
(* factor and a power factor x^alpha with alpha an algebraic number.   *)

(* Here is an overview how several services call each other to compute *)
(* the local solutions of a differential equation at some point:       *)
(* LocalSolutionsPoint: transseries solutions at an arbitrary point    *)
(* -> LocalSolutionsOrigin: transseries solutions at the origin        *)
(*    -> LogSeriesAlpha: log-power series with x^alpha factor          *)
(*       -> LogSeriesZero: log-power series solutions                  *)

INCLUDE "preamble.ml"

let title (_, _, _, _, _) = <:text<Local Solutions at the Origin>>


(* Compute the generalized exponents. This proc basically calls   *)
(* DEtools[gen_exp] and postprocesses its output a little bit,    *)
(* by simplifying RootOf expressions, extracting the ramification *)
(* indices, and flattening the lists. The final output is of the  *)
(* form [[e1, r1], [e2, r2], ...] where ej is a generalized       *)
(* exponent (i.e., exp(int(ej/x,x)) is an exponential factor) and *)
(* rj is the corresponding ramification index (an integer).       *)
let generalized_exponents eqn1 y1 x1 =
  <<
    (proc(eqn, y, x, $)
      local genexp, i, g, a, r, s, T, pf, ram;

      genexp := DEtools[gen_exp](eqn, y(x), T);

      # get rid of RootOf-s that have no index parameter
      genexp := map(allvalues, genexp);

      for i from 1 to nops(genexp) do
        g := genexp[i];

        # we expect the last element to be of the form a*T^r=x.
        if not(match(op(1, g[-1]) = a*T^r, T, 's')
          and evalb(op(0, g[-1]) = `=`)
          and evalb(op(2, g[-1]) = x))
        then error("Generalized exponent different than expected.")
        end if;

        # ram is the ramification index and pf the above prefactor a.
        ram := subs(s, r);
        pf := subs(s, a);
        g := subsop(-1 = NULL, g);
        g := map(a -> subs(T = x/a, g), [allvalues(RootOf(_Z^ram - pf))]);
        g := map(a -> [a, ram], ListTools:-Flatten(g));
        genexp[i] := g;
      end do;
      return ListTools:-FlattenOnce(genexp);
    end proc)($(eqn1), $(y1), $(x1))
  >>


(* This function is called from the main service in order to loop over  *)
(* all exponential factors. It returns a pair (doc, sols) that consists *)
(* of the documented trace of the calculations and an Ocaml list        *)
(* containing the solutions corresponding to this exponential factor.   *)
let deal_with_exp_factor parallel eqn y x u n exp_factor ram =

  (* Perform the transformation that corresponds to the ramification index. *)
  let trans_eqn =
    << collect(
      evala(gfun:-algebraicsubs($(eqn), $(y) - $(x)^$(ram), $(y)($(x)))),
      {diff, $(y)}, factor) >> in
  let ram_text = if <:bool< $(ram) = 1 >> then <:par<>> else
    <:par<
      The $(t_ent: Glossary.g "ramification index") for this exponential
      factor is <:isymb<$(ram)>>. Therefore, removing the exponential
      factor as it is would result in a differential equation with higher
      order and produce spurious solutions. In order to avoid this
      the substitution
      <:imath< $(symb: x) \mapsto $(symb: x)^$(symb: ram)>>
      is performed, which leads to the differential equation
      <:dmath< $(symb: trans_eqn) = 0. >>
    >> in

  (* Remove the exponential factor from the differential equation. *)
  let trans_eqn =
    << collect(
      evala(gfun:-`diffeq*diffeq`(
        $(trans_eqn),
        gfun:-holexprtodiffeq(1/$(exp_factor), $(y)($(x))), $(y)($(x)))),
      {diff, $(y)}, factor) >> in
  let exp_text = if <:bool< $(exp_factor) = 1 >> then <:par<>> else
    <:par<
      The exponential factor <:isymb< $(exp_factor) >> is removed
      by applying the transformation
      <:imath< << $(y)($(x)) >> \mapsto
        << simplify(1/$(exp_factor)) * $(y)($(x)) >> >>.
      This results in the differential equation
      <:dmath< $(symb: trans_eqn) = 0. >>
    >> in

  (* Compute all solutions of the form x^alpha * logseries. *)
  let series_text, lsers =
    if parallel then
      DC.inline_service (LogSeriesAlpha.descr (trans_eqn, y, x, u, n) None),
      LogSeriesAlpha.obj ((trans_eqn, y, x, u, n), ())
    else
      LogSeriesAlpha.doc_obj trans_eqn y x u n
  in
  let lsers = CommonTools.symb_list_of_symb lsers in

  (* Include the exponential factor into the log-power-series solutions, *)
  (* and undo the ramification substitution. *)
  List.iter
    (fun a ->
       <:unit<
         $(a):-exp_factor := subs($(x) = $(x)^(1/$(ram)), $(exp_factor));
         $(a):-alpha := $(a):-alpha / $(ram);
         $(a):-finite_part := subs($(x) = $(x)^(1/$(ram)), $(a):-finite_part);
         $(a):-ramification := $(ram);
       >>)
    lsers;

  (ram_text @:@ exp_text) @@@ series_text, lsers


(* Compute all transseries solutions of a differential equation at 0. *)
(* The functions u[i](n) are used to encode the coefficient sequences *)
(* in the series expansions. *)
let doc_obj eqn y x u n =

  let genexp = generalized_exponents eqn y x in

  (* Determine equivalence classes (same exponential) among the *)
  (* generalized exponents. *)
  let genexp =
    <<
      DDMF:-equivalence_classes(
        $(genexp),
        (a, b) -> not(has(simplify(a[1] - b[1]), $(x))) and evalb(a[2] = b[2]))
    >> in

  (* Now take some representative for each equivalence class. *)
  (* The sort is only to enforce deterministic behaviour of maple. *)
  let genexp = << sort(map(a -> op(1, a), $(genexp))) >> in

  (* Delete the constant part of each polynomial; this corresponds to a *)
  (* power of x only and should go into the alpha of LogSeriesAlpha. *)
  (* Of course, one could use this information directly, but we decided to *)
  (* perform the extraction of alpha ourselves using the indicial equation. *)
  let genexp =
    << map(a -> [a[1] - subs($(x) = infinity, a[1]), a[2]], $(genexp)) >> in

  (* Transform the Laurent polynomials into exponentials. *)
  let exp_factors =
    CommonTools.symb_list_of_symb
      << map(a -> exp(int(a[2] * a[1] / $(x), $(x))), $(genexp)) >> in

  (* Extract the ramification indices. *)
  let rams = CommonTools.symb_list_of_symb << map(a -> a[2], $(genexp)) >> in

  (* For each exponential factor compute the log-powerseries solutions. *)
  let parallel = List.length exp_factors > 1 in
  let fcall = deal_with_exp_factor parallel eqn y x u n in
  let contents, sol_list = List.split (List.map2 fcall exp_factors rams) in

  (* The exponential factors as they appear in the result. *)
  let exp_factors =
    List.map (fun a -> << $(List.hd a):-exp_factor >>) sol_list in
  let nexps = List.length exp_factors in

  (* If there are nontrivial exponential factors, write something about *)
  (* how they were obtained. *)
  let intro =
    if nexps = 1 && <:bool< $(List.hd exp_factors) = 1 >>
    then DC.ent_null
    else
      <:par<
        The modular algorithm described in $(t_ent:Bibliography.b "ClHo04")
        is employed to determine that the local solutions involve the
        exponential factor$(str: Wording.ending_of_int nexps)
      >> @:@ (Wording.enumeration_of_maples exp_factors) @:@ <:par<.>> in

  (* If there are more than just one exponential factor, create subsections. *)
  let content =
    if nexps = 1
    then List.hd contents
    else List.fold_left (@@@) DC.ent_null
      (List.map2
         (fun a b ->
           DC.section
             <:text<Solutions with exponential factor <:isymb< $(a) >>>> b)
         exp_factors
         contents) in

  intro @@@ content, CommonTools.symb_of_symb_list (List.flatten sol_list)


let_service LocalSolutionsOrigin
  (eqn : diffeq maple)
  (y : name maple)
  (x : name maple)
  (u : name maple)
  (n : name maple) :
  DC.sec_entities * any maple with { title = title } =
  doc_obj eqn y x u n

Generated by GNU Enscript 1.6.5.90.