LocalExpansion.ml

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

INCLUDE "preamble.ml"

(* Some exceptions used in this module. *)
exception BadTransf
exception BadInits

let title (_, point, _, _, _) = <:text<Expansion at <:isymb< $(point) >>>>


let contribution_from_alpha basis var pfac transf monomials =
  (* Identify basis elements that are compatible with pfac and alpha. *)
  let basis = CommonTools.symb_list_of_symb basis in
  let e1 = (List.hd monomials).SF.x_expon in
  let match_pfac b =
    <:bool<
      proc(b, v) local pfac_b, e, s;
        pfac_b := DDMF:-transform_back(b:-exp_factor,
          DDMF:-dynamow_to_maple(b:-variable), v);
        pfac_b := pfac_b * v ^ b:-alpha;
        test := combine($(pfac) * v ^ $(e1) / pfac_b, 'power', 'symbolic');
        return match(test = v^e, {v}, s) and
          type(b:-ramification * subs(s, e), integer);
      end proc($(b), $(var))
    >> in
  let basis = List.filter match_pfac basis in

  (* Investigate the range of the x_expon's that occur in expansion. *)
  (* This is needed to determine the order of the truncations later. *)
  let x_expons =
    CommonTools.symb_of_symb_list
      (List.map (fun a -> a.SF.x_expon) monomials) in
  let e1 = << $(x_expons)[1] >> in
  let min_expon = << $(e1) + min(map(a -> a - $(e1), $(x_expons))) >> in
  let max_diff = << max(map(a -> a - $(min_expon), $(x_expons))) >> in

  (* Compute truncated version of the basis elements up to max_diff. *)
  let truncate b =
    (* We need to correct max_diff by possible contributions that come *)
    (* from prefactor and alpha, since the entries in the database are *)
    (* not consistent in where to store powers of x (prefactor or x_expon). *)
    let x = << DDMF:-dynamow_to_maple($(b):-variable) >> in
    let pfac_b = << DDMF:-transform_back($(b):-exp_factor, $(x), $(var)) >> in
    let pfac_b = << $(pfac_b) * $(var) ^ $(b):-alpha >> in
    let corr1 =
      <<
        proc(expr, x) local e, s;
          match(expr = x^e, {x}, s);
          return subs(s, e);
        end proc(combine($(pfac) * $(var) ^ $(min_expon) / $(pfac_b)), $(var))
      >> in
    let corr2 = try <:int< $(b):-alpha >> with _ -> 0 in
    let ord = <:int< floor($(max_diff) + $(corr1)) >> + corr2 + 1 in
    let trunc =
      TruncatedSeries.obj
        ((<< DDMF:-encode_linear_combination([[1, $(b)]]) >>, ""), ord) in
    let trunc = << DDMF:-dynamow_to_maple($(trunc)) >> in
    let trunc =
      << DDMF:-transform_back(eval($(trunc), O = 0), $(x), $(var)) >> in
    << expand($(trunc) / $(pfac)) >> in
  let trunc_basis = List.map truncate basis in

  (* Extract the coefficient of x^x_expon * ln(x)^ln_expon from expr, *)
  (* where x_expon is some algebraic expression in the parameters and *)
  (* ln_expon is a nonnegative integer. *)
  let coeff_x_ln x x_expon ln_expon expr =
    let coeff_ln = << coeff($(expr), ln($(x)), $(ln_expon)) >> in
    << DDMF:-coeff_of_power_sum($(coeff_ln), $(x), $(x_expon)) >> in

  (* Extract the equation that corresponds to the particular initial *)
  (* condition given by mon, as a matrix row. *)
  let get_row mon =
    CommonTools.symb_of_symb_list
      ((List.map (coeff_x_ln var mon.SF.x_expon mon.SF.ln_expon) trunc_basis) @
       [ << -$(mon.SF.coeff) >> ]) in

  let mat =
    CommonTools.symb_of_symb_list (List.map get_row monomials) in
  let ns = << LinearAlgebra[NullSpace](Matrix($(mat))) >> in
  if <:bool< nops($(ns)) <> 1 or $(ns)[1][-1] = 0 >> then raise BadInits;
  let ns1 = << convert($(ns)[1], list) >> in
  let coeffs = << map(a -> a / $(ns1)[-1], subsop(-1 = NULL, $(ns1))) >> in
  << zip(`[]`, $(coeffs), $(CommonTools.symb_of_symb_list basis)) >>


let contribution_from_expansion basis var expansion =
  let mons = expansion.SF.monomials in
  let x_expons = List.map (fun a -> a.SF.x_expon) mons in
  let ram = << if $(basis) = [] then 1 else $(basis)[1]:-ramification end if >>
  and x_expons = CommonTools.symb_of_symb_list x_expons in
  let eqv =
    <<
      DDMF:-equivalence_classes($(x_expons),
        (a, b) -> type($(ram) * (a - b), integer))
    >> in
  let reps = CommonTools.symb_list_of_symb << map(a -> a[1], $(eqv)) >> in
  let monomials_in_class rep =
    List.filter
      (fun a -> <:bool< type($(ram) * ($(rep) - $(a.SF.x_expon)), integer) >>)
      mons in
  let contributions =
    List.map
      (contribution_from_alpha
         basis var expansion.SF.prefactor expansion.SF.transf)
      (List.map monomials_in_class reps) in
  << map(op, $(CommonTools.symb_of_symb_list contributions)) >>


let majorant_series notation def y x =
  let maj_series = << gfun:-NumGfun:-bound_diffeq($(def), $(y)($(x)), N) >> in
  let coeftayls = << selectfun($(maj_series), 'Coeftayl') >> in
  <:par<
    A $(t_ent:Glossary.g "majorant series") for <:imath< $(str: notation) >>
    in 0 is given by
  >> @:@
  if <:bool< $(coeftayls) <> {} >> then
    <:par<
      <:dsymb< subs(op($(coeftayls))=u[n], $(maj_series)) >>
      where <:imath< u_n >> denotes the coefficient of <:imath< z^n >> in
      the series expansion at the origin of
      <:dmath<<< op(1, op($(coeftayls))) >>.>>
    >>
  else <:par<<:dmath< << subs(op($(coeftayls))=u[n], $(maj_series)) >> . >>>>


let_service LocalExpansion
  (marshaled_sf : string)
  (point : any maple)
  (y : name maple)
  (notation : string)
  (parametrized : bool) :
  DC.sec_entities * any maple with {title = title} =

  let err = << FAIL >> and t = "No" in
  let content, result, terminology =
    try

      let sf = DB.unmarshal_sf marshaled_sf in
      let var = SF.var_of_t sf
      and ics = SF.get_initial_conditions sf point in
      let analytic = ics.SF.analytic in
      let terminology = if analytic then "Taylor" else "local"
      and def = Definition.definition sf y point in

      (* Check whether all expansions are given in terms of the same transf. *)
      let transfs = List.map (fun e -> e.SF.transf) ics.SF.expansions in
      if <:bool< nops({op($(CommonTools.symb_of_symb_list transfs))}) <> 1 >>
      then raise BadTransf;
      let transf = List.hd transfs in

      (* Find names for sequences and summation indices that do not conflict *)
      (* with the choices of the parameters. For this purpose, all symbolic  *)
      (* expressions in ics and the ode are assembled as testing expression. *)
      let symb_in_mon m = [m.SF.x_expon ; m.SF.ln_expon ; m.SF.coeff] in
      let symb_in_exp e =
        [e.SF.prefactor ; e.SF.transf] @
        (List.flatten (List.map symb_in_mon e.SF.monomials)) in
      let symb_in_ic i =
        [i.SF.point] @ (List.flatten (List.map symb_in_exp i.SF.expansions)) in
      let symbs = CommonTools.symb_of_symb_list (symb_in_ic ics) in
      let test_expr = << [ $(sf.SF.lode), $(symbs) ] >> in
      let seq_var, sum_var = CommonTools.unused_names_for_sequence test_expr in

      let basis =
        LocalSolutionsPoint.obj
          ((sf.SF.lode, y, var, transf, seq_var, sum_var), ()) in

      let conts =
        List.map (contribution_from_expansion basis var) ics.SF.expansions in
      let linear_comb = << map(op, $(CommonTools.symb_of_symb_list conts)) >> in
      let linear_comb = << select(a -> op(1,a) <> 0, $(linear_comb)) >> in

      (* Make all sequence names different. *)
      <:unit<
        proc(u) local i, ser;
          for i from 1 to nops($(linear_comb)) do
            ser := $(linear_comb)[i,2];
            ser:-sequence_name := u[i];
            ser:-recurrences := subs(u = u[i], ser:-recurrences);
            ser:-initial_conditions := subs(u = u[i], ser:-initial_conditions);
          end do;
        end proc($(seq_var))
      >>;

      (* For the following service calls, we need the encoded version. *)
      let linear_comb = << DDMF:-encode_linear_combination($(linear_comb)) >> in

      (* Expansion (in closed form, if possible). *)
      let closed_form, linear_comb =
        ClosedForm.obj ((linear_comb, notation), ()) in
      let cf_text =
        DC.inline_service (ClosedForm.descr (linear_comb, notation) None) in

      (* Recurrence relations. *)
      let funcs =
        << indets(DDMF:-dynamow_to_maple($(closed_form)), function) >> in
      let seqs = << select(a -> has(op(0, a), $(seq_var)), $(funcs)) >> in
      let s = Wording.ending_of_seq seqs in
      let rr_link =
        DC.link_service
          (RecurrenceRelations.descr
             (notation, linear_comb, def, point, y, var, seq_var,
              sum_var, analytic)
             None)
          <:text<recurrence relation$(str: s)>> in
      let rr_text =
        if <:bool< $(seqs) = {} >> then
          <:par<
            See the $(rr_link) for the coefficients of the
            $(str:terminology) expansion.
          >>
        else
          let a = if <:bool< nops($(seqs)) = 1 >> then "a " else "" in
          <:par<The coefficient sequence$(str: s) >> @:@
          (Wording.enumeration_of_maplelist << [op($(seqs))] >>) @:@
          <:par< can be computed by $(str: a)linear $(rr_link).>> in

      (* First terms. *)
      let decode = << DDMF:-decode_linear_combination($(linear_comb)) >> in
      let ft_text =
        if <:bool< has(map(a -> a[2]:-closed_form, $(decode)), Sum) >>
        then
          DC.inline_service
            (TruncatedSeries.descr (linear_comb, notation) None)
        else DC.ent_null in

      (* Taylor truncation. *)
      let tt_text =
        (* FIXME: This may in fact work in some non-analytic cases. *)
        if <:bool< $(point) = 0 >> && analytic && (not parametrized)
        then
          DC.inline_service
            (TaylorTruncation.descr (def, y, var, notation) None)
        else DC.ent_null in

      (* Majorant series. *)
      let ms_text =
        (* FIXME: This may in fact work in some non-analytic cases. *)
        if <:bool< $(point) = 0 >> && analytic && (not parametrized)
        then majorant_series notation def y var
        else DC.ent_null in

      (* Plot of Taylor expansions. *)
      let tp_text =
        if <:bool< $(point) = 0 >> && (not parametrized)
        then DC.inline_service (TaylorExpansionPlot.descr (sf.SF.rep, var) None)
        else DC.ent_null in

      DC.unordered_list
        [cf_text @@@ rr_text ; ft_text ; tt_text ; ms_text ; tp_text],
      linear_comb,
      String.capitalize terminology

    with
    | BadTransf ->
        DC.warning <:par<Different transformations for the same point.>>, err, t
    | BadInits ->
        DC.warning <:par<Bad initial conditions.>>, err, t
    | SF.UndefinedInitialCondition ->
        DC.warning <:par<No initial conditions found for this point.>>, err, t
    | _ ->
        DC.warning <:par<Some error occurred.>>, err, t
  in

  DC.section
    <:text<$(str: terminology) Expansion at <:isymb<$(point)>>>> content,
  result

Generated by GNU Enscript 1.6.5.90.