Lecture 1 & 2 Numerical accuracy and precision issues¶
Han Q Le (c)
-copyrighted
import numpy as np
import pandas as pd
import scipy as sp
import matplotlib.pyplot as plt
Section 1 - Floating point vs. integer¶
simple equation fails to hold numerically¶
x1=5; y1=np.sqrt(5);
z=x1-y1**2
z
z is not zero as it should be: numerical error due to finite digit representation
np.log2(np.abs(z))
See that abs(z) is exactly 1/2**(-50)
print([np.abs(z), 1/(2**(50))])
why is z exactly 2^integer? see lecture on binary representation of number in FPU Also, download related apps on apps webpage.
consider a range of numbers¶
u1=np.array([[x,x-(np.sqrt(x))**2,np.log2(np.abs(x-(np.sqrt(x))**2))] for x in range(50)])
u1
x1arr=u1[:,0]; z=np.abs(u1[:,1]);
plt.plot(x1arr,z)
create a range from 0 to 1,000,000
u2=np.array([[x,x-(np.sqrt(x))**2,np.log2(np.abs(x-(np.sqrt(x))**2))]
for x in np.arange(0,1000000,2503)])
x2arr=u2[:,0]; z2=np.abs(u2[:,1]);
plt.loglog(x2arr,z2)
Note that the error grows proportionally with th number magnitude: this is because the same error in the mantissa is multiplied with the exponent: intrinsic error in the mantissa
Section 2 - sinusoidal error with large argument¶
Consider this problem
plot EM wave at d=0¶
wl=0.5; k=2*np.pi/wl ; xarr2=np.arange(0,5,0.02);
d=0;
ef1=np.cos(k*(d+xarr2))
plt.plot(xarr2,ef1)
plot EM wave at d=384000 km¶
wl=0.5; k=2*np.pi/wl ; xarr2=np.arange(0,5,0.02);
d=384000e9;
ef3=np.cos(k*(d+xarr2))
plt.plot(xarr2,ef3)
for d=1000 km, error is small, but can be seen¶
wl=0.5; k=2*np.pi/wl ; xarr2=np.arange(0,5,0.02);
d=1000e9;
ef2=np.cos(k*(d+xarr2))
plt.plot(xarr2,ef2)
Plot the difference between ef2 and ef1, we can already see loss-of-precision numerical error
plt.plot(xarr2,ef1-ef2)
The numerical error of cosine grows with argument magnitude. The correct way is to use the periodic property of the function:
wl=0.5; k=2*np.pi/wl ; xarr2=np.arange(0,5,0.02);
d=384000e9;d0=np.fmod(d,wl);
efcorrect=np.cos(k*(d0+xarr2))
plt.plot(xarr2, efcorrect, 'g', xarr2, ef3, 'r')
Similar errors¶
x=np.arange(-1,1,0.05)
sinx=np.array([np.sin(2*np.pi*2*n*1000 + x*1.e-11) for n in range(1,4)])
plt.plot(x,x*1.e-11,'b')
[plt.plot(x,sinx[i],'r') for i in range(len(sinx))]
The plot below is Sin[n 2 Pi +x] vs x where n is an integer. The result should be a straight line (blue) because: Sin[n 2 Pi +x] = Sin[x] ~ x for small x (which is on the order of 10^-11). Yet for large n, ~ few 2000, one can see errors due to loss of precision as well as biases due to a large argument.
System epsilon, largest, smallest¶
small number add to large number¶
epsilon=0.123456789e-14
epsilon
y=1+epsilon
delta=y-1
print([delta, epsilon])
x=np.arange(-1,1,0.0125); xeps=x*1e-15
y=1+xeps
plt.plot(x,y-1,'r',x,xeps,'b')
largest number and smallest number¶
x=np.arange(0,2,0.5)*1e308
np.array([x,2*x])
soft (or slow) flush to zero
x=1e-322
for i in range(6):
x=0.5*x
print(x)