Interval Arithmetic and Solvability of Boundary Value Problems

Jiří Benedikt

Západočeská univerzita v Plzni, Czech Republic

This work has been supported by the Research Plan MSM 4977751301 of the Ministry of Education, Youth and Sports of the Czech Republic.

Machine precision

Let us start with a few examples of problems with computing in machine precision. The following expressions are expected to evaluate to 1:

a=10.^16;
a+1-a
a-a+1
a=10.^16+2;
a+1-a

a=Sqrt[2.];
SetPrecision[a+1,30]-SetPrecision[a,30]

The polynomial (x-1000)(x-1001)(x-1002)(x-1003) has roots 1000, 1001, 1002, 1003:

Expand[(x-1000.)(x-1001.)(x-1002.)(x-1003.)];
x/.Solve[%==0,x]
%-1000

The polynomial (x-10000)(x-10001)(x-10002)(x-10003) has roots 10000, 10001, 10002, 10003:

Expand[(x-10000.)(x-10001.)(x-10002.)(x-10003.)];
x/.Solve[%==0,x]

Expand[(x-10000.)(x-10001.)(x-10002.)(x-10003.)]

If we do not expand the polynomial, we get a correct answer:

(x-10000.)(x-10001.)(x-10002.)(x-10003.)
x/.Solve[%==0,x]

Sine at multiples of π is expected to be zero:

Table[Sin[10.^12*Pi*n],{n,1000000,1000010}]

IEEE 754 standard

The above errors appear because when we compute in machine precision, the operations are performed by CPU (FPU) and their results are rounded to the nearest machine number representable in IEEE 754 standard. In double precision, machine number is a number of magnitude greater or equal to 2^-1022 and less than 2^1024, and with at most 53 significant binary digits. Although these numbers are often called real numbers, it is only a finite selection of special rational numbers. Notice that, e.g., all rational numbers in the interval (0,1) cannot be represented in memory of any computer since it is an infinite set and a computer memory has only a finite number of states.

As for the first example, the numbers 10^16 and 10^16+2 are machine numbers, but 10^16+1 is not, because it has 54 significant binary digits. Hence the sum 10^16+1 must be rounded either down to 10^16 or up to 10^16+2. In both cases, the error is 1. FPU chooses rounding down. Subtracting 10^16, we then get 0.

a=10.^16;
TableForm[Table[{i,BaseForm[10^16+i,2],BaseForm[Round[a+i],2],Round[a+i-a]},{i,0,8}],TableHeadings->{None,{"i","10^16 + i in binary form","nearest machine number in binary form","nearest machine nubmer - 10^16"}}]

The built-in constant $MachineEpsilon is the gap between 1 and the next machine number. It equals 2^-52.

$MachineEpsilon
2.^-52

1+$MachineEpsilon/2 is not a machine number, it is in the middle between 1 and 1+$MachineEpsilon. It evaluates to 1 (the result is rounded down).

SetPrecision[1+$MachineEpsilon,30]
SetPrecision[1+$MachineEpsilon/2,30]

The gap between the machine numbers left to 1 (and right to 1/2) is twice smaller, between 2 and 4 it is twice bigger, between 4 and 8 it doubles again, and so on.

SetPrecision[1-$MachineEpsilon/2,30]

The following graph show the difference between x and the nearest machine number (rounding error) for x between 1-5*10^-16 and 1+5*10^-16.

Plot[(1+x*10^-16)-SetPrecision[SetPrecision[1+x*10^-16,MachinePrecision],30],{x,-5,5},PlotRange->{-1.1*10^-16,1.1*10^-16},WorkingPrecision->30,PlotPoints->1000]

The reason why sine in the example above did not evaluate to 0 is that, e.g., the difference between 10^18π and the nearest machine number is more than 40.

SetPrecision[10.^12*Pi*1000000,30]
N[10^18*Pi,30]
%%-%

Tracking the relative error

We can track the error in Mathematica using the notation a`b where b is minus decadic logarithm of the relative error. For example, 10.`1 is 10 with the relative error 10^-1=0.1, i.e., a number between 9 and 11. The error is tracked during the computation. When we compute the square of 10.`1, the result is 100 with double relative error (0.2), i.e., a number between 80 and 120. But the square of a number between 9 and 11 need not be between 80 and 120 (it can be even 121)! In other words, this error tracking gives only a first-order estimate of the relative error, not the worst case. Hence, it is not guaranteed that the error is not greater.

a=10.`1;
FullForm[a]
10*(1+{-1,1}*10^-1.)
FullForm[a*a]
100*(1+{-1,1}*10^-.698970004336)

Let us now consider the number 1000π with the relative error 10^-3, i.e., a number between 999π and 1001π. Sine of such number can be anywhere between -1 and 1, but we get 0 with the relative error 10^-0=1.

