Sollya User's Manual

Name:

taylorform computes a rigorous polynomial approximation (polynomial, interval error bound) for a function, based on Taylor expansions.

Usage:

taylorform(f, n, x0, I, errorType) : (function, integer, constant, range, absolute|relative) -> list taylorform(f, n, x0, I, errorType) : (function, integer, range, range, absolute|relative) -> list taylorform(f, n, x0, errorType) : (function, integer, constant, absolute|relative) -> list taylorform(f, n, x0, errorType) : (function, integer, range, absolute|relative) -> list

Parameters:

f is the function to be approximated.
n is the order of the Taylor form, meaning n-1 is the degree of the polynomial that must approximate f.
x0 is the point (it can be a real number or an interval) where the Taylor exansion of the function is to be considered.
I is the interval over which the function is to be approximated. If this parameter is omitted, the behavior is changed (see detailed description below).
errorType (optional) is the type of error to be considered. See the detailed description below. Default is absolute.

Description:

taylorform computes an approximation polynomial and an interval error bound for function f. More precisely, it returns a list L=[p, coeffErrors, Delta] where:
- p is an approximation polynomial of degree n-1 which is roughly speaking a numerical Taylor expansion of f at the point x0.
- coeffsErrors is a list of n intervals. Each interval coeffsErrors[i] contains an enclosure of all the errors accumulated when computing the i-th coefficient of p.
- Delta is an interval that provides a bound for the approximation error between p and f. Its significance depends on the errorType considered.
The polynomial p and the bound Delta are obtained using Taylor Models principles.
Please note that x0 can be an interval. In general, it is meant to be a small interval approximating a non representable value. If x0 is given as a constant expression, it is first numerically evaluated (leading to a faithful rounding xapprox0 at precision prec), and it is then replaced by the (exactly representable) point-interval [xapprox0, xapprox0]. In particular, it is not the same to call taylorform with x0 = pi and with x0 = [pi], for instance. In general, if the point around which one desires to compute the polynomial is not exactly representable, one should preferably use a small interval for x0.
More formally, the mathematical property ensured by the algorithm may be stated as follows. For all xi0 in x0, there exist (small) values eps[i] in coeffsErrors[i] such that:
If errorType is absolute, for all x in I, there exists delta in Delta such that f(x) - p(x-xi0) = sum{i=0...n-1} eps[i]*(x-xi0)^i + delta.
If errorType is relative, for all x in I, there exists delta in Delta such that f(x) - p(x-xi0) = sum{i=0...n-1} eps[i]*(x-xi0)^i + delta*(x-xi0)^n.
It is also possible to use a large interval for x0, though it is not obvious to give an intuitive sense to the result of taylorform in that case. A particular case that might be interesting is when x0=I in relative mode. In that case, denoting by p_i the coefficient of p of order i, the interval p_i + coeffsError[i] gives an enclosure of f^(i)(I)/i!. However, the command autodiff is more convenient for computing such enclosures.
When the interval I is not given, the approximated Taylor polynomial is computed but no remainder is produced. In that case the returned list is L=[p, coeffErrors].
The relative case is especially useful when functions with removable singularities are considered. In such a case, this routine is able to compute a finite remainder bound, provided that the expansion point given is the problematic removable singularity point.
The algorithm does not guarantee that by increasing the degree of the approximation, the remainder bound will become smaller. Moreover, it may even become larger due to the dependecy phenomenon present with interval arithmetic. In order to reduce this phenomenon, a possible solution is to split the definition domain I into several smaller intervals.
The command taylor also computes a Taylor polynomial of a function. However it does not provide a bound on the remainder. Besides, taylor is a somehow symbolic command: each coefficient of the Taylor polynomial is computed exactly and returned as an expression tree exactly equal to theoretical value. It is henceforth much more inefficient than taylorform and taylorform should be prefered if only numercial (yet safe) computations are required. The same difference exists between commands diff and autodiff.

Example 1:

   > TL=taylorform(sin(x)/x, 10, 0, [-1,1], relative);
   > p=TL[0];
   > Delta=TL[2];
   > errors=TL[1];
   > for epsi in errors do epsi;
   [0;0]
   [0;0]
   [0;5.3455294201843912922810729343029637576303937602101e-51]
   [0;0]
   [-3.3409558876152445576756705839393523485189961001313e-52;3.3409558876152445576756705839393523485189961001313e-52]
   [0;0]
   [-1.04404871487976392427364705748104760891218628129103e-53;1.04404871487976392427364705748104760891218628129103e-53]
   [0;0]
   [-1.63132611699963113167757352731413688892529106451724e-55;1.63132611699963113167757352731413688892529106451724e-55]
   [0;0]
   [-1.91171029335894273243465647732125416670932546623114e-57;1.91171029335894273243465647732125416670932546623114e-57]
   > p; Delta;
   1 + x^2 * (-0.16666666666666666666666666666666666666666666666667 + x^2 * (8.3333333333333333333333333333333333333333333333333e-3 + x^2 * (-1.98412698412698412698412698412698412698412698412698e-4 + x^2 * (2.75573192239858906525573192239858906525573192239859e-6 + x^2 * (-2.50521083854417187750521083854417187750521083854419e-8)))))
   [-1.6135797443886066084999806203254010793747502812764e-10;1.6135797443886066084999806203254010793747502812764e-10]

Example 2:

   > TL=taylorform(exp(x), 10, 0, [-1,1], absolute);
   > p=TL[0];
   > Delta=TL[2];
   > p; Delta;
   1 + x * (1 + x * (0.5 + x * (0.16666666666666666666666666666666666666666666666667 + x * (4.1666666666666666666666666666666666666666666666667e-2 + x * (8.3333333333333333333333333333333333333333333333333e-3 + x * (1.38888888888888888888888888888888888888888888888889e-3 + x * (1.98412698412698412698412698412698412698412698412698e-4 + x * (2.4801587301587301587301587301587301587301587301587e-5 + x * (2.75573192239858906525573192239858906525573192239859e-6 + x * 2.7557319223985890652557319223985890652557319223986e-7)))))))))
   [-2.31142719641187619441242534182684745832539555102969e-8;2.7312660755642474420206278018039434042553645532164e-8]

Example 3:

   > TL1 = taylorform(exp(x), 10, log2(10), [-1,1], absolute);
   > TL2 = taylorform(exp(x), 10, [log2(10)], [-1,1], absolute);
   > TL1==TL2;
   false

Example 4:

   > TL1 = taylorform(exp(x), 3, 0, [0,1], relative);
   > TL2 = taylorform(exp(x), 3, 0, relative);
   > TL1[0]==TL2[0];
   true
   > TL1[1]==TL2[1];
   true
   > length(TL1);
   3
   > length(TL2);
   2

Example 5:

   > f = exp(cos(x)); x0 = 0;
   > TL = taylorform(f, 3, x0);
   > T1 = TL[0];
   > T2 = taylor(f, 3, x0);
   > print(coeff(T1, 2));
   -1.35914091422952261768014373567633124887862354684999
   > print(coeff(T2, 2));
   -(0.5 * exp(1))

WARNING