Chapter 4 - Approximations and errors in computing

Example No. 4_01 Pg No. 63

In [1]:
#Greatest precision

a = '4.3201'
b = '4.32'
c = '4.320106'
na = len(a)-(a.index('.')+1)
print '\n %s has a precision of 10**-%d\n'%(a,na)
nb = len(b)-(b.index('.')+1)
print '\n %s has a precision of 10**-%d\n'%(b,nb)
nc = len(c)-c.index('.')-1
print '\n %s has a precision of 10**-%d\n'%(c,nc)
e = max(na,nb,nc)
if e ==na:
    print '\n The number with highest precision is %s\n'%(a)
elif e == nb:
    print '\n The number with highest precision is %s\n'%(b)
else:
    print '\n The number with highest precision is %s\n'%(c)
 4.3201 has a precision of 10**-4


 4.32 has a precision of 10**-2


 4.320106 has a precision of 10**-6


 The number with highest precision is 4.320106

Example No. 4_02 Pg No. 63

In [2]:
#Accuracy of numbers

def sd(x):
    nd = x.index('.')  #position of point
    num = [ord(c) for c in x] #str2code(x)
    if (nd)==None and num(length(x)) == 0:
        print 'Accuracy is not specified\n'
        n = 0 #
    else:
        if num[0]>= 1 and (nd==None):
            n = len(x)
        elif num[1] >= 1 and not((nd==None)):
                n = len(x) - 1
        else:
            for i in range(0,len(x)):
                if num[i] >= 1 and num[i] <= 9:
                    break;
                
            
            n = len(x)- i + 1
    return n
        
    

a = '95.763'
na = sd(a)
print '%s has %d significant digits\n'%(a,na)
b = '0.008472'
nb = sd(b)
print '%s has %d significant digits.The leading or higher order zeros are only place holders\n'%(b,nb)
c = '0.0456000'
nc = sd(c)
print '%s has %d significant digits\n'%(c,nc)
d = '36.0'
nd = sd(d)
print '%s has %d significant digits\n'%(d,nd)
e = '3600.0'
sd(e)
f = '3600.00'
nf = sd(f)
print '%s has %d significant digits\n'%(f,nf)
95.763 has 5 significant digits

0.008472 has 7 significant digits.The leading or higher order zeros are only place holders

0.0456000 has 8 significant digits

36.0 has 3 significant digits

3600.00 has 6 significant digits

Example No. 4_03 Pg No. 64

In [3]:
from __future__ import division
from math import floor
from numpy import arange,zeros
a = 0.1
b = 0.4
afrac=[];bfrac=[];
for i in range(0,8):
    afrac.append(floor(a*2))
    a = a*2 - floor(a*2)
    bfrac.append(floor(b*2))
    b = b*2 - floor(b*2)

afrac_s = '0' + '.' +(str(afrac)) #string form binary equivalent of a i.e 0.1
bfrac_s = '0' + '.' +(str(bfrac))
print '\n 0.1_10 = %s \n 0.4_10 = %s \n '%( afrac_s , bfrac_s)
    
summ=zeros(8)
for j in arange(7,0,-1):
    summ[j] = afrac[j] + bfrac[j]
    if summ[j] > 1:
        summ[j] = summ[j]-2
        afrac[(j-1)] = afrac[(j-1)] + 1
summ_dec = 0
for k in arange(7,0,-1):
    summ_dec = summ_dec + summ[k]
    summ_dec = summ_dec*1.0/2 

print 'sum =',summ_dec
print 'Note : The answer should be 0.5, but it is not so.This is due to the error in conversion from decimal to binary form.'    
 0.1_10 = 0.[0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 1.0] 
 0.4_10 = 0.[0.0, 1.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0] 
 
sum = 0.9921875
Note : The answer should be 0.5, but it is not so.This is due to the error in conversion from decimal to binary form.

Example No. 4_04 Pg No. 66

In [4]:
#Rounding-Off

