Sources of Numerical Errors

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.

$\text{Standard Python \texttt{float}: $15$ significant digits} \\ \text{Numpy \texttt{float64}: 15 significant digits} \\ \text{Numpy \texttt{float32}: 6 significant digits}$

Standard Python float

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

Numpy float, Numpy finfo

  • Other name for numpy.float64: float in python, double in C/C++
  • Other name for numpy.float32: float in C/C++
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

$\frac{v^2}{c^2} = 4.4444444444444444 \times 10^{-17}\\ 1 - (1 - \frac{v^2}{c^2}) = 0.0$

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$:

$\text{For } (v, c, m) = (2, 3\times 10^{8}, 1), \text{ we get:}\\ E_{k_{1}} = {\color{red} 0}.0 \\ E_{k_{2}} = {\color{red} 2}.0\\ E_{k_{3}} = {\color{red} 2}.0\\ \quad \\ \text{For } (v, c, m) = (1\times 10^{5}, 3\times 10^{8}, 1), \text{ we get:}\\ E_{k_{1}} = 5000000{\color{red} 416.0} \\ E_{k_{2}} = 5000000{\color{red} 000.0} \\ E_{k_{3}} = 5000000{\color{red} 416.666705}$

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:

$0.3+0.2-0.5 = 0.0 \text{. This is zero.}\\ 0.1+0.2-0.3 = 5.551115123125783\times 10^{-17} \text{. This is {\color{red}not} zero. Beware of round of errors!}$

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

$n_2 - n_1 = 0.000000038578796379695745599747$$1697628498077392578125000$

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:

$\text{For } x = 20, \frac{1 - \cosh(x)}{1 + \cosh(x)} = -0.9999999917553856\\ \text{This doesnot produce overflow error.}\\\quad\\ \text{For } x = 800, \frac{1 - \cosh(x)}{1 + \cosh(x)} = \texttt{nan}\\ \text{This {\color{red}produces} overflow error.} \\\texttt{RuntimeWarning: overflow encountered }$$ \texttt{in cosh} \\ \texttt{RuntimeWarning: invalid value }$$ \texttt{encountered in double\_scalars} \\\quad\\ \text{For } x = 800, \frac{1 - \cosh(x)}{1 + \cosh(x)} \equiv \frac{\frac{1}{\cosh(x)} - 1}{\frac{1}{\cosh(x)} + 1} = -1.0\\ \text{This {\color{red}produces} overflow error.} \\ \texttt{RuntimeWarning: overflow encountered }$$ \texttt{in cosh}$

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.")

Permalink at https://www.physicslog.com/cs-notes/nm/sources-of-numerical-errors

Published on Feb 13, 2023

Last revised on Feb 17, 2023

References