# DEMOSTDFCTS Interval standard functions in INTLAB

## Contents

As explained in "DEMOINTVAL", standard functions with interval argument compute an inclusion of the true function value or the true function range.

## Input arguments I

Special care is necessary if an input argument is not a floating-point number. For example,

```format long _
x = intval([2.5 -0.3]);
y = erf(x)
```
```intval y =
0.99959304798255  -0.32862675945912
```

computes the true range of the error function erf(2.5), but not necessarily of erf(-0.3). To obtain true results use

```x = intval('2.5 -0.3');
y = erf(x)
```
```intval y =
0.99959304798255
-0.32862675945912
```

Note that intval converts a character string representing a vector always into a column vector. To see the accuracy of the result, the infsup-notation may be used:

```infsup(y)
relerr(y)
```
```intval y =
[   0.99959304798255,   0.99959304798256]
[  -0.32862675945913,  -0.32862675945912]
ans =
1.0e-15 *
0.055533750803183
0.422296341011770
```

## Input arguments II

Similarly, the range of a function is included by

```x = intval('[2.5,2.500001]')
y = erf(x)
```
```intval x =
2.500001________
intval y =
0.99959305______
```

Note that the large diameter of the output is due to the large diameter of the input.

## Accuracy of standard functions

The complementary error function erfc(x) is defined by 1-erf(x). It rapidly approximates 1 for larger x. In this case erfc(x) is much more accurate than 1-erf(x).

```x = intval(5);
y1 = erfc(x)
y2 = 1 - erf(x)
infsup([y1;y2])
```
```intval y1 =
1.0e-011 *
0.15374597944280
intval y2 =
1.0e-011 *
0.15375_________
intval  =
1.0e-011 *
[   0.15374597944280,   0.15374597944281]
[   0.15374368445009,   0.15375478668034]
```

Sometimes the built-in standard functions are not very accurate:

```x = linspace(1+1e-14,1+1e-5,100000);
y = acoth(x); Y = acoth(intval(x));
close, semilogy( x,abs((y-Y.mid)./(y+Y.mid)), x,relerr(Y) )
``` ## Nonlinear equations involving standard functions I

The zero of a nonlinear function may be approximated by some Newton procedure. Consider

```x = linspace(0,10);
f = @(x) erf(x)-sin(x)
close
plot(x,f(x),x,0*x)
```
```f =
function_handle with value:
@(x)erf(x)-sin(x)
``` The graph shows the function erf(x)-sin(x) between 0 and 10. Besides the obvious root zero there seem to be two small roots near 1 and 2, and mathematically it follows that there must be two roots near 8.

To approximate the roots, a Newton procedure may be used. For example

```x = 1;
for i=1:5
y = f(gradientinit(x));
x = x - y.x/y.dx
end
```
```x =
1.009823156064284
x =
1.009829234350912
x =
1.009829234354570
x =
1.009829234354571
x =
1.009829234354570
```

## Inclusion of the solution of nonlinear systems with standard functions

As expected, the iteration converges very rapidly. The final value is a very good approximation of a root of f, and an inclusion is computed by

```Y = verifynlss(f,x)
```
```intval Y =
1.00982923435457
```

## The range of a function and overestimation

The computed inclusion of the range of a function is a true estimate, but may be not sharp:

```x = infsup(0,1);
yinf = f(x)
```
```intval yinf =
[  -0.84147098480790,   0.84270079294972]
```

The inclusion can be improved using affine interval arithmetic:

```yaff = f(affari(x))
```
```affari yaff =
[  -0.05999375863532,   0.10111067175095]
```

## Nonlinear equations involving standard functions II

More interesting is the root cluster near x=8. An attempt to compute inclusions fails:

```Y1 = verifynlss(f,7.5)
Y2 = verifynlss(f,8.5)
```
```intval Y1 =
NaN
intval Y2 =
NaN
```

The reason is that numerically this is a double zero because erf(8) is very close to 1, as can be checked by erfc(8) = 1-erf(8) :

