# DEMOINTVAL Interval computations in INTLAB

## Contents

- How to define an interval I
- How to define an interval II
- How to define an interval III
- How to define an interval IV
- How to define an interval V
- Output formats of intervals I
- Output formats of intervals II
- Rigorous output
- What you see is correct
- Output formats of intervals III
- Changing interval output permanently
- Display with uncertainty
- Newton iteration
- Bit representation
- Invoking interval operations
- Interval matrix operations
- Sharp interval multiplication
- Fast interval multiplication
- Acceleration by vector/matrix notation
- Overestimation of interval operations
- Interval standard functions
- Complex interval standard functions
- Standard functions with argument out of range
- A common misuse of interval arithmetic
- Rigorous solution of linear systems
- Accuracy of rigorous linear system solving: Hilbert matrices
- Extremely ill-conditioned linear systems
- Verification of regularity of extremely ill-condtioned matrices
- Structured linear systems
- Sparse linear systems
- Inclusion of eigenvalues and eigenvectors
- Eigenvalue pairs and invariant subspaces
- Eigenvalues of structured matrices
- Nonlinear systems of equations, polynomials, etc.
- Enjoy INTLAB

A key to interval operations in INTLAB is changing the rounding mode. Following we ensure "rounding to nearest".

setround(0)

## How to define an interval I

There are four possibilities to generate an interval, the first is a simple typecast of a real or complex quantity, for example a matrix. It uses Matlab conversion, i.e. the first component does **not** contain "2.3". This is because Matlab first converts "2.3" into binary format, and then the type cast is called.

format compact long infsup u = intval( [ 2.3 -4e1 ; 3 0 ] )

intval u = [ 2.29999999999999, 2.30000000000000] [ -40.00000000000000, -40.00000000000000] [ 3.00000000000000, 3.00000000000000] [ 0.00000000000000, 0.00000000000000]

## How to define an interval II

The second possibility is to use INTLAB conversion of constants. In this case the argument is a string and INTLAB verified conversion to binary is called rather than Matlab conversion.

```
u = intval( '0.1 -3.1 5e2 .3e1' )
```

intval u = 1.0e+002 * [ 0.00099999999998, 0.00100000000001] [ -0.03100000000001, -0.03099999999998] [ 5.00000000000000, 5.00000000000000] [ 0.03000000000000, 0.03000000000000]

The first component, for example, definitely contains "0.1". Since 0.1 is not finitely representable in binary format, the radius of the first component must be nonzero.

u.rad

ans = 1.0e-15 * 0.013877787807814 0.444089209850063 0 0

Generating an interval by an input string always produces a column vector. To change "u" into a 2 x 2 matrix, use

reshape(u,2,2)

intval ans = 1.0e+002 * [ 0.00099999999998, 0.00100000000001] [ 5.00000000000000, 5.00000000000000] [ -0.03100000000001, -0.03099999999998] [ 0.03000000000000, 0.03000000000000]

Note that arrays are stored columnwise in Matlab.

## How to define an interval III

The third possibility is by specification of midpoint and radius.

u = midrad( [ -3e1+2i ; .4711 ; 3 ] , 1e-4 )

intval u = [ -30.00010000000001 + 1.99989999999999i, -29.99989999999999 + 2.00010000000001i] [ 0.47099999999998 - 0.00010000000001i, 0.47120000000001 + 0.00010000000001i] [ 2.99989999999999 - 0.00010000000001i, 3.00010000000001 + 0.00010000000001i]

## How to define an interval IV

The fourth possibility is by specification of infimum and supremum

u = infsup( [ 1 2 ] , [ 4 7 ] )

intval u = [ 1.00000000000000, 4.00000000000000] [ 2.00000000000000, 7.00000000000000]

## How to define an interval V

Yet another possibility is to use [infimum,supremum] or midpoint,radius or giving significant digits in a string:

```
u = intval('[3,4] <5,.1> 1.23_')
```

intval u = [ 3.00000000000000, 4.00000000000000] [ 4.89999999999998, 5.10000000000001] [ 1.21999999999999, 1.24000000000001]

## Output formats of intervals I

The output format can be changed using the different Matlab formats, for example

format midrad long e X = midrad( [ -3e1+2i ; .4711 ; 3 ] , 1e-4 )

