Mathematical topics covered:
Programming concepts covered:
mpmathNote: This material is not covered in the various textbooks.
If $A$ is a set of real numbers, we write $a \in A$ to denote the statement "$a$ is in $A$".
In mathematics, we can extend the usual operations on real numbers to sets of real numbers.
For example, the addition operation is $$A+B = \{a+b:~\text{$a \in A$ and $b \in B$}\}.$$ So, for instance $$\{1,3\}+\{5,7\}= \{6,8,10\}.$$
Basically interval arithmetic implements these operations for intervals. Most operations between intervals yield new intervals. For example, $$[a,b] + [c,d] = [a+c, b+d].$$
The one operation that is an issue is division. Since for example, $$[1,2] \div [-1,1]$$ is not an interval. But then one of the numbers in the resulting set should be $1 \div 0$ which is undefined. So, when implemented as a computer program, this operation results in an error instead of a strange set.
Recall that most numbers you want to work with can't be stored exactly in Python. When we do arithmetic, we round at every step. As we do more operations, the value stored slowly diverges from the result of the true computation.
Interval arithmetic gives a partial resolution to this. If $x$ is a real number, we use a (hopefully small) interval containing $x$ to represent $x$. When we do computations involving $x$, we round outward (making the intervals bigger) to be sure that the real value is somewhere inside. This guarantees that the true value that would result from the computation must be inside the interval.
We will demonstrate the Arbitrary-precision interval arithmetic built into the mpmath package. (Warning: The package states that interval arithmetic support is still experimental.)
The inteval arithmetic package can be loaded with:
from mpmath import iv
We will work with $5$ digits of precision just so we can see what is going on easily. You would normally want to set this higher.
iv.dps = 5 # Use 5 digits of precision
The following repesents the interval $[1, 1]$. That is $1$ is the only point in the interval.
x = iv.mpf(1)
x
Recall that $1/3$ can not be stored as a floating point number. Here is $1/3$ as an interval:
y = x/3
y
Observe that $1/3$ is in the interval above. Also if we add $y$ three times, we get an interval containing one:
y + y + y
The number $\pi$ can't be stored exactly and so is represented as an interval.
pi = iv.pi
pi
The standard functions are included and are applied to intervals with rounding outwards. For example, computing $cos(\frac{\pi}{2})$ results in a interval containing $0$.
iv.cos(pi/2)
Observe that this rigorously guarantees that the true value has absolute value less than $1.3 \times 10^{-6}$. This isn't the best example, because we know the answer is zero. But we don't know the true value of $\cos(1)$.
iv.cos(1)
The above guarantees that $$\cos(1)=0.54030\ldots$$ is correct to the decimal places displayed (and the next digit is either a $2$ or a $3$).
The following command constructs an interval with endpoints $2$ and $3$.
iv.mpf([2,3])
Intervals are only equal if they have the same endpoints. (Overlapping is not good enough.) This does create some issues. For example $$\frac{2}{7}+\frac{1}{5} = \frac{17}{35}$$ but:
x = iv.mpf("2/7") + iv.mpf("1/5") # 2/7 + 1/5
print(x)
y = iv.mpf("17/35")
print(y)
x==y
(Above we made use of conversion from string to an interval.)
Less than and greater than comparisons return True, False, or None. Here None is reserved for the case that the answer can not be determined (usually because the intervals overlap.)
print( iv.mpf([2,3]) < iv.mpf([4,5]) )
print( iv.mpf([2,3]) < iv.mpf([0,1]) )
print( iv.mpf([2,3]) < iv.mpf([3,4]) )
print( iv.mpf([2,3]) <= iv.mpf([3,4]) )
The number $x = \sqrt[3]{3}$ is very close to $y=26639450/18470763$. Suppose we want to know for sure which is bigger. Below we use floats, but what if there is a round off error issue?
x = 3.0**(1/3)
x
y = 26639450/18470763
y
The following makes the check with intervals:
x = iv.mpf(3)**iv.mpf("1/3")
x
y = iv.mpf("26639450/18470763")
y
print(x < y)
The intervals overlap so the above check is inconclusive. We can figure out for sure by checking higher and higher precision. (Quantities computed at higher precision will have smaller interval lengths.)
result = (x < y)
while result is None:
iv.prec = 2 * iv.prec # Double the bits of precision.
x = iv.mpf(3)**iv.mpf("1/3") # Compute x to the new precision.
y = iv.mpf("26639450/18470763") # Compute y to the new precision.
result = (x < y)
print("The statement x<y is {}.".format(result))
print("We had to use {} bits of precision.".format(iv.prec))
iv.dps = 5 # set the precision back to 5 digits
Lets say we really wanted to know $$x = \cos\big(\sqrt{3}-\sqrt{2}\big)$$ to $30$ significant digits.
To get a rough idea of what is going on, we can compute the quantity using the math library:
import math as m
m.sqrt(3) - m.sqrt(2)
Observe that the most significant digit is the tenths place (value $10^{-1}$). The 30th most significant digit then has value $10^{-30}$, so we want to find an approximation within $10^{-30}$ of the true value.
The following computes an interval containing the true value at $5$ digits of precision.
iv.dps = 5
x = iv.cos(iv.sqrt(3) - iv.sqrt(2))
x
The .delta property of an interval gives the width of the interval (as an interval).
x.delta
So we can see that we do not have enough precision. Below we increase the precision to $53$ bits of precision (double arithemetic) and try again.
iv.prec = 53
x = iv.cos(iv.sqrt(3) - iv.sqrt(2))
x.delta
We still only have 16 digits of precision or so. We need around $2$ times more precision.
iv.prec = 120
x = iv.cos(iv.sqrt(3) - iv.sqrt(2))
x.delta
This means that $x$ should be correct to 36 digits of precision. Perfect.
Here we convert the interval to a string with 30 digits of precision. This rounds each of the endpoints.
s = iv.nstr(x,30)
s
We can extract the number with s[1:33]. This will be the substring starting in position $1$, up to and not including the $32$nd position. (The $0$th position is [, so the number begins at index $1$. The first character in the number 0 is not a significant digit, and neither is the decimal point. So, we chose $33$ since it is $1 + 2 + 30$, where $1$ is the index where the number starts, $2$ is the number of characters that aren't significant digits, and $30$ is the number of significant digits.)
string_representation = s[1:33]
string_representation
The endpoints can be accessed as intervals using x.a and x.b:
print("Left: {}".format(x.a))
print("Right: {}".format(x.b))
You can access the endpoints as floating point numbers at the same precision using mpmath's mp module. You should be careful to set it to use the same precision.
from mpmath import mp
mp.prec = iv.prec
left = mp.mpf(x.a)
left
right = mp.mpf(x.b)
right
Then we can print 30 significant digits for the average of left and right to get a best approximating value. For this we use mp.nstr.
mp.nstr( (left+right)/2, n=30)
Mainly, interval arithmetic is useful to rule out the possibility of accumulation of errors due to rounding in operations. It also is useful for computing mathematical quantities to a desired accuracy.