```y = erfc(intval(8))
```
```intval y =
1.0e-028 *
0.11224297172982
```

At least an inclusion of the extremum of f near 8 can be computed. Clearly it must be near to 2.5*pi.

```Y = verifylocalmin(f,8)
y = 2.5*pi
```
```intval Y =
7.85398163397448
y =
7.853981633974483
```

For interval enthusiasts an inclusion of 2.5*pi can be computed as well. Note, however, that this is only close to the extremum.

```y1 = 2.5*intval('pi')
```
```intval y1 =
7.85398163397448
```

## Verification of two roots

In the previous example, the function f has a nearly double root. It can be shown that a perturbed function g has a true double root:

```Y = verifynlss2(f,8);
root = Y(1)
perturbation = Y(2)
```
```intval root =
7.85398163397448
intval perturbation =
1.0e-015 *
0.______________
```

That shows that the function g(x) := f(x)-e with e included in the computed perturbation has a double root contained in Y(1).

## Accuracy of the gamma function I

Usually, Matlab standard functions produce very accurate approximations. For example,

```x = pi
y = gamma(pi)
Y = gamma(intval('pi'))
```
```x =
3.141592653589793
y =
2.288037795340032
intval Y =
2.28803779534003
```

Here Y is a true inclusion of the value of the gamma function at the transcendental number pi: X=intval('pi') is an inclusion of the mathematical pi, and the gamma function for an interval argument produces an inclusion of all values, in particular for the transcendental number pi.

## Accuracy of the gamma function II

For negative arguments something strange happens in Matlab. Consider

```x1 = -171.5;
x2 = -172.5;
y1 = gamma(x1)
y2 = gamma(x2)
```
```y1 =
Inf
y2 =
Inf
```

Regardless of the fact that for larger negative argument the gamma function is very small in absolute value, this cannot be true because the gamma function changes the sign when passing a negative integer:

``` close
axis([-5 5 -10 10])
hold on
for k=-5:5
x = linspace(k,k+1,10000);
plot(x,gamma(x))
end
plot([-5 5],[0 0],':',[0 0],[-10 10],'b:')
hold off
``` The true values of the gamma function at x1=-171.5 and x2=-172.5 are indeed very small:

```Y1 = gamma(intval(x1))
Y2 = gamma(intval(x2))
```
```intval Y1 =
1.0e-309 *
0.19316265431712
intval Y2 =
1.0e-311 *
-0.1119783503283_
```

## Accuracy of the gamma function III

In previous releases of Matlab, the accuracy of the gamma function near negative integers was very poor, sometimes with an error of 100 per cent. That was fixed from Matlab release 2018b on.

However, other problems with the gamma function remain. Consider

```x = -171.25
y = gamma(x)
Y = gamma(intval(x))
```
```x =
-1.712500000000000e+02
y =
Inf
intval Y =
1.0e-309 *
0.98910291950447
```

Obviously gamma(x) must be very small in absolute value. In contrast, the Matlab result is infinity. Moreover,

```x = -172.25
y = gamma(x)
Y = gamma(intval(x))
```
```x =
-1.722500000000000e+02
y =
Inf
intval Y =
1.0e-311 *
-0.5742252072588_
```

still produces +infinity in Matlab, whereas the sign must change from each integer interval to the next.

## Inverse gamma function

Using automatic differentiation, values of the inverse gamma function are easily approximated. For example, compute x such that gamma(x)=3:

```x = 4;
for i=1:5
y = gamma(gradientinit(x)) - 3;
x = x - y.x/y.dx
end
```
```x =
3.601948119538654
x =
3.430611982891113
x =
3.406291166087256
x =
3.405870109547128
x =
3.405869986309577
```

This value is approximate, the true value can be verified as follows:

```X = verifynlss(@(x) gamma(x)-3 , 4)
```
```intval X =
3.40586998630956
```

## Minimum of the gamma function

Similarly, the minimum on the real axis of the gamma function is enclosed by

```mu = verifylocalmin(@gamma,2)
```
```intval mu =
1.46163214496836
```

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