intval X = < -3.000000000000000e+01 + 2.000000000000000e+00i, 1.000000000000001e-004> < 4.711000000000000e-01 + 0.000000000000000e+00i, 1.000000000001111e-004> < 3.000000000000000e+00 + 0.000000000000000e+00i, 1.000000000000001e-004>

or

format infsup long X

intval X = [ -30.00010000000001 + 1.99989999999999i, -29.99989999999999 + 2.00010000000001i] [ 0.47099999999998 - 0.00010000000001i, 0.47120000000001 + 0.00010000000001i] [ 2.99989999999999 - 0.00010000000001i, 3.00010000000001 + 0.00010000000001i]

or with a common exponent

1e4*X

intval ans = 1.0e+005 * [ -3.00001000000001 + 0.19998999999999i, -2.99998999999998 + 0.20001000000001i] [ 0.04709999999999 - 0.00001000000001i, 0.04712000000001 + 0.00001000000001i] [ 0.29998999999999 - 0.00001000000001i, 0.30001000000001 + 0.00001000000001i]

## Output formats of intervals II

With two arguments the functions infsup and midrad define an interval, with one argument they control the output of an interval:

```
format short
u = infsup( [ 1 2 ] , [ 4 7 ] );
infsup(u), midrad(u)
```

intval u = [ 1.0000, 4.0000] [ 2.0000, 7.0000] intval u = < 2.5000, 1.5000> < 4.5000, 2.5000>

## Rigorous output

Note that output in INTLAB is rigorous. That means,

left <= ans <= right for inf/sup notation ans in mid+/-rad for mid/rad notation

where ans is the true (real or complex) answer, and left,right, mid,rad are the numbers corresponding to the *displayed* decimal figures.

## What you see is correct

Consider

format short midrad x = midrad(-5,1e-20)

intval x = < -5.0000, 0.0001>

One might get the impression that this is a very crude display. However, the radius of the interval is nonzero. Hence, in four decimal places, the displayed radius 0.0001 is the best possible. Changing the display format gives more information:

```
format long
x
```

intval x = < -5.00000000000000, 0.00000000000001>

## Output formats of intervals III

A convenient way to display narrow intervals is the following:

x=midrad(pi,1e-8); format short, infsup(x), midrad(x), disp_(x) format long, infsup(x), midrad(x), disp_(x) format short

intval x = [ 3.1415, 3.1416] intval x = < 3.1416, 0.0001> intval x = 3.1415 intval x = [ 3.14159264358979, 3.14159266358980] intval x = < 3.14159265358979, 0.00000001000001> intval x = 3.1415927_______

Mathematically the following is true: Form an interval of the displayed midpoint and a radius of 1 unit of the last displayed decimal figure, then this is a correct inclusion of the stored interval.

## Changing interval output permanently

The interval output format can be changed permanently, for example, to infimum/supremum notation:

```
u = midrad( [ -3e1+2i ; .4711 ; 3 ] , 1e-4 );
format infsup
u
```

intval u = [ -30.0002 + 1.9998i, -29.9998 + 2.0002i] [ 0.4709 - 0.0002i, 0.4713 + 0.0002i] [ 2.9998 - 0.0002i, 3.0002 + 0.0002i]

or to midpoint/radius notation:

```
format midrad
u
```

intval u = < -30.0000 + 2.0000i, 0.0002> < 0.4711 + 0.0000i, 0.0002> < 3.0000 + 0.0000i, 0.0002>

or to display with uncertainties depicted by "_":

```
format _
u
```

intval u = -30.0000 + 2.0000i 0.4711 + 0.000_i 3.0000 + 0.000_i

## Display with uncertainty

Display with uncertainty makes only sense for sufficiently narrow intervals. If the accuracy becomes too poor, INTLAB changes automatically to inf-sup or mid-rad display for real or complex intervals, respectively:

for k=-5:-1 disp_(midrad(pi,10^k)) end

intval = 3.1416 intval = 3.142_ intval = 3.14__ intval = 3.1___ [ 3.0415, 3.2416]

## Newton iteration

The following code is an interval Newton iteration to include sqrt(2).

format long _ f = @(x)(x^2-2); % Function f(x) = x^2-2 X = infsup(1.4,1.7); % Starting interval for i=1:4 xs = X.mid; % Midpoint of current interval Y = f(gradientinit(X)); % Gradient evaluation of f(X) X = intersect( X , xs-f(intval(xs))/Y.dx ) % Interval Newton step with intersection end

