# 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.