DEMOTAYLOR Short demonstration of Taylor toolbox

Some sample applications of the Taylor toolbox

The Taylor toolbox computes Taylor coefficients of a univariate function in forward mode, which is conveniently to use by the Matlab operator concept. It works much in the spirit of the gradient and Hessian toolbox, so I recommend to visit the gradient and Hessian demo first, see "demo toolbox intlab".

Initialization of Taylor variables

In order to use the automatic Taylor toolbox, the independent variable need to be identified and a value has to be assigned. This is performed by the function "taylorinit", for example

```format compact short _
u = taylorinit(2.89)
```
```Taylor u.t =
2.8900    1.0000         0         0         0
```

Operations between Taylor variables

If at least one operand is of type Taylor, operations are executed as Taylor operations. For example,

```x = taylorinit(3.5);
y = sin(3*x-sqrt(x+5))
```
```Taylor y.t =
0.9639    0.7530   -3.8545   -1.0178    2.5661
```

Note that first component of y is f(x), followed by f'(x), f''(x)/2!, f'''(x)/3!, etc.

Interval Taylor variables

If arguments are of type intval, an inclusion of the true value is computed:

```format short
y = sin(3*x-sqrt(x+5))
```
```intval Taylor y.t =
0.9639    0.7529   -3.8545   -1.0178    2.5661
```

For f(x):=exp(3*x-sqrt(x)), the result y contains in y.t the function value f(3.5) and the first 4 derivatives:

```y.t
```
```intval ans =
0.9639    0.7529   -3.8545   -1.0178    2.5661
```

Affari Taylor variables

If arguments are of type affari, an inclusion of the true value is computed using affine arithmetic:

```x = taylorinit(affari(midrad(3.5,1e-8)));
y = sin(3*x-sqrt(x+5))
```
```affari Taylor y.t =
0.9639    0.7529   -3.8545   -1.0178    2.5661
```

In this case the interval result is very narrow, so there is no difference.

Taylor vector variables

Note that the Taylor toolbox accepts one independent variable. One may initialize a Taylor variable of a vector argument; this is the same as initializing each component as the independent variable (with a different value). It is convenient for function evaluations with many arguments:

```f = inline('sin(3*x-sqrt(x+5))')
x = taylorinit([-3 0.1 3.5]')
y = f(x)
```
```f =
Inline function:
f(x) = sin(3*x-sqrt(x+5))
Taylor x.t =
-3.0000    1.0000         0         0         0
0.1000    1.0000         0         0         0
3.5000    1.0000         0         0         0
Taylor y.t =
0.8357   -1.4533   -2.9508    1.6048    1.8148
-0.9258   -1.0500    3.5700    1.3794   -2.2864
0.9639    0.7530   -3.8545   -1.0178    2.5661
```

Complex arguments

When evaluating the expression for another argument, use the same statement as before with new values. Here we assign the Taylor variable to carry 2 derivatives (the default is 4):

```x = taylorinit(-3.5+.2i,4);
y = sin(3*x-sqrt(x))
```
```Taylor y.t =
Columns 1 through 4
1.7385 + 0.7030i  -2.8592 + 4.2242i  -7.1868 - 4.5292i   5.4110 - 5.5865i
Column 5
4.8130 + 4.3813i
```

The Taylor coefficients are accessed by {}, so that y{0} is the function value and y{k} denotes the k-th Taylor coefficient f^(k)(x)/k! :

```y{0}
y{1:3}
```
```ans =
1.7385 + 0.7030i
ans =
-2.8592 + 4.2242i  -7.1868 - 4.5292i   5.4110 - 5.5865i
```

When initializing a Taylor vector, the individual vector components are accessed by () and Taylor coefficients by {}. For example,

```f = inline('sin(3*x-sqrt(x+5))')
x = taylorinit([-3 0.1 3.5]')
y = f(x)
y(1)
y{2}(3)
```
```f =
Inline function:
f(x) = sin(3*x-sqrt(x+5))
Taylor x.t =
-3.0000    1.0000         0         0         0
0.1000    1.0000         0         0         0
3.5000    1.0000         0         0         0
Taylor y.t =
0.8357   -1.4533   -2.9508    1.6048    1.8148
-0.9258   -1.0500    3.5700    1.3794   -2.2864
0.9639    0.7530   -3.8545   -1.0178    2.5661
Taylor ans.t =
0.8357   -1.4533   -2.9508    1.6048    1.8148
ans =
-3.8545
```

