Consider the quadratic equation
$$\tag*{1} a x^{2} + b x + c = 0.$$This has two roots that can be calculated in two different ways
$$\tag*{2} x_{1,2} = \frac{-b \pm \sqrt{b^{2} -4ac}}{2a}\quad\mbox{or}\quad x_{1,2}' = \frac{-2c}{b \pm \sqrt{b^{2} -4ac}}.$$By keeping a and b constant while decreasing the value of c, and investigating the ability for the computer to evaluate the algorithms above, one can see the effect of subtractive cancellation error propagation. The square root and its preceding term almost cancel out but according to error propagation, the error of this value gets large. Consider the fractional error,
$$\tag*{3} \begin{align} \frac{r_{c} - r} {r} \simeq \frac{s}{r}\; \mbox{max}(|\epsilon_{s}|,|\epsilon_{t}|). \end{align}$$If r ends up small, but s is relatively large, the fractional error is also large. The first equation (2) is therefore calculating a small number with error order 1 by 1, and therefore doesn't suffer too much from the effects of subtractive cancellation but the second equation involves division of a small number by a smaller number with error order 1 and therefore will diverge.
In [6]:
from prettytable import PrettyTable
import numpy as np
def root_finder():
a = 1
b = 1
crray = np.geomspace(0.1, (10 ** (-200)), num=200, endpoint=False)
t = PrettyTable(['root1', 'root2', 'root3', 'root4'])
for c in crray:
root = np.sqrt((b * b) - (4 * a * c))
#print(root)
r1 = (-b + root) / (2 * a)
r2 = (-b - root) / (2 * a)
r3 = -2 * c / (b + root)
if b - root == 0:
r4 = 'div0'
else:
r4 = -2 * c / (b - root)
t.add_row([r1, r2, r3, r4])
print(t)
root_finder()
The code above shows this effect (compare columns 2 and 4). Notice how the the last column diverges. Additionally the second equation appears to be better at calculating the 'non-cancelling' root, returning a small value rather than zero, because of the nature of error propagation when dividing.
Next consider three ways of representing an expression in terms of sums.
$$\tag*{4} S_{N}^{(1)} = \sum_{n=1}^{2N} (-1)^{n} \frac{n}{n+1}.$$$$\tag*{5} S_{N}^{(2)} = - \sum_{n=1}^{N} \frac{2n-1}{2n} + \sum_{n=1}^{N} \frac{2n}{2n+1} .$$$$\tag*{6} S_{N}^{(3)} = \sum_{n=1}^{N} \frac{1}{2n(2n+1)} .$$If we assume the last one is correct and calculate how the other two differ, we run into subtractive cancellation issues, causing increasingly large errors for large N.
In [10]:
import numpy as np
import matplotlib.pyplot as plt
def sum_1(n):
outsum = 0
outarray = np.zeros((n,))
ii = (i for i in range(1, 2*n+1) if i % 2 != 0)
for i in ii:
outsum += ((-1) ** i) * i / (i + 1)
outsum += ((-1) ** (i+1)) * (i+1) / (i+1 + 1)
b = (int((i+1)/2))-1
outarray[b] = outsum
return outarray
def sum_2(n):
outsum = 0
outarray = np.zeros((n,))
for i in range(1, n+1):
outsum -= (2 * i - 1) / (2 * i)
outsum += (2 * i) / (2 * i + 1)
outarray[(i-1)] = outsum
return outarray
def sum_3(n):
outsum = 0
outarray = np.zeros((n,))
for i in range(1, n+1):
outsum += 1 / (2 * i * (2 * i + 1))
outarray[(i-1)] = outsum
return outarray
n = 1000000
method_1 = sum_1(n)
method_2 = sum_2(n)
method_3 = sum_3(n)
error_1 = abs(np.subtract(method_1, method_3))
error_2 = abs(np.subtract(method_2, method_3))
#print(method_1)
#print(method_2)
#print(method_3)
relerror_1 = np.divide(error_1, method_3)
relerror_2 = np.divide(error_2, method_3)
mil = np.arange(1, n+1, 1)
#print(mil)
logerror_1 = np.log10(relerror_1)
logerror_2 = np.log10(relerror_2)
logmil = np.log10(mil)
plt.xlabel = ('logN')
plt.ylabel = ('SN')
#plt.plot(logmil, logerror_1, '-b', lw=1)
plt.plot(logmil, logerror_2, '-r', lw=1)
plt.grid(True)
plt.show()
The code above generates log-log graphs for N against relative error for the first two sums. Plotting them together makes it hard to distinguish between the two so I kept them separate (one is currently commented out). Generally both graphs show that the log error is linearly dependent on log N. The error seems to be more unstable for large N for the second sum. This demonstrates the dangers that come with calculating expressions via large sums. If the summation includes an alternating sign, we are effectively calculating the sum of many differences, each difference with a relative error discussed before. As N increases, we are taking the difference between larger and larger numbers.
This is most easily seen in the first sum. The difference is small, however the fractional error is equal to the product of the errors in storing the input numbers and from equation (3),
$$ \frac{s}{r}$$s is of order 1 and r is small.
In [ ]: