TaylorTruncProof.ml

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

INCLUDE "preamble.ml"


let title _ = <:text<Taylor Polynomial Approximation>>


exception Warning of DC.sec_entities

let symbolic_r str =
  if Str.string_match (Str.regexp "\\.[0-9]+\\|[0-9]+\\.?[0-9]*") str 0 then
    << convert(parse($(str:str)), 'rational', 'exact'): >>
  else
    raise (Warning <:par< <:imath< r >> should be a positive number. >>)


let_service TaylorTruncProof
  (deq : diffeq maple)
  (y : name maple)
  (x : name maple)
  (notation : string)
  (r0 : string)
  (prec : int) :
  DC.sec_entities * int option with { title = title } =
  try

    (* for Arctan and its friends *)
    let deq = <<gfun:-diffeqtohomdiffeq($(deq), $(y)($(x)))>> in

    let rad = symbolic_r r0 in
    let infsing = << gfun:-NumGfun:-utilities:-diffeq_infsing($(deq),
                                                              $(y)($(x))) >>
    in
    if <:bool< signum($(rad) - abs($(infsing))) >= 0 >> then
      <:par< It is not easy to determine whether
        <:imath< $(str: notation) >> can be evaluated on the disk
        <:imath< |$(symb:x)| \leq $(symb:rad) >> using this power series
        expansion, since the said disk contains <:imath< $(symb:infsing) >>,
        which is a singular point of the differential equation above. >>,
      None
    else


      (* 1. Sketch of proof *)


      let intro_p1 =
        <:par< The aim is to find an integer <:imath< n >> such that the
          <:imath< n >>th order Taylor polynomial
          of <:imath< $(str: notation) >> approximates it uniformly on the
          disk <:imath< \mathopen|$(symb:x)\mathclose| \leq r = $(symb:rad) >>
          with absolute error less
          than <:imath< \epsilon = 10^{ $(int:-prec) } >>.
          To achieve this goal, one looks for a
          $(t_ent:Glossary.g "majorant series") of <:imath< $(str: notation) >>
          in order to bound the tail of its Taylor expansion. >> in

      (* list of coefficients + order *)
      let c = << polyapprox:-extractcoeffs($(deq), $(y)($(x))) >> in
      let coeffs = << $(c)[1] >> in
      let r = <:int< $(c)[2] >> in

      (* non-zero coefficients *)
      let nz = << polyapprox:-nzcoeffs($(coeffs), $(int: r)) >> in

      (* equation to display (without zero coefficients) *)
      let equa = << diff($(y)($(x)), [$(x)$$(int: r)]) = add(``($(coeffs)[k]) *
                                      diff($(y)($(x)), [$(x)$k]), k in $(nz)) >>
      in

      let intro_p2 =
        <:par< The starting point is the differential equation
          <:dmath< <<$(equa)>>, >>
          satisfied by <:imath< $(str: notation) >>, together with initial
          conditions
          <:dmath< <<seq((D@@i)($(y))(0) = subs(select(type, $(deq), `=`),
            (D@@i)($(y))(0)) ,i=0..$(int:r-1))>>. >> >> in

      let intro_p3 =
        if <:int< nops(nz) >> = 1 then
          <:par< In the following section, we find a majorant
            series <:isymb< a[op(1,$(nz))]($(x)) >> of the
            coefficient <:isymb< $(coeffs)[op(1,$(nz))] >> of the equation.>>
        else
          <:par< In the following section, we find majorant series
            <:isymb< seq(a[k]($(x)), k in $(nz)) >> of the respective
            coefficients <:isymb< seq($(coeffs)[k], k in $(nz)) >>
            of the equation. >> in

      (* majorant equation (without zero coefficients) *)
      let majequa =
        << diff(phi($(x)), [$(x)$$(int: r)]) =
             add(a[k]($(x)) * diff(phi($(x)), [$(x)$k]), k in $(nz)) >>
      in

      let intro_p4 =
        <:par< The approximation method is then based on the following remark:
          if <:imath< \phi($(symb:x)) >> is a solution of the majorant equation
          <:dmath< $(symb:majequa), >> and if
          <:dmath< <<
            seq(abs((D@@i)($(y))(0)) <= (D@@i)(phi)(0), i = 0 .. $(int: r-1))
          >>, >>
          then <:imath< \phi($(symb:x)) >> is a majorant series of
          <:imath< $(str: notation) >>. >> in

      let sketch = DC.section
          <:text<Sketch of proof>>
          (intro_p1 @@@ intro_p2 @@@ intro_p3 @@@ intro_p4) in


      (* 2. Majorant series of the coefficients *)


      (* compute a radius 1/a between r and the distance to the first
       * singularity + and evaluation of the (relative) error *)
      let sing = << polyapprox:-boundeqsing($(coeffs), $(x), $(int:r),$(rad)) >>
      in

      let a = << polyapprox:-pretty_print($(sing)[1]) >>
      and g = << polyapprox:-pretty_print($(sing)[2]) >> in

      (* majorant series for the coefficients *)
      let ms  = << polyapprox:-boundcoeffs($(coeffs),$(x),$(int:r),$(g),$(a))>>
      in

      (* display them *)

      let majorants = ref [] in

      for i = 0 to r - 1 do

        if <:bool< evalb($(coeffs)[$(int:i)] <> 0) >> then

          (* link to the proof for this majorant series *)

          let ldescr = RatMajSeries.descr
              (<:symb< $(coeffs)[$(int:i)] >>, x, g, a, r - i) None in
          let rlink = DC.link_service ldescr <:text<here>> in

          (* append it to the list *)

          let mi = << polyapprox:-pretty_print($(ms)[$(int:i)]) >> in

          let thismaj =
            <:par< <:imath< a_{$(int:i)}($(symb:x)) =
                \frac{M_{$(int: i)}}{(1 - \alpha $(symb:x))^{$(int: r - i)}} >>,
              where <:imath< M_{$(int:i)} = $(symb:mi) >>, is a majorant series
              of <:isymb< $(coeffs)[$(int:i)] >>. More details >>
          in

          majorants :=
            (thismaj @:@ <:par<$(rlink)>> @:@ <:par<.>>) :: !majorants

    done;

    let introtext =
      if <:bool< evalb($(infsing) = infinity) >> then
        <:par< Since the coefficients of the differential equation we started
        from have no poles, the function <:imath< $(str: notation) >> has no
        finite singularity, and the radius of convergence of its power series
        expansion is infinite. This allows us to look for majorant series
        of the form <:imath< \frac M {(1 - \hat\alpha $(symb:x))^m}>>, where
        <:imath< \hat\alpha < r^{-1} = <:symb< 1/$(rad) >> >>
        and <:imath< m >> can be chosen as we
        please. Somewhat arbitrarily, we set <:isymb< alpha = $(a) >>. Then: >>
      else
        <:par< Now let <:imath< \hat\alpha >> denote the inverse of the
        minimum modulus of a singularity of <:imath< $(str: notation) >>. One
        can show that the exponential behaviour of the Taylor coefficients of
        a (generic) solution of the differential equation satisfied
        by <:imath< $(str: notation) >> is controlled
        by <:imath< \hat\alpha >>. This invites to look for majorant series
        of the form <:imath< \frac M {(1 - \hat\alpha $(symb:x))^m}>>.
        It is not easy to compute <:imath< \hat\alpha >>, but using the
        Dendelin-Graeffe iteration, one can compute an
        approximation <:imath< \alpha > \hat\alpha >>,
        namely <:isymb< alpha = $(a) >>. Then: >>
    in

    let majseries =
      DC.section
        <:text<Majorant Series of the Coefficients>>
        (introtext @@@ DC.unordered_list !majorants)
    in


    (* 3. Majorant equation and its resolution *)


    (* majorant equation (without zero coefficients) *)
    let majeq_symb =
      << diff(phi($(x)), [$(x)$$(int: r)]) =
         add(``(M[k] /
           (1 - alpha * $(x)) ^ ($(int:r) - k)) * diff(phi($(x)), [$(x)$k]),
           k in $(nz)) >>
    in

    let majeq_p1 =
      <:par< According to the results of the previous sections, the following
        equation is a majorant equation for <:imath< $(str: notation) >>:
        <:dmath< $(symb: majeq_symb). >> >>
    in

    (* majorant series for our function *)
    let maj = << polyapprox:-bounddffun($(deq), $(y)($(x)), $(g), $(a)) >> in
    let cstA = << polyapprox:-pretty_print($(maj)[1]) >> in
    let lambda = << polyapprox:-pretty_print($(maj)[2]) >> in

    let majeq_p2 =
      <:par< This is an Euler equation, which admits a solution in the form
        <:dmath< \frac{\hat A}{(1 - \alpha $(symb:x))^{\hat\lambda}}, >> where
        <:imath< \hat\lambda >> is the unique positive solution of the equation
        <:isymb< ($(a) * $(g)) ^ $(int: r) * expand(pochhammer(X, $(int: r))) -
          add(($(a) * $(g)) ^ j * expand(pochhammer(X, j)),
            j = 0 ..($(int: r) - 1)) = 0 >>. >> @@@

      <:par< Hence, according to the result stated in the first section, it
        suffices to choose <:imath< \hat A >> so that the inequalities on the
        initial conditions hold. Numerically, we get
        <:imath< \hat A \leq $(symb: cstA) =: A>> and
        <:imath< \hat\lambda \geq $(symb: lambda) =: \lambda >>. >> @@@

      <:par< The function <:imath< \tilde\phi($(symb:x)) >> defined by
        <:dmath< \tilde\phi($(symb:x)) =
          \frac A {(1 - \alpha $(symb:x))^\lambda}, >>
        is then a majorant series of
        <:imath<\frac{\hat A}{(1 - \alpha $(symb:x))^{\hat\lambda}}>>,
        which is itself a majorant series of <:imath< $(str: notation) >>. >>
    in

    let majeq =
      DC.section
        <:text< Resolution of the majorant equation>>
        (majeq_p1 @@@ majeq_p2) in


    (* 4. Bound on the tail of the majorant series *)


    (* truncation order, first estimation *)
    let order =
      <:int< polyapprox:-taylortrunc($(deq), $(y)($(x)), $(g), $(a),
                                     $(rad), 10 ^ (-$(int: prec)) / 2) >>
    in

    let tailbound = DC.section <:text< Bound on the tail>>
      (<:par< Note that the remainder <:imath< R_n($(symb:x)) >> of the Taylor
         expansion of <:imath< $(str: notation) >> is bounded by the remainder
         <:imath< \tilde R_n(|$(symb:x)|) >> of the Taylor expansion of
         <:imath< \tilde\phi(|$(symb:x)|) >> for every <:imath< n >>
         and <:imath< $(symb:x) >>. We will look for an integer <:imath< n >>
         such that <:imath< \tilde R_n(|$(symb:x)|) < \epsilon / 2 >> on the
         disk <:imath< \mathopen|$(symb:x)\mathclose| \leq r = $(symb:rad) >>,
         with <:imath< \epsilon = 10 ^{-$(int: prec)} >>. >> @@@

      <:par< Let <:imath< \eta = \frac 1 2 \left(r + \frac 1 \alpha\right) >>
        and <:imath< M = \max_{\theta \in [0, 2\pi]} |\phi(\eta e^{i\theta})|>>.
        It is easy to prove that
        <:imath< M = \frac A { | 1 - \alpha \eta | ^ \lambda } >>.
        Then, applying the Cauchy inequality, we have:
        <:dmath< |R_{n-1}($(symb:x))| \leq \tilde R_{n-1}(|$(symb:x)|) \leq
          \sum_{k \geq n} \frac M {\eta^k} |$(symb:x)|^k \leq
          M \left(\frac {|$(symb:x)|} \eta\right)^n
          \frac 1 {1-\frac {|$(symb:x)|} \eta}, >>
        whenever <:imath< |$(symb:x)| \leq r < \eta >>. Thus, if
        <:dmath< n \geq $(int: order), >>
        then
        <:dmath< |R_{n-1}($(symb:x))| \leq M \left(\frac r \eta\right)^n
          \frac 1 {1 - \frac r \eta} \leq \epsilon / 2. >> >>)
      in


    (* 5. Degree improvement *)


    (* much better truncation order *)
    let improved = <:int< polyapprox:-optimize($(deq), $(y)($(x)),
$(int: order), $(rad), 10 ^ (-$(int: prec))) >> in

    let optim =
      DC.section
        <:text<Degree improvement>>
        (<:par< Let
           <:dmath< \sum_{n \geq 0} u(n) $(symb:x)^n >>
           be the Taylor expansion of <:imath< $(str: notation) >>. We have
           proved that
           <:dmath< \left| \sum_{n=0}^{$(int: order)}
             u(n) $(symb:x)^n - $(str: notation) \right| <
             \epsilon / 2, >>
           for <:imath< \mathopen| $(symb:x) \mathclose| \leq r =
           $(symb:rad) >> and <:imath< \epsilon = 10 ^{-$(int: prec)} >>. >> @@@

      <:par< Computing explicit values for the coefficients <:imath<u(n)>>
        by the recurrence relation, one gets, numerically,
        <:dmath<
          \sum_{k = $(int: improved + 1)}^{$(int: order)} |u(n)| r^n
          \leq \epsilon / 2, >>
        so that finally, by the triangular inequality
        <:dmath<
          \left|\sum_{n=0}^{$(int: improved)} u(n) $(symb:x)^n -
            $(str: notation) \right| < \epsilon, >>
        whenever
        <:imath< \mathopen| $(symb:x) \mathclose| \leq r = $(symb:rad) >> and
        <:imath< \epsilon = 10 ^{-$(int: prec)} >>. >>) in

    DC.section
      (title ()) (sketch @@@ majseries @@@ majeq @@@ tailbound @@@ optim),
    Some improved


  with (Warning msg) -> (DC.warning msg, None)

Generated by GNU Enscript 1.6.5.90.