accesses the Taylor value f(-3) and the second Taylor coefficient of f(3.5), respectively.

An example: Taylor series

Define

```f = inline('sinh(x-exp(2/x))')
```
```f =
Inline function:
f(x) = sinh(x-exp(2/x))
```

Then the Taylor coefficients 0..7 of f at x=1.234 are computed by

```kmax = 7;
x = 1.234;
y = f(taylorinit(x,kmax))
```
```Taylor y.t =
1.0e+05 *
Columns 1 through 7
-0.0002    0.0017   -0.0089    0.0371   -0.1357    0.4525   -1.4057
Column 8
4.1252
```

The Taylor coefficients y{k} satisfy f(x+e) = sum[0..k]( y{k}*e^k ) + O(e^(k+1)) :

```format long
e = 1e-3;
v = f(x+e)
yapprox = sum( y{0:kmax} .* e.^(0:kmax) )
```
```v =
-22.682513324779098
yapprox =
-22.682513324779070
```

Inclusion of a function value by Taylor series

For an inclusion of the function value we may calculate the Taylor coefficients in interval arithmetic and add the error term:

```format long _
x = intval('1.234');
Y = f(taylorinit(x,kmax));
e = intval('1e-3');
Y_ = f(taylorinit(x+hull(0,e),kmax+1));
for k=0:kmax
Yincl = sum( Y{0:k} .* e.^(0:k) ) + Y_{k+1}*e.^(k+1)
end
```
```intval Yincl =
-22.68____________
intval Yincl =
-22.68251_________
intval Yincl =
-22.6825133_______
intval Yincl =
-22.682513325_____
intval Yincl =
-22.682513324779__
intval Yincl =
-22.682513324779__
intval Yincl =
-22.682513324779__
intval Yincl =
-22.682513324779__
```

Note how nicely the linear convergence can be observed by the "_"-notation. Also note that this is a true inclusion of f(1.234+1e-3)=f(1.235) because both arguments x=1.234 and e=1e-3 are intervals including the decimal numbers 1.234 and 0.001 (both are not floating-point numbers).

An Application: integration

Consider

```f = @(x)(sin(pi*x)-sin(x)); a = 0; b = 20;
x = linspace(a,b,1000); close, plot(x,f(x),x,0*x)
```

It is easy to see that for the transcendental number pi, the true value of the integral of f from a to b is cos(b)-1:

```cos(b)-1
```
```ans =
-0.591917938186608
```

There is a rudimentary integration routine "verifyquad" using Romberg's rule based on the Taylor toolbox. It calculates

```ApproxIncl = verifyquad(f,a,b)
infsup(ApproxIncl)
```
```intval ApproxIncl =
-0.591918________
intval ApproxIncl =
[  -0.59191828548002,  -0.59191759089296]
```

Integration using transcendental constants

This is a true inclusion of the integral with "pi" denoting the floating-point approximation of the transcendental number pi. To calculate an inclusion of the function with the true transcendental number pi, we use the following program [Note we would better use the inclusion intval('pi'), this is just an example]:

```% function y = testfuntaylor(x)
%   if isintval(x)
%     Pi = 4*atan(intval(1));
%   else
%     Pi = pi;
%   end
%   y = sin(Pi*x)-sin(x);
%
```

The result, however, does not change very much:

```TrueIncl = verifyquad(@testfuntaylor,a,b)
infsup(TrueIncl)
```
```intval TrueIncl =
-0.591918________
intval TrueIncl =
[  -0.59191828548006,  -0.59191759089293]
```

A comparison to the Matlab function "quad"

For this particular function the approximate routine may get problems if we specify a little more accuracy:

```e = 1e-12;
```
```Approx =
-0.591917938186609
Elapsed time is 0.048475 seconds.
intval Incl =
-0.591917938187__
Elapsed time is 0.054425 seconds.
```

Note that the verification routine is calculates an inclusion of the 'true' function (with the transcendental number pi). Insisting on even more accuracy make things worse:

```e = 1e-14;
```
```Approx =
0.984600203152193
Elapsed time is 0.037034 seconds.
intval Incl =
-0.591917938187__
Elapsed time is 0.072119 seconds.
```

Note that the approximate value has no correct digit, but the Matlab routine "quad" gives no error message.

Enjoy INTLAB

INTLAB was designed and written by S.M. Rump, head of the Institute for Reliable Computing, Hamburg University of Technology. Suggestions are always welcome to rump (at) tuhh.de