intval X = 1.4_____________ intval X = 1.4142__________ intval X = 1.414213562_____ intval X = 1.41421356237309

The "format _" output format shows nicely the quadratic convergence.

The last displayed result (which is in fact an interval) proves that the true value of sqrt(2) is between 1.41421356237308 and 1.41421356237310. Indeed, sqrt(2)=1.41421356237309504...

The format "long e" in Matlab displays the most figures. With this we see that the internal accuracy of the final X is in fact even better, the width is only 2 units in the last place.

format long e X format short

intval X = 1.414213562373095e+000

## Bit representation

The bit pattern of floating-point numbers and of intervals as well is displayed as follows.

getbits(X)

ans = 4×62 char array 'infimum ' ' +1.0110101000001001111001100110011111110011101111001011 * 2^0' 'supremum ' ' +1.0110101000001001111001100110011111110011101111001101 * 2^0'

That shows indeed that the left and right bound of X coincide up to 2 bits.

## Invoking interval operations

An operation uses interval arithmetic if at least one of the operands is of type intval. For example, in

u = intval(5); y = 3*u-exp(u)

intval y = -133.4131

the result y is an inclusion of 15-exp(5). However, in

u = intval(5); z = 3*pi*u-exp(u)

intval z = -101.2892

the first multiplication "3*pi" is a floating point multiplication. Thus it is not guaranteed that the result z is an inclusion of 15pi-exp(5).

The following ensures a correct inclusion of the mathematical 3 \pi u - exp(u):

Z = 3*intval('pi')*u - exp(u) format long infsup(Z) format short

intval Z = -101.2892 intval Z = 1.0e+002 * [ -1.01289269298730, -1.01289269298729]

## Interval matrix operations

INTLAB is designed to be fast. Case distinctions in interval multiplication can slow down computations significantly due to interpretation overhead. Therefore, there is a choice between 'fast' and 'sharp' evaluation of interval matrix products. This applies only to 'thick' intervals, i.e. intervals with nonzero diameter.

## Sharp interval multiplication

In the following example, c is a real random matrix, C is an interval matrix with diameter zero (a thin interval matrix), and CC is an interval matrix with nonzero diameter (a thick interval matrix), all of dimension nxn for n=1000. First we measure the computing time with option 'SharpIVmult'.

```
n = 1000;
c = randn(n);
C = intval(c);
C_ = midrad(c,.1);
intvalinit('SharpIVmult')
tic, scc = c*c; toc
tic, sCC = C*C; toc
tic, sCC = C*C_; toc
tic, sCC__ = C_*C_; toc
```

===> Slow but sharp interval matrix multiplication in use Elapsed time is 0.053658 seconds. Elapsed time is 0.116197 seconds. Elapsed time is 0.177183 seconds. Elapsed time is 14.014845 seconds.

As can be seen, there is not much penalty if not both matrices are thick interval matrices; for both thick intervals, however, computation is slowed down significantly.

## Fast interval multiplication

```
intvalinit('FastIVmult')
tic, fcc = c*c; toc
tic, fCC = C*C; toc
tic, fCC = C*C_; toc
tic, fCC__ = C_*C_; toc
max(max(diam(fCC__)./diam(sCC__)))
```

===> Fast interval matrix multiplication in use (maximum overestimation factor 1.5 in radius, see FAQ)Elapsed time is 0.042636 seconds. Elapsed time is 0.106939 seconds. Elapsed time is 0.162284 seconds. Elapsed time is 0.212915 seconds. ans = 1.0615

As can be seen there is again not much penalty if not both matrices are thick. However, the 'fast' implementation is much faster than the 'sharp' at the cost of a little wider output. If intervals are very wide and any overestimation cannot be afforded (as in global optimization), the option 'SharpIVmult' may be used. It is shown in

S.M. Rump: Fast and parallel interval arithmetic. BIT Numerical Mathematics, 39(3):539-560, 1999

that the maximum (componentwise) overestimation by the option 'FastIVmult' compared to 'SharpIVmult' is a factor 1.5, for real and complex intervals. For not too thick intervals, however, it is shown that the overestimation is negligible.

## Acceleration by vector/matrix notation

It is advisable to use vector/matrix notation when using interval operations. Consider

n = 10000; x = 1:n; y = intval(x); tic for i=1:n y(i) = y(i)^2 - y(i); end t1 = toc

