LogSeriesZero.ml

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

(* The service LogSeriesZero computes all log-power series solutions   *)
(* of a differential equation at 0, i.e., solutions of the form        *)
(*   \sum_{n=0}^\infty \sum_{j=0}^{t-1} u_j(n) \ln(x)^j x^n.           *)
(* For this purpose it is assumed that at least one such solution      *)
(* exists, i.e., that the origin is not an irregular singular point,   *)
(* and that the equation has been transformed such that the smallest   *)
(* integer root of the indicial equation is 0.                         *)

(* 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<Log-Power Series at the Origin>>


(* This function is called from the main service in order to loop over  *)
(* all canonical solutions. Each canonical solution is characterized by *)
(* case = << [k, d] >> where k is the index of the root                 *)
(* n_0 = int_roots[k] and d is an index labeling the different          *)
(* solutions that correspond to this root (i.e., whose series expansion *)
(* starts at n_0).                                                      *)
(* It returns a pair (doc, sol) that consists of the documented trace   *)
(* of the calculations and the solution as a Maple Record.              *)
let deal_with_case x u n bigQ q int_roots mults case =
  let k = << $(case)[1] >>
  and d = << $(case)[2] >> in
  let n0 = << $(int_roots)[$(k)] >>
  and m = << $(mults)[$(k)] >> in

  (* Compute t early for a mathematically sound local expansion. *)
  (* (The origin of t is explained later.) *)
  let t =
    << $(d) + 1 + add($(mults)[i], i = $(k) + 1 .. nops($(int_roots))) >> in
  let solution_form =
    <:par<
      <:imath< << ln($(x))^$(d) * $(x)^$(n0) >> +
        O(<< ln($(x))^($(t) - 1) * $(x)^($(n0) + 1) >>) >>\
    >> in

  let intro =
    <:par<Since the indicial equation has the integer root <:isymb< $(n0) >> >>
    @:@
    (if <:bool< $(m) = 1 >>
    then <:par<, there exists a $(t_ent:Glossary.g "formal logarithmic sum") >>
    else
      <:par<
        with multiplicity <:isymb< $(m) >>, there exist <:isymb< $(m) >>
        linearly independent $(t_ent:Glossary.g
                                 ~display:(Some "formal logarithmic sums")
                                 "formal logarithmic sum")
      >>)
    @:@
    <:par<
      whose lowest exponent of <:isymb< $(x) >> is exactly <:isymb< $(n0) >>.
    >> @:@
    (if <:bool< $(m) = 1 >>
    then <:par<>>
    else
      <:par<This subsection deals with the >> @:@
      (Wording.ordinal <:int< $(d) + 1 >>) @:@
      <:par< of these solutions, which takes the form >> @:@
      solution_form @:@ <:par<.>>) in

  (* Compute t and describe what's going on. *)
  let m1 =
    CommonTools.symb_list_of_symb
      << [seq($(mults)[i], i = $(k)+1 .. nops($(mults)))] >> in
  let sum_mults = List.map (fun a -> "+" ^ (Maple.serialization_of_t a)) m1 in
  let sum = List.fold_left (^) (Maple.serialization_of_t d) sum_mults in
  let t_text = if <:bool< $(t) = 1 >> then <:par<>> else
    <:par<
      There, the highest possible exponent of <:isymb< ln($(x)) >>
      is determined by the sum of the current index
      (which is <:isymb< $(d) >>)
      and the multiplicities of all integer roots of the indicial
      equation that are greater than the current root, <:isymb< $(n0) >>;
      this gives <:imath< $(str: sum) = << $(t) - 1 >> >>.
    >> in

  (* This is basically formula (3.9), *)
  (* where u[k](n) has been substituted by k!*u[k](n-n0) *)
  let s = << $(n) = $(n) + $(n0) >> in
  (* Since the Q's contain the parameters of the given special function, *)
  (* use of local variables is essential here! *)
  let recsys = <<
    proc(bigQ, u, t, s, n) local i, j, l;
      [seq(
        factor(subs(s, bigQ[1,1])) * u[i](n) =
          add(factor(subs(s, -bigQ[1,j-i+1]) * binomial(j, i)) *
            u[j](n), j = i+1 .. t-1) +
          add(add(factor(subs(s, -bigQ[l+1,j-i+1]) * binomial(j, i)) *
            u[j](n-l), j = i .. t-1), l = 1 .. $(q)),
        i = 0 .. t-1)]
    end proc($(bigQ), $(u), $(t), $(s), $(n)) >> in

  (* Explain and display the recurrences. *)
  let recs1 = CommonTools.symb_list_of_symb recsys in
  let recs2 =
    List.map
      (fun a -> <:par<<:dmath< $(symb: a), >>>>)
      (List.rev (List.tl recs1)) in
  let recs3 = recs2 @ [ <:par<<:dmath< << $(List.hd recs1) >>.>>>> ] in
  let recs = List.fold_left (@:@) <:par<>> recs3 in
  let smnd = << add($(u)[i]($(n))*ln($(x))^i, i = 0 .. $(t)-1) * $(x)^$(n) >> in
  let ansatz =
    if <:bool< $(n0) = 0 >>
    then <:par<<:dmath< \sum_{$(symb: n)=0}^\infty $(symb: smnd) >>>>
    else
      <:par<<:dmath< << $(x)^$(n0) >>
        \sum_{$(symb: n)=0}^\infty $(symb: smnd) >>>> in
  let recsys_text =
    <:par<Applying the differential equation to the \
      $(t_ent:Glossary.g "ansatz")>>
    @:@ ansatz @:@
    <:par<and comparing equal powers of <:isymb< $(x) >> >>
    @:@
    (if <:bool< nops($(recsys)) = 1 >>
     then <:par<leads to the recurrence >>
     else <:par<and <:isymb< ln($(x)) >> leads to the system of recurrences >>)
    @:@ recs @:@
    <:par<An alternative expression for >> @:@
    (if <:bool< nops($(recsys)) = 1 >>
     then <:par<this recurrence>>
     else <:par<these recurrences>>) @:@
    <:par< \
      is in terms of the polynomials <:imath< Q_j($(symb: n)) >> computed above
      and reads
      <:dmath< \sum_{\ell=0}^$(symb: q) \sum_{j=i}^<< $(t)-1 >>
        {j\choose i} Q_{\ell}^{(j-i)}(<< $(n) + $(n0) >>) \,
        $(symb: u)_j($(symb: n)-\ell) = 0 \quad (0 \leq i < $(symb: t)).>>
    >> in

  (* Initial values *)
  let inits = <<
    [
      seq($(u)[j](0) = 0, j = 0 .. $(d)-1),
      $(u)[$(d)](0) = 1,
      seq($(u)[j](0) = 0, j = $(d)+1 .. $(t)-1)
   ] >> in
  let inits_text1 =
    <:par<Since the solution should be of the form >> @:@
    solution_form @:@
    <:par<\
      , the initial value$(str: Wording.ending_of_int <:int< $(t) >>) for
      <:imath< $(symb: n) = 0 >> must be chosen as follows:
    >> @:@
    (Wording.enumeration_of_maplelist inits) @:@
    <:par<. >> in
  let arb_values =
    << [seq(seq($(u)[j]($(int_roots)[i] - $(n0)),
      j = 0 .. $(mults)[i]-1),
      i = $(k)+1 .. nops($(int_roots))) ] >> in
  let narb = <:int< nops($(arb_values)) >> in
  let inits_text2 = if narb = 0 then <:par<>> else
    <:par<A closer inspection of the recurrences above reveals that >> @:@
    (Wording.enumeration_of_maplelist arb_values) @:@
    <:par< never occur>> @:@
    (if narb = 1 then <:par<s. Its value>> else <:par<. Their values>>) @:@
    <:par< therefore may be chosen arbitrarily, <:imath<0>> for example. >> in
  let inits = << [op($(inits)), op(map(a -> a = 0, $(arb_values)))] >> in

  (* Compute the remaining initial values using the recurrences and *)
  (* the assumption that the u[i]'s at negative arguments are zero. *)
  let eqns = << [seq(op(subs($(n) = i, $(recsys))),
    i = 1 .. max($(q)-1, max($(int_roots)) - $(n0)))] >> in
  let eqns = << subs($(inits), $(eqns)) >> in
  let eqns = << subs(
    [seq(seq($(u)[j](i) = 0, j = 0 .. $(t)-1), i = 1-$(q) .. -1)],
    $(eqns)) >> in
  let vars =
    << indets($(eqns), specfunc(integer, {seq($(u)[j], j = 0 .. $(t)-1)})) >> in
  let sols = << solve($(eqns), [op($(vars))])[1] >> in
  let is =
    << {seq(i, i = 0 .. $(q)-1), op(map(a -> a - $(n0), $(int_roots)))} >> in
  let sols = << select(a -> member(op(1, op(1, a)), $(is)), $(sols)) >> in
  let inits = << [op($(inits)), op($(sols))] >> in

  (* If we give additional initial values, display them and say where they *)
  (* come from. *)
  let u0orj =
    if <:bool< $(t) = 1 >>
    then << $(u)[0]($(n)) >>
    else << $(u)[j]($(n)) >> in
  let inits_text3 = if <:bool< nops($(sols)) = 0 >> then <:par<>> else
    <:par<
      But for the sake of completeness and convenience, all initial values
      <:isymb< $(symb: u0orj) >>,
    >> @:@
    (if <:bool< $(t) = 1 >>
    then <:par<>>
    else <:par<<:imath< 0 \leq j \leq << $(t) - 1 >>>> and >>)
    @:@
    <:par<<:imath< 0 \leq << $(n) >> < $(symb: q) >> must be given>> @:@
    (if <:bool< $(t) = 1 >> then <:par<. >> else
      <:par<,
        as well as those <:isymb< $(symb: u0orj) >> for which the leading
        coefficient of the respective recurrence vanishes.
      >>) @:@
    <:par<
      The missing values, however, can be obtained from the previous ones
      using the recurrence$(str: Wording.ending_of_seq recsys)
      and the assumption that
      <:imath< $(symb: u0orj) = 0 >> for <:imath< $(symb: n) < 0 >>.
      This produces <:dmath< << op($(sols)) >>. >>
    >> in

  (* Discard sequences whose initial values are all 0. *)
  let oldt = t
  and inits, recsys, t = ref inits, ref recsys, ref t in
  while
    <:bool<
      {0, op(map(a -> op(2, a), select(has, $(!inits), $(u)[$(!t)-1])))} = {0}
    >>
  do
    inits := << remove(has, $(!inits), $(u)[$(!t)-1]) >>;
    recsys := << eval(subsop(-1 = NULL, $(!recsys)), $(u)[$(!t)-1] = 0) >>;
    t := << $(!t) - 1 >>
  done;
  let inits, recsys, t = !inits, !recsys, !t in
  let zero_seq_text = if <:bool< $(t) = $(oldt) >> then <:par<>> else
    <:par<It turns out that all initial values of the sequence>> @:@
    <:par<$(str: Wording.ending_of_int <:int< $(oldt) - $(t) >>) >> @:@
    (Wording.enumeration_of_maplelist
       << [seq($(u)[i]($(n)), i = $(t) .. $(oldt)-1)] >>) @:@
    <:par< are zero. Hence >> @:@
    (if <:bool< $(oldt) - $(t) = 1 >>
    then <:par<this sequence is>>
    else <:par<these sequences are>>)
    @:@
    <:par< identically zero and can be discarded. >> in

  let solution = <<
    Record(
     'exp_factor' = 1,
     'alpha' = $(n0),
     'finite_part' = 0,
     'closed_form' = undefined,
     'variable' = $(x),
     'ramification' = 1,
     'log_powers' = [seq(i, i = 0 .. $(t)-1)],
     'sequence_name' = $(u),
     'index_name' = $(n),
     'recurrences' = ListTools:-Reverse($(recsys)),
     'initial_conditions' = sort($(inits))) >> in

  intro @:@ t_text @:@ recsys_text @:@ inits_text1 @:@
  inits_text2 @:@ inits_text3 @:@ zero_seq_text,
  solution



(* Compute all power series solutions of eqn at the origin of the form *)
(*   \sum_{n=0}^\infty \sum_{j=0}^{t-1} u_j(n) \ln(x)^j x^n            *)
(* where it is assumed that 0 is the smallest integer root             *)
(* of the indicial equation.                                           *)
(* The names u and n have to be given as arguments.                    *)
(* Equation numbers and notation are according to                      *)
(* Joris van der Hoeven: Fast Evaluation of Holonomic Functions Near   *)
(*   and in Regular Singularities (JSC 31, 717-743, 2001).             *)
let doc_obj eqn y x u n =

  let ord = << DDMF:-ode_order($(eqn)) >> in

  (* Transform eqn into Euler operator notation delta = x*D_x. *)
  (* In view of the correspondence (delta(f))_n = n*u(n) we use delta = n. *)
  let eqn_euler =
    << subs([$(y)($(x)) = 1, seq(diff($(y)($(x)), $(x) $ j) =
      $(x)^(-j) * add(combinat[stirling1](j, i) * $(n)^i, i = 1..j),
      j = 1 .. $(ord))], $(eqn)) >> in
  let eqn_euler =
    << expand($(x)^(-ldegree(expand($(eqn_euler)), $(x))) * $(eqn_euler)) >> in
  let eu = << theta[$(x)] >> in
  let euler_text =
    <:par<
      It is obtained by rewriting the
      differential equation into operator notation
      in terms of the Euler operator
      <:imath< << $(eu)(f($(x))) >> = $(symb: x)f'($(symb: x)) >>:
      <:dmath< << collect(subs($(n) = $(eu), $(eqn_euler)), $(eu), factor) >>.>>
    >> in

  (* Compute the polynomials Q that are introduced in Formula (3.7) *)
  (* and their derivatives: Q[i+1,j+1] corresponds to Q_i^{(j)}. *)
  (* Q_0 = Q[1,1] is the indicial polynomial. *)
  let bigQ =
    << [seq(collect(subs($(n) = $(n) - i, coeff($(eqn_euler), $(x), i)), $(n)),
        i = 0..degree($(eqn_euler), $(x)))] >> in
  let bigQ =
    << map(a -> [seq(diff(a, [$(n) $ j]), j = 0 .. degree($(bigQ)[1], $(n)))],
      $(bigQ)) >> in
  let q = << nops($(bigQ)) - 1 >> in
  let qs =
    List.map
      (fun a ->
        <:par<<:dmath<
          Q_$(int: a-1)($(symb: n)) = << $(bigQ)[$(int: a), 1] >>,
        >>>>)
      (CommonTools.range <:int< $(q) + 1 >>) in
  let bigQ_text =
    <:par<
      From this form, it can be easily read off how the differential equation
      acts on the Taylor coefficients of
      <:imath< f($(symb: x)) = \sum_{$(symb: n)=0}^\infty
        f_$(symb: n) $(symb: x)^$(symb: n) >>.
      The rules
      <:dmath< [$(symb:x)^$(symb:n)] \, ($(symb:eu) f($(symb:x))) =
        $(symb:n) f_$(symb:n), \quad
        [$(symb:x)^$(symb:n)] \, ($(symb:x) f($(symb:x))) = f_{$(symb:n)-1} >>
      imply that the <:isymb< $(n) >>th Taylor coefficient of the series
      obtained by substituting <:imath< f($(symb: x)) >> into the differential
      equation is given by
    >> @:@
    (match <:int< $(q) >> with
    | 0 -> <:par<<:dmath< Q_0($(symb: n)) f_$(symb: n) >>>>
    | 1 -> <:par<<:dmath< Q_0($(symb: n)) f_$(symb: n) +
          Q_1($(symb: n)) f_{$(symb: n)-1} >>>>
    | 2 -> <:par<<:dmath< Q_0($(symb: n)) f_$(symb: n) +
          Q_1($(symb: n)) f_{$(symb: n)-1} +
          Q_2($(symb: n)) f_{$(symb: n)-2} >>>>
    | _ -> <:par<<:dmath< Q_0($(symb: n)) f_$(symb: n) + \dots +
          Q_$(symb: q)($(symb: n)) f_{$(symb: n)-$(symb: q)} >>>>)
    @:@
    <:par<
      for some polynomials <:imath< Q_j($(symb: n)) >>,
      <:imath< 0 \leq j \leq $(symb: q) >>.
      (For the notation <:imath< [$(symb:x)^$(symb:n)] \, (\dots) >>, see
      $(t_ent:Glossary.g "formal power series").)
      Each <:imath< Q_j($(symb: n)) >>
      is obtained as the coefficient of <:imath< $(symb: x)^j >> in the
      differential operator above, where <:isymb< $(eu) >> has been replaced
      by <:imath< $(symb: n) - j >>, leading to:
    >> @:@
    (List.fold_left (@:@) <:par<>> qs) @:@
    <:par<
      and it is seen in this way \
      that the indicial polynomial is <:imath< Q_0 >>.
    >> in

  (* The integer roots of the indicial equation and their multiplicities. *)
  let all_roots = << [solve($(bigQ)[1,1] = 0, $(n))] >> in
  let int_roots = << sort(select(type, $(all_roots), integer)) >> in
  let int_roots = << DDMF:-equivalence_classes($(int_roots), `=`) >> in
  let mults = << map(nops, $(int_roots)) >>
  and int_roots = << map(a -> op(1, a), $(int_roots)) >>
  and sorted_roots = << convert(sort(convert($(all_roots), set)), list) >> in

  (* Display the indicial equation and its root(s). *)
  let ind_text =
    <:par<
      The indicial equation is
      <:dmath< << $(bigQ)[1,1] >> = 0 , >>
      whose
    >> @:@
    (if <:bool< nops($(sorted_roots)) = 1 >>
    then <:par<only root is >>
    else <:par<roots are >>)
    @:@ (Wording.enumeration_of_maplelist sorted_roots) @:@ <:par<.>> in

  (* All pairs [k, d] of indices for which a canonical solution f^{[n0,d]} *)
  (* exists where n0 = int_roots[k]. *)
  let all_cases = CommonTools.symb_list_of_symb
    <<
      [seq(seq([k, d], d = 0 .. $(mults)[k] - 1), k = 1 .. nops($(int_roots)))]
    >>
  in

  let contents, sols =
    List.split
      (List.map (deal_with_case x u n bigQ q int_roots mults) all_cases) in

  (* Create subsections only if there are several solutions. *)
  let sec_titles =
    List.map
      (fun a -> <:text<The canonical solution
        <:imath< f^{[ << $(int_roots)[$(a)[1]] >>, << $(a)[2] >> ]} >>>>)
      all_cases in
  let content =
    if List.length sols = 1
    then List.hd contents
    else List.fold_left
        (@@@)
        DC.ent_null
        (List.map2 DC.section sec_titles contents) in

  (ind_text @:@ euler_text @:@ bigQ_text) @@@ content,
  CommonTools.symb_of_symb_list sols


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