LogSeriesAlpha.ml

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

(* The service LogSeriesAlpha computes all log-power series solutions  *)
(* of a differential equation at 0 including possibly an x^alpha       *)
(* factor, i.e., solutions of the form        *)
(*   x^alpha * \sum_{n=0}^\infty \sum_{j=0}^{t-1} u_j(n) \ln(x)^j x^n, *)
(* where alpha is an algebraic number.                                 *)
(* It is assumed that at least one such solution exists, i.e., that 0  *)
(* is not an irregular singular point of the differential equation.    *)

(* 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<Series Solutions with Power Factor at the Origin>>


(* This function is called from the main service in order to loop over  *)
(* all shifts that are determined by the indicial equation. It returns  *)
(* a pair (doc, sols) that consists of the documented trace of the      *)
(* calculations and an Ocaml list sols containing the solutions that    *)
(* correspond to the given shift.                                       *)
let deal_with_shift parallel eqn y x u n shift =

  (* Remove the power factor. *)
  let trans_eqn =
    << collect(
      evala(gfun:-`diffeq*diffeq`(
        $(eqn),
        $(x) * diff($(y)($(x)), $(x)) + $(shift) * $(y)($(x)), $(y)($(x)))),
      {diff, $(y)}, factor) >> in
  let fac_text = if <:bool< $(shift) = 0 >> then <:par<>> else
    <:par<
      The power factor <:isymb< $(x)^$(shift) >> is removed by means of the
      transformation
      <:imath< << $(y)($(x)) >> \mapsto << $(x)^$(shift) * $(y)($(x)) >> >>
      which results in the differential equation
      <:dmath< $(symb: trans_eqn) = 0. >>
    >> in

  (* Compute the log-power-series solutions. *)
  let content, psers =
    if parallel then
      DC.inline_service (LogSeriesZero.descr (trans_eqn, y, x, u, n) None),
      LogSeriesZero.obj ((trans_eqn, y, x, u, n), ())
    else
      LogSeriesZero.doc_obj trans_eqn y x u n
  in
  let psers = CommonTools.symb_list_of_symb psers in

  (* Include the power factor into the log-power-series solutions. *)
  List.iter (fun a -> <:unit< $(a):-alpha := $(shift) + $(a):-alpha >>) psers;

  fac_text @@@ content, psers


(* Compute all local solutions of a differential equation at 0, which *)
(* are of the type                                                    *)
(*   x^\alpha \sum_{n=0}^\infty \sum_{j=0}^{t-1} u_j(n) \ln(x)^j x^n  *)
(* where alpha is an algebraic number.                                *)
(* The names u and n have to be given as arguments.                   *)
let doc_obj eqn y x u n =

  let intro =
    <:par<
      Studying the $(t_ent: Glossary.g "indicial equation")
      of the differential equation above explains what exponents
      of <:isymb< $(x) >> can be involved in local series solutions.
    >> in

  (* Compute the roots of the indicial equation and group them into *)
  (* classes determined by integer shift equivalence.               *)
  (* The equivalence classes are sorted by their first elements,    *)
  (* particularly for enforcing deterministic behaviour of maple.   *)
  let ind_eqn = << evala(gfun:-indicialeq(
    gfun:-formatdiffeq([$(eqn), $(y)($(x))]), $(x), $(n))) >> in
  let ind_roots = << [solve($(ind_eqn) = 0, $(n))] >> in
  let ecs =
    << DDMF:-equivalence_classes(
      $(ind_roots), (a, b) -> type(a - b, integer)) >> in
  let ecs = << map(sort, $(ecs)) >> in
  let ecs = CommonTools.symb_list_of_symb
    << sort($(ecs), (a, b) -> [a[1], b[1]] = sort([a[1], b[1]])) >> in
  let necs = List.length ecs in

  (* Create the corresponding content. *)
  let ind_eqn_text =
    if necs = 1 && <:bool< $(List.hd ecs)[1] = 0 >>
    then
      (* The indicial equation is always displayed by the LogSeriesZero. *)
      (* In this case the same equation would be displayed twice in a row. *)
      <:par<>>
    else
      <:par<
        The indicial equation in this case is
        <:dmath< $(symb: ind_eqn) = 0 >> and its
      >> @:@
      (if <:bool< nops($(ind_roots)) = 1 >>
      then <:par<only root is >>
      else <:par<roots are >>)
      @:@
      (Wording.enumeration_of_maples (CommonTools.symb_list_of_symb ind_roots))
      @:@ <:par<.>> in
  let ecs_text = if necs <= 1 then <:par<>> else
    <:par<
      These roots are grouped into classes determined by integer shift
      equivalence, i.e., all elements in an equivalence class differ by
      an integer only. Thus the equivalence classes are
    >> @:@ (Wording.enumeration_of_maples ecs) @:@ <:par<,
      each of which is treated separately in one of the following subsections.
    >> in

  (* Extract the "smallest" element s from each equivalence class, so that *)
  (* all others are of the form s+n with n being a positive integer. *)
  let shifts =
    List.map (fun a -> << $(a)[1] + min(map(b -> b - $(a)[1], $(a))) >>) ecs in

  let parallel = List.length shifts > 1 in
  let fcall = deal_with_shift parallel eqn y x u n in
  let contents, sols = List.split (List.map fcall shifts) in

  (* Create subsections if not all indicial roots are are *)
  (* in the same equivalence class. *)
  let content =
    if necs = 1
    then List.hd contents
    else List.fold_left (@@@) DC.ent_null
      (List.map2
         (fun a b ->
           DC.section
             <:text<Solutions with power factor <:isymb< $(x)^$(a) >>>> b)
         shifts
         contents) in

  (intro @:@ ind_eqn_text @:@ ecs_text) @@@ content,
  CommonTools.symb_of_symb_list (List.flatten sols)


let_service LogSeriesAlpha
  (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.