t1 = 4.7558

This simple initialization takes considerable computing time. Compare to

tic y = intval(x); y = y.^2 - y; t2 = toc ratio = t1/t2

t2 = 0.0048 ratio = 991.7781

Sometimes code looks more complicated, a comment may help. It is worth it.

## Overestimation of interval operations

Note that the main principle of interval arithmetic is that for given intervals A,B and an operation o, the result a o b is included in the interval result A o B for all a in A and all b in B. Since the result must be an interval, overestimations cannot be avoided in many situations. For example, in

close, kmax = 40; i = sqrt(-1); a=midrad(2,.7); b=midrad(1-i,1); plotintval(3*a-exp(b)); hold on phi = linspace(0,2*pi,kmax); [A,B] = meshgrid( mid(a)+rad(a)*exp(i*phi) , mid(b)+rad(b)*exp(i*phi) ); plot(3*A-exp(B)) hold off

the blue circle is the result of the interval operations, whereas the many circles approximate the power set operation (see also the DEMOINTLAB). Another reason for overestimation are dependencies, see below.

## Interval standard functions

Interval standard functions in INTLAB are rigorous. For a given interval X and a function X let Y be the computed value of f(X). Then f(x) is in Y for all x in X. For example

```
x = intval(1e10); format long
sin(x)
```

intval ans = -0.48750602508751

Note that the result is rigorous (try sin(2^1000) or similar). For timing comparison consider

```
format short
n=10000; x=random(n,1); X=intval(x);
tic, exp(x); tapprox = toc
tic, exp(X); trigorous = toc
ratio = trigorous/tapprox
```

tapprox = 0.0025 trigorous = 0.0061 ratio = 2.4648

## Complex interval standard functions

Complex interval standard functions are rigorous as well, for example

```
format long
Z = midrad(3+4i,1e-7);
R = sin(Z)
```

intval R = 3.85374_________ - 27.01681_________i

It is mathematically correct, that sin(z) is an element of R for every complex z with Eucledian distance less than or equal to 1e-7 from 3+4i.

## Standard functions with argument out of range

When entering a real argument leading to a complex value of a standard function, there are three possibilities to be specified by intvalinit:')

intvalinit('RealStdFctsExcptnNan'); sqrt(intval(-2)) intvalinit('RealStdFctsExcptnWarn'); sqrt(intval(-2)) intvalinit('RealStdFctsExcptnAuto'); sqrt(intval(-2))

===> Result NaN for real interval input out of range intval ans = NaN ===> Complex interval stdfct used automatically for real interval input out of range, but with warningintval ans = -0.00000000000000 + 1.41421356237309i ===> Complex interval stdfct used automatically for real interval input out of range (without warning)intval ans = -0.00000000000000 + 1.41421356237309i

## A common misuse of interval arithmetic

The dependency problem is the most serious problem of (naive) interval arithmetic. The following procedure:

" Take some numerical algorithm and replace every operation by its corresponding interval operation. Then the computed interval result(s) definitely contain the true result which would be obtained without the presence of rounding errors. "

will most certainly fail in practice. Although a true statement (if no exception like divide by a zero interval occurs), the computed result interval(s) will, for very modest problem size, most certainly be of huge diameter and thus useless.

Consider, for example, the triangular matrix T where all elements on and below the diagonal are equal to 1, and take a randomly generated right hand side. The following lines do this for dimension n=50:

n = 50; T = tril(ones(n)); b = randn(n,1);

Then perform a standard forward substitution to compute an inclusion T\b. Note that X is defined to be an interval vector, so all operations are indeed interval operations (see above section "Invoking interval operations").

X = intval(zeros(n,1)); for i=1:n X(i) = b(i) - T(i,1:i-1)*X(1:i-1); end X

intval X = 1.54496886565095 -1.73641760729346 0.92018854401184 0.65197328202597 -2.33827755592769 0.16629614976347 2.43421718014301 -1.1757752311920_ -0.3244045981527_ -0.3109015016353_ 0.8243081711104_ 0.062040815522__ -0.257822709060__ -1.33661157686___ 1.18458604457___ -0.33469668576___ 0.11636807926___ -0.1864900273____ 1.0265758952____ -1.1168262990____ 0.3330179902____ 0.397103103_____ -2.959729555_____ 3.799835771_____ -0.87738401______ -2.15995104______ 2.00125210______ -0.7093691_______ 0.1979432_______ 0.1816328_______ -0.793923________ 1.278688________ -1.115279________ 0.108663________ 2.59745_________ -1.95887_________ 1.58669_________ 0.7210__________ -1.4453__________ -1.4330__________ 0.951___________ -1.588___________ 1.736___________ -2.084___________ 1.10____________ 0.06____________ -1.63____________ 1.2_____________ 0.6_____________ 0.4_____________