a=SetPrecision[1000Pi,3]
SetPrecision[a*(1+{-1,1}*10^-3),6]
Sin[a]//FullForm

Interval arithmetic

The previous examples are solved using Interval.

Interval[{9,11}]^2
Sin[Interval[{999Pi,1001Pi}]]

If we plug a machine number x into Interval, we get an interval I with machine end-points, such that any real number for which the nearest machine number is x, is contained in I.

a=Interval[1.];
SetPrecision[a,30]
SetPrecision[{1-2*2^-53,1+2^-52},30]

The square of the interval a with machine end-points is again an interval with machine end-points which contains squares of all real numbers in a. Notice that the right end-point of a^2 is the third machine number after 1, because the square of the first machine number after 1 is bigger than the second machine number after 1 (the difference is only 2^-104, but nothing is neglected).

SetPrecision[a^2,30]
SetPrecision[{1-5*2^-53,1+3*2^-52},30]

SetPrecision[1+2*2^-52,40]
SetPrecision[(1+2^-52)^2,40]

The following example shows that the sum of two intervals in Mathematica need not contain all sums of arbitrary two numbers from these intervals. For example, if the right end-point of b is positive, but so small that the sum of the right end-point of a and the right end-point of b rounds down back to the right end-point of a, then the right end-point of a+b is equal to the right end-point of a. But the sum of the right end-points is greater. It violates the expectation that the result of an operation with intervals is an interval containing all possible results of the operation with arbitrary numbers from these intervals.

a=Interval[1.];
SetPrecision[a,30]
b=Interval[{0.,1.1*10^-16}]
SetPrecision[a+b,30]

b=Interval[{0.,1.2*10^-16}]
SetPrecision[a+b,30]

SetPrecision[Max[a]+1.1*10^-16,30]
SetPrecision[Max[a]+1.2*10^-16,30]

Interval Riemann integral and boundary value problems

Integrate[Sin[x],{x,0,Pi}]

Using interval arithmetic, we can easily compute lower and upper Riemann sum on a given partition of the interval. We get an interval that certainly contains the Riemann integral, if it exists.

left=0.;
right=N[Pi];
f[x_]=Sin[x];
n=100;
Sum[f[Interval[{left+i*(right-left)/n,left+(i+1)*(right-left)/n}]]*(right-left)/n,{i,0,n-1}]

We can even estimate integrals of functions that cannot be handled by Integrate or NIntegrate.

Integrate[Sin[x*Interval[{0.99,1.01}]],{x,0,Pi}]
NIntegrate[Sin[x*Interval[{0.99,1.01}]],{x,0,Pi}]

Integrate[Min[Sin[.99*x],Sin[1.01*x]],{x,0,Pi}]

left=0.;
right=N[Pi];
f[x_]=Sin[x*Interval[{0.99,1.01}]];
n=10000;
Sum[f[Interval[{left+i*(right-left)/n,left+(i+1)*(right-left)/n}]]*(right-left)/n,{i,0,n-1}]

We finish by showing an application to solvability of a boundary value problem. The well-known Fredholm alternative implies that the boundary value problem for the second-order linear ODE u''(x)+u(x)=f(x) where f is continuous on [0,π], and the Dirichlet boundary conditions u(0)=u(π)=0, has a classical solution if and only if (∫_0^π)f(x)sin(x)dx=0. Let us consider f(x)=sin((x^2))-1/4.

Integrate[Sin[x]*(Sin[x^2]-1/4),{x,0,Pi}]
NIntegrate[Sin[x]*(Sin[x^2]-1/4),{x,0,Pi}]

We cannot be sure if the result is zero or not, because Integrate gives a very complicated expression and NIntegrate gives a number near zero and we do not have a guaranteed upper estimate of the error.

left=0.;
right=N[Pi];
f[x_]=Sin[x]*(Sin[x^2]-1/4);
n=10000;
Sum[f[Interval[{left+i*(right-left)/n,left+(i+1)*(right-left)/n}]]*(right-left)/n,{i,0,n-1}]

Using Interval, we get that the integral lies in an interval that does not contain zero. Hence it is certainly nonzero and the boundary value problem is not solvable. However, NDSolve solves it with no warning:

s[x_]=u[x]/.NDSolve[{u''[x]+u[x]==Sin[x^2]-1/4,u[0]==0,u[Pi]==0},u[x],{x,0,Pi}][[1]]

Plot[s[x],{x,0,Pi}]

Apparently, the "solution" resembles a big multiple of sine. Any multiple of sin(x) is a solution of the homogeneous Dirichlet problem u''(x)+u(x)=0, u(0)=u(π)=0. Since the right-hand side of our nonhomogeneous equation is very small comparing to u(x) and u''(x), it vanishes in the rounding errors. This is a nice example of using a numerical computation to prove that another numerical computation (both in machine precision) is meaningless.