fx = 0.7526
E =3
gx = 0.835
d = E - (-1)
#Chopping Method
Approx_x = fx*10**E
Err = gx*10**(E-d)
print '\n Chooping Method : \n Approximate x = %.4f*10**%d \n Error = %.4f \n '%(fx,E,Err)
#Symmetric Method
if gx >= 0.5:
    Err = (gx -1)*10**(-1)
    Approx_x = (fx + 10**(-d))*10**E
else:
    Approx_x = fx*10**E
    Err = gx * 10**(E-d)

print '\n Symmetric Rounding :\n Approximate x = %.4f*10**%d \n Error = %.4f \n '%(fx + 10**(-d),E,Err)
 Chooping Method : 
 Approximate x = 0.7526*10**3 
 Error = 0.0835 
 

 Symmetric Rounding :
 Approximate x = 0.7527*10**3 
 Error = -0.0165 
 

Example No. 4_05 Pg No. 68

In [5]:
from scipy.misc import factorial
#Truncation Error

x = 1.0/5
#When first three terms are used
Trunc_err = x**3/factorial(3) + x**4/factorial(4) + x**5/factorial(5) + x**6/factorial(6)
print '\n a) When first three terms are used \n Truncation error = %.6E \n '%(Trunc_err)

#When four terms are used
Trunc_err = x**4/factorial(4) + x**5/factorial(5) + x**6/factorial(6)
print '\n b) When first four terms are used \n Truncation error = %.6E \n '%(Trunc_err)

#When Five terms are used
Trunc_err = x**5/factorial(5) + x**6/factorial(6)
print '\n c) When first five terms are used \n Truncation error = %.6E \n '%(Trunc_err)
 a) When first three terms are used 
 Truncation error = 1.402756E-03 
 

 b) When first four terms are used 
 Truncation error = 6.942222E-05 
 

 c) When first five terms are used 
 Truncation error = 2.755556E-06 
 

Example No. 4_06 Pg No. 68

In [6]:
from scipy.misc import factorial

#Truncation Error

x = -1.0/5
#When first three terms are used
Trunc_err = x**3/factorial(3) + x**4/factorial(4) + x**5/factorial(5) + x**6/factorial(6)
print '\n a) When first three terms are used \n Truncation error = %.6E \n '%(Trunc_err)

#When four terms are used
Trunc_err = x**4/factorial(4) + x**5/factorial(5) + x**6/factorial(6)
print '\n b) When first four terms are used \n Truncation error = %.6E \n '%(Trunc_err)

#When Five terms are used
Trunc_err = x**5/factorial(5) + x**6/factorial(6)
print '\n c) When first five terms are used \n Truncation error = %.6E \n '%(Trunc_err)
 a) When first three terms are used 
 Truncation error = -1.269244E-03 
 

 b) When first four terms are used 
 Truncation error = 6.408889E-05 
 

 c) When first five terms are used 
 Truncation error = -2.577778E-06 
 

Example No. 4_07 Pg No. 71

In [7]:
#Absolute and Relative Errors

h_bu_t = 2945 #
h_bu_a = 2950 #
h_be_t = 30 #
h_be_a = 35 #
e1 = abs(h_bu_t - h_bu_a)
e1_r = e1/h_bu_t
e2 = abs(h_be_t - h_be_a)
e2_r = e2/h_be_t
print '\n For Building : \n Absolute error,  e1 = %d \n Relative error , e1_r = %.2f percent \n '%(e1,e1_r*100) 
print '\n For Beam : \n Absolute error,  e2 = %d \n Relative error , e2_r = %.2G percent \n '%(e2,e2_r*100)
 For Building : 
 Absolute error,  e1 = 5 
 Relative error , e1_r = 0.17 percent 
 

 For Beam : 
 Absolute error,  e2 = 5 
 Relative error , e2_r = 17 percent 
 

Example No. 4_08 Pg No. 72

In [8]:
from math import log10
#Machine Epsilon

def Q(p):
    q = 1 + (p-1)*log10(2)
    return q
p = 24
q = Q(p)
print 'q = %.1f \n We can say that the computer can store numbers with %d significant decimal digits \n '%(q,q)
q = 7.9 
 We can say that the computer can store numbers with 7 significant decimal digits 
 