The result is displayed with uncertainty, perfectly making visible the exponential loss of accuracy. This is due to one of the most common misuses of interval arithmetic, also called "naive interval arithmetic". For more details and examples cf.

S.M. Rump: Verification methods: Rigorous results using floating-point arithmetic. Acta Numerica, 19:287-449, 2010.

to be downloaded from "www.ti3.tuhh.de/rump". See, in particular,

Chapter 10.1: The failure of the naive approach: interval Gaussian elimination (IGA)

Note that the linear system above is very well-conditioned:

cond(T)

ans = 64.270085531579554

By the well-known rule of thumb of numerical analysis we expect at least 14 correct digits in a floating-point approximation T\b. Using a proper (non-naive) method, an inclusion of this quality is indeed achieved:

verifylss(T,b)

intval ans = 1.54496886565095 -1.73641760729346 0.92018854401184 0.65197328202597 -2.33827755592769 0.16629614976346 2.43421718014301 -1.17577523119201 -0.32440459815271 -0.31090150163526 0.82430817111040 0.06204081552205 -0.25782270906023 -1.33661157686160 1.18458604456945 -0.33469668576355 0.11636807926340 -0.18649002731532 1.02657589515373 -1.11682629899407 0.33301799020482 0.39710310294512 -2.95972955511010 3.79983577121854 -0.87738400810208 -2.15995104006627 2.00125209821602 -0.70936905983630 0.19794316597518 0.18163284792087 -0.79392266302791 1.27868838309226 -1.11527942977767 0.10866280278142 2.59745200007108 -1.95886795922010 1.58668718772608 0.72102363349527 -1.44532990095592 -1.43303475689456 0.95128027871488 -1.58813522029171 1.73628335123737 -2.08378329283865 1.10008919761228 0.05734129838531 -1.62759197715764 1.24373722971505 0.59753040062128 0.38015759066768

Such methods are called "self-validating methods" or "verification methods". For some details see the reference above or

S.M. Rump: Self-validating methods. Linear Algebra and its Applications (LAA), 324:3-13, 2001.

Due to an improved evaluation of the residual (default option "intvalinit('ImprovedResidual')" , see also function "lssresidual.m") 15 correct decimal digits of the result are computed.

## Rigorous solution of linear systems

The INTLAB linear system solver can be called with "\" or "verifylss".' For example, [bare with me, I am often in Japan where the backslash appears like japanese Yen.]

n = 100; A = randn(n); b = A*ones(n,1); X = verifylss(A,b);

generates and solves a randomly generated 100x100 linear system. The inclusion and its quality is checked by

X([1:3 98:100]) max( X.rad ./ abs(X.mid) )

intval ans = 1.00000000000000 0.99999999999999 0.99999999999999 0.99999999999999 1.00000000000000 0.99999999999999 ans = 5.551115123125781e-16

which calculates the maximum relative error of the inclusion radius with respect to the midpoint. The same is done by

max(relerr(X))

ans = 4.440892098500623e-16

## Accuracy of rigorous linear system solving: Hilbert matrices

For estimating accuracy, try

format long e n = 10; H = hilb(n); b = ones(n,1); X = verifylss(H,b)

intval X = -9.99830__________e+000 9.89853__________e+002 -2.375688_________e+004 2.402116_________e+005 -1.261125_________e+006 3.783408_________e+006 -6.726110_________e+006 7.000691_________e+006 -3.937911_________e+006 9.237120_________e+005

The notoriously ill-conditioned Hilbert matrix is given by H_ij := 1/(i+j-1). To estimate the accuracy, we use the symbolic toolbox to compute the perturbation of the solution when perturbing only the (7,7)-element of the input matrix by 2^(-52):

```
Hs = sym(H,'f');
Hs(7,7) = Hs(7,7)*(1+sym(2^(-52)));
double( Hs \ b )
```

