Page - 154 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python 3.6, Volume Second Edition

154 6 ComputingIntegralsandTestingCode Rounding Errors When we compute with real numbers, these numbers are inaccuratelyrepresentedonthecomputer,andarithmeticoperationswith inaccurate numbers lead tosmall roundingerrors in thefinal results. Dependingon the typeof numericalalgorithm, the roundingerrorsmayormaynotaccumulate. Testing with aTolerance If wecannotmake tests like0.1 + 0.2 == 0.3,what should we then do? The answer is that we must accept some small inaccuracy and makea test witha tolerance.Here is the recipe: In [1]: a = 0.1; b = 0.2; expected = 0.3 In [2]: computed = a + b In [3]: diff = abs(expected - computed) In [4]: tol = 1E-15 In [5]: diff < tol Out[5]: True Here we have set the tolerance for comparison to 10−15, but calculating0.3 - (0.1 + 0.2) shows that it equals-5.55e-17, so a lower tolerance could be used in this particular example. However, in other calculations we have little idea about howaccurate theanswer is (therecouldbeaccumulationofroundingerrors inmore complicated algorithms), so 10−15 or 10−14 are robust values. As we demonstrate below, these tolerancesdependonthemagnitudeof thenumbersin thecalculations. AbsoluteandRelativeDifferences Doinganexperimentwith10k+0.3−(10k+ 0.1 + 0.2) for k = 1,.. .,10 shows that the answer (which should be zero) is around 1016−k. This means that the tolerance must be larger if we compute with larger numbers. Setting a proper tolerance therefore requires some experiments to see what level of accuracy one can expect. A way out of this difficulty is to work with relative instead of absolute differences. In a relative difference we divide by oneof theoperands,e.g., a=10k+0.3, b= (10k+0.1+0.2), c= a−b a . Computing this c for various k shows a value around 10−16. A safer procedure is thus touse relative differences. We mayexemplify this ina quicksession,usingk=10, In [1]: a = 10**10 + 0.3 In [2]: b = 10**10 + 0.1 + 0.2 In [3]: diff = a-b In [4]: diff Out[4]: -1.9073486328125e-06 In [5]: rel_diff = (a-b)/a