Example No. 4_09 Pg No. 75

In [9]:
#Propagation of Error

x = 0.1234*10**4
y = 0.1232*10**4
d = 4
er_x = 10**(-d + 1)/2
er_y = 10**(-d + 1)/2
ex = x*er_x
ey = y*er_y
ez = abs(ex) + abs(ey)
er_z = abs(ez)/abs(x-y)

print '\n |er_x| <= %.2f o/o\n |er_y| <= %.2fo/o \n ex = %.3f \n ey = %.3f \n |ez| = %.3f \n |er_z| = %.2fo/o \n'%(er_x *100,er_y*100,ex,ey,ez,er_z*100)
 |er_x| <= 0.05 o/o
 |er_y| <= 0.05o/o 
 ex = 0.617 
 ey = 0.616 
 |ez| = 1.233 
 |er_z| = 61.65o/o 

Example No. 4_10 Pg No. 77

In [10]:
#Errors in Sequence of Computations

x_a = 2.35 #
y_a = 6.74 #
z_a = 3.45 #
ex = abs(x_a)*10**(-3+1)/2
ey = abs(y_a)*10**(-3+1)/2
ez = abs(z_a)*10**(-3+1)/2
exy = abs(x_a)*ey + abs(y_a)*ex
ew = abs(exy) + abs(ez)
print '\n ex = %.5f \n ey = %.5f \n ez = %.5f \n exy = %.5f \n ew = %.5f \n'%(ex,ey,ez,exy,ew)
 ex = 0.01175 
 ey = 0.03370 
 ez = 0.01725 
 exy = 0.15839 
 ew = 0.17564 

Example No. 4_11 Pg No. 77

In [11]:
from math import floor
#Addition of Chain of Numbers

x = 9678 #
y = 678 #
z = 78 #
d = 4 #  #length of mantissa
fx = x/10**4
fy = y/10**4
fu = fx + fy
Eu = 4
if fu >= 1:
    fu = fu/10
    Eu = Eu + 1

#since length of mantissa is only four we need to maintain only four places in decimal, so
fu = floor(fu*10**4)/10**4
u = fu * 10**Eu
w = u + z
n = len(str(w))
w = floor(w/10**(n-4))*10**(n-4) #To maintain length of mantissa = 4
print 'w =',w
True_w = 10444
ew = True_w - w
er_w = (True_w - w)/True_w
print 'True w =',True_w
print 'ew =',ew
print 'er,w =',er_w 
w = 10000.0
True w = 10444
ew = 444.0
er,w = 0.0425124473382

Example No. 4_12 Pg No. 77

In [12]:
#Addition of chain Numbers

x = 9678 #
y = 678 #
z = 78 #
d = 4 #  #length of mantissa
n = max(len( str(y) ) , len(str(z)))
fy = y/10**n
fz = z/10**n
fu = fy + fz
Eu = n
if fu >= 1:
    fu = fu/10
    Eu = Eu + 1

u = fu * 10**Eu
n = max(len( str(x) ) , len(str(u)))
fu = u/10**4
fx = x/10**4
fw = fu + fx
Ew = 4
if fw >= 1:
    fw = fw/10
    Ew = Ew + 1

#since length of mantissa is only four we need to maintain only four places in decimal, so
fw = floor(fw*10**4)/10**4
w = fw*10**Ew
print 'w =',w
True_w = 10444
ew = True_w - w
er_w = (True_w - w)/True_w
print 'True w =',True_w
print 'ew =',ew 
print 'er,w =',er_w
w = 10430.0
True w = 10444
ew = 14.0
er,w = 0.00134048257373

Example No. 4_14 Pg No. 79

In [13]:
from math import sqrt
from scipy.misc import derivative
#Absolute & Relative Errors

xa = 4.000
def f(x):
    f = sqrt(x) + x
    return f