ans = -9.997198050336207e+00 9.897576960772726e+02 -2.375483680969619e+04 2.401930526008937e+05 -1.261036061017173e+06 3.783164425266260e+06 -6.725710140930751e+06 7.000304229698808e+06 -3.937707813532462e+06 9.236673822756851e+05

The statement "sym(H,'f')" makes sure that no conversion error occurs when changing H into symbolic format.

## Extremely ill-conditioned linear systems

By default, all computations in INTLAB are, like in Matlab, performed in double precision. This restricts treatable linear systems to a maximum condition number of roughly below 10^16.

Starting with INTLAB Version 7, I rewrote my linear system solver completely. Now, although only double precision is used, linear systems with larger condition numbers are solvable. Consider

format long _ n = 20; A = invhilb(n); condA = norm(double(inv(sym(A))))*norm(A)

condA = 5.478702657981728e+27

The common rule of thumb tells that for a condition number 10^k, an algorithm in double precision should produce 16-k correct digits. In our case this means roughly 16-27=-11 "correct" digits, namely none. For a random right hand side Matlab computes

b = A*randn(n,1); x = A\b

x = 1.0e+05 * -3.526561759042908 -1.454567045706571 -0.592520417296815 -0.235078095024219 -0.089144109711128 -0.031474527357279 -0.009882856050908 -0.002503614144822 -0.000398423032250 0.000031205126401 0.000038665477879 0.000006113262176 -0.000029375253182 0.000009822814340 0.000032598385605 0.000010223756347 0.000008182354189 -0.000003888579149 0.000057657456205 0.000191873032237

A corresponding warning indicates the difficulty of the problem. Note that in this case the Matlab guess of the condition number is pretty good.

A verified inclusion of the solution is computed by

```
X = verifylss(A,b,'illco')
```

intval X = 1.0e+005 * -3.79143_________ -1.65610894______ -0.709759661_____ -0.2950853244____ -0.1169928675649_ -0.04311353750180 -0.01409374776675 -0.00368058769965 -0.00054223078157 0.00011892565814 0.00010397698909 0.00001270523760 -0.00005403043476 -0.00001774136504 0.00001222576888 -0.00001114755428 -0.00003401120403 -0.00009069689194 -0.00009515901791 -0.00004152613564

As expected the Matlab approximation differs significantly from the true values, for some components the sign is incorrect. The maximum relative error of the components of the computed inclusion is

max(relerr(X))

ans = 3.511080212154039e-07

so that each component is accurately included with at least 7 correct figures.

## Verification of regularity of extremely ill-condtioned matrices

For accurate dot product computations some reference implementations are used which suffer from interpretation overhead. This makes it time consuming to solve linear systems with extremely ill-conditioned matrix.

However, checking regularity of very extremely ill-conditioned matrices can be done very fast using mod-p arithmetic. Consider

```
n = 50;
cnd = 1e200;
digits(400)
A = randmat(n,cnd);
tic, Cnd = norm(A,inf)*norm(double(inv(sym(A,'f'))),inf), toc
```

Cnd = 7.086324698036915e+226 Elapsed time is 4.833116 seconds.

As can be seen the matrix A is indeed extremely ill-conditioned. However, checking regularity in mod-p arithmetic is rigorous and very fast:

tic, isregular(A), toc

ans = 1 Elapsed time is 0.010447 seconds.

With a chance of about 1/p regularity cannot be checked. Since p is of the order 1e8, this chance is slim. It can be further diminished by using several primes, see "isregular".

Note that verification of regularity, i.e. a sufficient criterion, is simple, whereas a necessary and sufficient check of singularity is much more involved.

## Structured linear systems

In general, intervals are treated as independent quantities. If, however, there are dependencies, then taking them into account may shrink the solution set significantly. An example is

format short n = 4; e = 1e-3; intvalinit('displayinfsup'); A = midrad( toeplitz([0 1+e e 1+e]),1e-4); b = A.mid*ones(n,1); Amid = A.mid X1 = verifylss(A,b)

===> Default display of intervals by infimum/supremum (e.g. [ 3.14 , 3.15 ]) Amid = 0 1.0010 0.0010 1.0010 1.0010 0 1.0010 0.0010 0.0010 1.0010 0 1.0010 1.0010 0.0010 1.0010 0 intval X1 = [ 0.3327, 1.6673] [ 0.3327, 1.6673] [ 0.3327, 1.6673] [ 0.3327, 1.6673]

First the matrix has been treated as a general interval matrix without dependencies. Recall that only the midpoint is displayed above; all entries of the interval matrix have a uniform tolerance of 1e-4.

Several structural information may be applied to the input matrix, for example,

X2 = verifystructlss(structure(A,'symmetric'),b); X3 = verifystructlss(structure(A,'symmetricToeplitz'),b); X4 = verifystructlss(structure(A,'generalToeplitz'),b); X5 = verifystructlss(structure(A,'persymmetric'),b); X6 = verifystructlss(structure(A,'circulant'),b); res = [ X1 X2 X3 X4 X5 X6 ]; rad(res)

ans = 0.6672 0.5062 0.1689 0.3376 0.5062 0.0003 0.6672 0.5062 0.1688 0.3375 0.5062 0.0003 0.6672 0.5062 0.1688 0.3375 0.5062 0.0003 0.6672 0.5062 0.1689 0.3376 0.5062 0.0003

Here only the radii of the inclusions are displayed. Note that the inclusion may become much narrower, in particular treating the input data as a circulant matrix.

## Sparse linear systems

The following generates a random sparse system with about 9 nonzero elements per row.

```
format short
n = 10000; A = sprand(n,n,2/n)+speye(n); A = A*A'; b = ones(n,1);
```

The linear system is generated to be symmetric positive definite. Before calling the verified linear system solver, the fill-in should reduced. The original matrix looks like

```
spy(A)
title('sparsity pattern of A')
```

whereas after minimum degree reordering the matrix looks like

```
p = symamd(A);
spy(A(p,p))
title('sparsity pattern of renumbered A')
```

The timing for the built-in (floating point) solver compared to the verified solver is as follows:

tic, x = A(p,p)\b(p); toc

Elapsed time is 0.525967 seconds.

tic, X = verifylss(A(p,p),b(p)); toc

Elapsed time is 2.297256 seconds.

## Inclusion of eigenvalues and eigenvectors

To compute verified inclusions of eigenvalue/eigenvector pairs of simple or multiple eigenvalues, consider, for example, the famous Wilkinson(21) matrix. Following inclusions of the last four eigenvalue/eigenvector pairs are displayed. Those are the most ill-conditioned.

format long _ W = wilkinson(21); % generation of the matrix [V,D] = eig(W); % eigenvalue/eigenvector approximations for k=18:21 [L,X] = verifyeig(W,D(k,k),V(:,k)) % inclusions for the small eigenvalues end

intval L = 9.21067864730491 intval X = -0.38247_________ 0.30189_________ 0.44607_________ 0.23816_________ 0.080419________ 0.020042________ 0.0039719_______ 0.00065440______ 0.000092313_____ 0.0000112430____ -0.00000000000___ -0.000011243042__ -0.000092313004__ -0.0006543963623_ -0.00397193251082 -0.02004206756034 -0.08041877341334 -0.23815677108033 -0.44606931512504 -0.30188982395948 0.38246757537813 intval L = 9.21067864736133 intval X = 0.38246757533800 -0.30188982390622 -0.44606931509072 -0.23815677111720 -0.08041877354260 -0.02004206794300 -0.00397193399399 -0.0006544037082_ -0.000092357143__ -0.00001155397___ -0.00000250882___ -0.0000115540____ -0.000092357_____ -0.00065440______ -0.0039719_______ -0.020042________ -0.080419________ -0.23816_________ -0.44607_________ -0.30189_________ 0.38247_________ intval L = 10.74619418290332 intval X = 0.55939864230126 0.41742001280922 0.16949775589362 0.04805373844102 0.01052087952088 0.00188039874002 0.0002842567806_ 0.000037252700__ 0.00000430986___ 0.0000004422____ -0.000000000_____ -0.00000044______ -0.0000043_______ -0.000037________ -0.00028_________ -0.00188_________ -0.0105__________ -0.048___________ -0.170___________ -0.42____________ -0.56____________ intval L = 10.74619418290339 intval X = 0.56____________ 0.42____________ 0.170___________ 0.048___________ 0.0105__________ 0.00188_________ 0.00028_________ 0.000037________ 0.0000043_______ 0.00000045______ 0.000000084_____ 0.0000004509____ 0.00000431088___ 0.000037252833__ 0.0002842568009_ 0.00188039874370 0.01052087952171 0.04805373844126 0.16949775589369 0.41742001280919 0.55939864230118

