Sources of Numerical Errors

Contents

Join the Discussion
← Preface
→ Convergence and Stability

Highlights

• Cancellation error
• Round of error
• Overflow error

We’ll start our journey of numerical methods with an introduction to numerical error. It’s a most important topic you’ll always encounter while using numerical algorithms. There are various sources of numerical errors. Some of the notables are:

Cancellation error

It’s also called accumulation or estimated error.

Standard double-precision float (i.e. FP64 or float64) can handle only up to 15-17 significant digits.

import sys
print("Standard Python float:", sys.float_info.dig, "significant digits")

import numpy as np
print("Numpy float64:", -int(np.log10(np.finfo(np.float64).resolution)), "significant digits")
print("Numpy float32:", -int(np.log10(np.finfo(np.float32).resolution)), "significant digits")

Let’s see an useful example:

The accurate equation for relativistic kinetic energy is $$E_{k_1} = \frac{mc^2}{\sqrt{1-\frac{v^2}{c^2}}}-mc^{2}$$

If the ratio of $\frac{v^2}{c^2} < 10^{-15}$ then, we will observe an arithmetic underflow. Given $v = 2$ and $c = 3\times 10^{8}$, we get numerically

v, c = 2, 3e8

print("v^2/c^2 =", v**2/c**2)
print("1 - (1-v^2/c^2) =", (1-v**2/c**2) - 1)

We expect from $1 - (1 - \frac{v^2}{c^2}) = 4.\overline{44}\times 10^{-17}$. This means the dependence of $v$ is lost. In this example, if we say there is no dependency of $v$ (i.e. $v$ is automatically becomes 0) then, the equation reduces to $$E_{k_1}\to\frac{mc^2}{1}-mc^2=0$$ which is painfully wrong.

Mitigation: Due to the cancellation errors sometimes it is better to use approximation rather than the full formula $$E_{k_2} = mc^2\left(1-\frac{v^2}{c^2}\right)^{-1/2}-mc^2\approx mc^2\left(1-\frac{v^2}{2c^2}-1\right)=\frac12mv^2$$ But this approximation is wrong for large values of $v$. Is there are any ways to have one formula which works in both cases?

Solution: With a bit of magic (short multiplication formulas) $$E_{k_3} = mc^2\left(\frac{c}{\sqrt{c^2-v^2}}-1\right) = mc^2 \left[\frac{\left(\frac{c}{\sqrt{c^2-v^2}}\right)^{2}-1}{\frac{c}{\sqrt{c^2-v^2}}+1} \right]$$ $$ = mc^2 \left(\frac{\frac{c^2}{{c^2-v^2}}-1}{\frac{c}{\sqrt{c^2-v^2}}+1} \right) = \frac{mv^2}{\sqrt{1-\frac{v^2}{c^2}} + 1-\frac{v^2}{c^2}} $$ The above formula is mathematically equivalent with the definition but cancellation does not produce errors for small values of $v$.

Let’s see how the above formula behaves numerically for different values of $v$:

import numpy as np

def E_k_1(m, c, v):
  return m*c**2/np.sqrt(1-v**2/c**2)-m*c**2

def E_k_2(m, c, v):
  return 0.5*m*v**2

def E_k_3(m, c, v):
  g = 1-v**2/c**2
  return m*v**2/(np.sqrt(g)+g)

print("For (v, c, m) = (2, 3x10^8, 1), we have:")
v, c, m = 2, 3e8, 1
print(r"E_k_1 = ", E_k_1(m, c, v), "\nE_k_2 = ", E_k_2(m, c, v), "\nE_k_3 = ", E_k_3(m, c, v), "\n")

print("For (v, c, m) = (1x10^5, 3x10^8, 1), we have:")
v, c, m = 1e5, 3e8, 1
print(r"E_k_1 = ", E_k_1(m, c, v), "\nE_k_2 = ", E_k_2(m, c, v), "\nE_k_3 = ", E_k_3(m, c, v), "\n")

Round of error

Let me show you a prime example:

print("0.3+0.2-0.5 =",  0.3+0.2-0.5, ". This is zero.")
print("0.1+0.2-0.3 =",  0.1+0.2-0.3, ". This is not zero. Beware of round of errors!")

If we subtract same values but represented in different floating point values then, they give stupid result. For example: Given a number 1.2345678 such that $n_1$ is in float32 and $n_2$ is in float64. Then, the difference between them is

import numpy as np

n_1 = np.float32(1.2345678)
n_2 = np.float64(1.2345678)

print("n_2 - n_1 = {:.55f}".format(n_2 - n_1))

But, we were expected to get 0.0. So, we should not use two different floating point values. Instead use the same floating point value for all the calculations.

Overflow error

Let me show you a prime example:

import numpy as np

x = 20
print("For x = 20, (1-cosh(x))/(1+cosh(x)) =", (1-np.cosh(x))/(1+np.cosh(x)), "\nThis doesnot produce overflow error.")

x = 800
print("For x = 800, (1-cosh(x))/(1+cosh(x)) =", (1-np.cosh(x))/(1+np.cosh(x)), "\nThis produces overflow error.")
print("For x = 800, (1/cosh(x)-1)/(1/cosh(x)+1) =", (1/np.cosh(x)-1)/(1/np.cosh(x)+1), "\nThis produces overflow error.")