#Assuming x is correct to 4 significant digits
ex = 0.5 * 10**(-4 + 1)
df_xa = derivative(f,4)
ef = ex * df_xa
er_f = ef/f(xa)
print '\n ex = %.0E \n df(xa) = %.2f \n ef = %.2E \n er,f = %.2E \n'%( ex,df_xa,ef,er_f)
 ex = 5E-04 
 df(xa) = 1.25 
 ef = 6.26E-04 
 er,f = 1.04E-04 

Example No. 4_15 Pg No. 80

In [14]:
#Error Evaluation

x = 3.00 #
y = 4.00 #
def f(x,y):
    f = x**2 + y**2
    return f
def df_x(x):
    df_x = 2*x
    return df_x
def df_y(y):
    df_y = 2*y
    return df_y
ex = 0.005
ey = 0.005
ef = df_x(x)*ex + df_y(y)*ey
print 'ef =',ef
ef = 0.07

Example No. 4_16 Pg No. 82

In [15]:
#Condition and Stability

C1 = 7.00 #
C2 = 3.00 #
m1 = 2.00 #
m2 = 2.01 #
x = (C1 - C2)/(m2 - m1)
y = m1*((C1 - C2)/(m2 - m1)) + C1
print 'x =',x
print 'y =',y
print 'Changing m2 from 2.01 to 2.005'
m2 = 2.005
x = (C1 - C2)/(m2 - m1)
y = m1*((C1 - C2)/(m2 - m1)) + C1
print '\n x = %d \n y = %d \n From the above results we can see that for small change in m2 results in almost 100 percent change in the values of x and y.Therefore, the problem is absolutely ill-conditioned \n'%(x,y)
x = 400.0
y = 807.0
Changing m2 from 2.01 to 2.005

 x = 800 
 y = 1607 
 From the above results we can see that for small change in m2 results in almost 100 percent change in the values of x and y.Therefore, the problem is absolutely ill-conditioned 

Example No. 4_18 Pg No. 84

In [16]:
from math import sqrt, floor
#Difference of Square roots

x = 497.0 #
y = 496.0 #
sqrt_x = sqrt(497)
sqrt_y = sqrt(496)
nx = len( str( floor( sqrt_x ) ) )
ny = len( str( floor( sqrt_y ) ) )
sqrt_x = floor(sqrt_x*10**(4-nx))/10**(4-nx)
sqrt_y = floor(sqrt_y*10**(4-ny))/10**(4-ny)
z1 = sqrt_x - sqrt_y
print 'z = sqrt(x) - sqrt(y) =',z1
z2 = ( x -y)/(sqrt_x + sqrt_y)
if z2 < 0.1:
    z2 = z2*10**4
    nz = len(str(floor(z2)))
    z2 = floor(z2*10**(4-nz))/10**(8-nz)

print  'z = ( x-y )/( sqrt(x) + sqrt(y) ) =',z2
z = sqrt(x) - sqrt(y) = 0.0
z = ( x-y )/( sqrt(x) + sqrt(y) ) = 0.022

Example 4_21 Pg No. 85

In [17]:
from __future__ import division
from math import exp,floor
x = -10
T_act=[1]
T_trc=[1]
e_x_cal = 1
TE=[]
for i in range(0,100):
    T_act.append(T_act[i]*x/(i+1))
    T_trc.append(floor(T_act[i+1]*10**5)/10**5)
    TE.append(abs(T_act[i+1]-T_trc[i+1]))
    e_x_cal = e_x_cal + T_trc[i+1]
e_x_act = exp(-10)
Sum=0
for s in TE:
    Sum+=s
print 'Truncation Error =',Sum
print 'calculated e**x using roundoff =',e_x_cal
print 'actual e**x = ',e_x_act
print 'Here we can see the difference between calculated e**x and actual e**x this is due to trucation error (which is greater than final value of e**x ), so the roundoff error totally dominates the solution'
Truncation Error = 0.000505399929331
calculated e**x using roundoff = -0.000459999999793
actual e**x =  4.53999297625e-05
Here we can see the difference between calculated e**x and actual e**x this is due to trucation error (which is greater than final value of e**x ), so the roundoff error totally dominates the solution