## Eigenvalue pairs and invariant subspaces

The largest eigenvlues are 10.74619418290332 and 10.74619418290339, where all displayed digits are verified to be correct. Invariant subspaces of nearby eigenvalues are in general ill-conditioned. Nearby eigenvalues can also be regarded as clusters. From the inclusions above we can judge how narrow the eigenvalues are. So one of the approximations can be used as an approximation of the pair.

[L,X] = verifyeig(W,D(18,18),V(:,18:19)) % inclusion of the 18/19 eigenvalue pair [L,X] = verifyeig(W,D(20,20),V(:,20:21)) % inclusion of the 20/21 eigenvalue pair

intval L = 9.2106786473____ + 0.0000000000____i intval X = -0.38246721566493 0.38246757533800 0.30188954003016 -0.30188982390622 0.44606889559399 -0.44606931509072 0.23815654709237 -0.23815677111720 0.08041869777897 -0.08041877354260 0.02004204871064 -0.02004206794300 0.00397192877519 -0.00397193399399 0.00065439574687 -0.00065440370821 0.00009231291679 -0.00009235714338 0.00001124303112 -0.00001155397350 -0.00000000000118 -0.00000250882017 -0.00001124304199 -0.00001155396293 -0.00009231300365 -0.00009235705656 -0.00065439636234 -0.00065440309275 -0.00397193251082 -0.00397193025836 -0.02004206756034 -0.02004204909331 -0.08041877341334 -0.08041869790822 -0.23815677108033 -0.23815654712923 -0.44606931512504 -0.44606889555967 -0.30188982395948 -0.30188953997691 0.38246757537813 0.38246721562480 intval L = 10.7461941829034_ + 0.0000000000000_i intval X = 0.55939864230126 0.53926516857120 0.41742001280922 0.40239653183024 0.16949775589362 0.16339731453128 0.04805373844102 0.04632422283758 0.01052087952089 0.01014221959041 0.00188039874009 0.00181272078417 0.00028425678094 0.00027402603484 0.00003725270198 0.00003591205802 0.00000430988244 0.00000415574012 0.00000044236679 0.00000043485205 0.00000000151024 0.00000008241241 -0.00000042613744 0.00000045076775 -0.00000415472850 0.00000431085765 -0.00003591192480 0.00003725283040 -0.00027402601454 0.00028425680050 -0.00181272078050 0.00188039874363 -0.01014221958959 0.01052087952169 -0.04632422283735 0.04805373844125 -0.16339731453120 0.16949775589369 -0.40239653183027 0.41742001280919 -0.53926516857129 0.55939864230118

Note that interval output with uncertainty ("_") is used, so all displayed decimal places of the bases of the invariant subspaces are verified to be correct. As explained in section "Output formats of intervals III", the inclusion 10.7461941829034_ of the two smallest eigenvlues reads [10.7461941829033,10.7461941829035], thus including the true eigenvalues as displayed above.

The mathematical statement is that the displayed intervals for the cluster contain (at least) two eigenvalues of the Wilkinson matrix W. The size of the cluster is determined by the number of columns of the invariant subspace approximation.

## Eigenvalues of structured matrices

As for linear systems, the interval input matrix may be structured. Taking into account such structure information may shrink the inclusion. As an example consider

format short e = 1e-3; A = midrad( toeplitz([0 1+e -e/2 1+e]),1e-4); [v,d] = eig(A.mid); xs = v(:,2:3); lambda = d(2,2); X1 = verifyeig(A,lambda,xs); X2 = verifystructeig(structure(A,'symmetric'),lambda,xs); X3 = verifystructeig(structure(A,'symmetricToeplitz'),lambda,xs); X4 = verifystructeig(structure(A,'generalToeplitz'),lambda,xs); X5 = verifystructeig(structure(A,'persymmetric'),lambda,xs); X6 = verifystructeig(structure(A,'circulant'),lambda,xs); res = [ X1 X2 X3 X4 X5 X6 ]; rad(res)

ans = 1.0e-03 * 0.7401 0.6455 0.3378 0.4850 0.5700 0.4000

As for linear systems, only the radii of the inclusions are displayed.

## Nonlinear systems of equations, polynomials, etc.

For inclusions of systems of nonlinear equations, of roots of polynomials etc. cf. the corresponding demos.

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