Theoretical Concepts and Simulations
S. Denisov MNTF, Uni Augsburg
1 / 27
Machine precision (continuation)
Determine your machine precision
I Determine experimentally the precision of single-precision floats
I Determine experimentally the precision of double-precision floats
2 / 27
Machine precision (continuation)
Determine your machine precision: an example (Fortran)
The program (left panel) and the results (righ panel) The results: step (first column), variable ’one−1’ (second cloumn), variable ’eps’ (third column).
3 / 27
Machine precision (continuation)
Determine your machine precision: an example (Fortran)
Almost the same program as above but now
’10000·(one−1)’ is calculated
The results: step (first column), variable ’10000·(one−1)’
(second cloumn), variable ’eps’ (third column).
4 / 27
Machine precision (continuation)
An example: summing series Consider an infinite series for sin(x):
sin(x) =x−x2 3! +x5
5! − x7 7! +...
The task is to caclulate sin(x) for a givenx with an absolute error 10−8. An algorithm is to use the finite sum
sin(x)w
N
X
n=1
(−1)n−1x2n−1 (2n−1)!
How do we decide when to stop summing?
5 / 27
Machine precision (continuation)
First idea: Calculate (−1)n−1x2n−1 and then divide it by (2n−1)!
It iswrong idea: Both are quickly grow very large even though their quotient does not.
Goodidea: both functions are products so the recursion can be used. Namely,
6 / 27
Machine precision (continuation)
To obtain an absolute error of 10−8, we stop the calculation when
In general, you are free to pick any tolerance level you desire, although if it is too close to, or smaller than, machine precision, your calculation may not be able to attain it. Below is the pseudo code for performing the summation
7 / 27
Errors in computations
Substractive cancellation
Lets model the computer representationxc of the exact number x as
xc wx(1 +x)
Herex is the relative error in xc, which we expect to be of a similar magnitude to the machine precisionm. If we apply this notation to the simple subtractiona=b−c, we obtain
a=b−c ⇒ac wbc−cc 'b(1 +b)−c(1 +c)
⇒ ac
a w1 +bb a −c
ac.
Error inais essentially a weighted average of the errors inb and c, with no assurance that the last two terms will cancel. The worst case isb 'c: we are subtracting off the most significant parts of both numbers and leaving the error-prone least-significant parts.
8 / 27
Errors in computations
Substractive cancellation (continuation) ac
a
def= 1 +a'1 +b
a(b−c)'1 +b
amax(|b|,|c|).
Multiplication by the large numberb/asignificantly magnify the error.
If you subtract two large numbers and end up with a small one, there will be less significance, and possibly a lot less significance, in the small one.
Example: Expansione−x '1−x+x2/2!−x3/3! +...for largex.
9 / 27
Errors in computations
Round-off error in a single division a= b
c ⇒ac = bc cc
= b(1 +b) c(1 +c),
⇒ ac
a = b c
(1 +b) (1 +c) w b
c(1 +b)(1−c)w b
c['1 +b−c]
⇒ ac
a = b
c[1 +|b|+|c|].
Multiplication by the large numberb/asignificantly magnify the error.
We can generalize this model to estimate the error in the evaluation of a general functionf(x),
E = f(x)−f(xc)
f(x) ' df(x)/dx
f(x) (x−xc)
10 / 27
Errors in computations
Round-off accumulation after many steps
An error on each step can be considered as a literal ˆastepˆa in a one dimensional random walk, that is, a walk for which each step is in a random direction. The total distance covered inN steps of lengthr, is, on the average,
R '√ Nr
The total relative errortot arising after N calculational steps each with machine precision errorm is, on the average,
tot '√ Nm
In some cases there may be no cancellation, and the error may increase asNm. Even worse, in some recursive algorithms, where the error generation is coherent, such as the upward recursion for Bessel functions (considered below), the error increases asN!
11 / 27
Errors in computations
Round-off accumulation after many steps (continuation) A fast computer may complete 1010 floating-point operations per second. This means that a program running for 3 h performs approximately 1014 operations. Therefore, if round-off error accumulates randomly, after 3 h we expect a relative error of 107m. For the error to be smaller than the answer, we need m <107, which requires double precision and a good algorithm. If we want a higher-precision answer, then we will need avery good algorithm.
12 / 27
Errors in computations
An example: Errors in calculations spherical Bessel functions There are two family of functions, the spherical Bessel of the first kind,jl(x), and of the second kind, yl(x) (also called ’Neumann’s functions’), which satisfy the same differential equation
Two families behave differently at the originx = 0 (nonsingular and singular). The spherical Bessel functions occur in many physical problems, such as the expansion of a plane wave into spherical partial waves
13 / 27
Errors in computations
An example: Errors in calculations spherical Bessel functions For the first twol values, their explicit forms are
Approximate values for the spherical Bessel functions of the orders 3, 5, and 8.
14 / 27
Errors in computations
An example: Errors in calculations spherical Bessel functions
Bessel functions jl(x) (left panel) and yl(x) (right panel) as functions of x.
Notice that for smallx, the values of for increasing l become progressively smaller forjl(x), and bigger for yl(x).
15 / 27
Errors in computations
An example: Errors in calculations spherical Bessel functions Numerical recursion relations
Both equations are the same relation, one written for upward recurrence from small to largel values, and the other for downward recurrence to smalll values.
To recur upward inl for fixed x, we start with the known forms for j0 and j1 (see previous slides) and use the first recurrence.
16 / 27
Errors in computations
An example: Errors in calculations spherical Bessel functions
Realization of the upward iteration idea.
17 / 27
Errors in computations
An example: Errors in calculations spherical Bessel functions Forx '2 it starts good but then fails. Why?
Lets look at the plot on the previous slide (two slides back). We see that as we recurjl up to largerl values, we are taking the difference of two “large” functions to produce a “small” value for jl. This process suffers from subtractive cancellation and always reduces the precision. As we continue recurring, we take the difference of two small functions with large errors and produce a smaller function with yet a larger error. After a while, we are left with only round-off error.
18 / 27
Errors in computations
An example: Errors in calculations spherical Bessel functions Even if we start with exactjl, after a short while the computerˆas lack of precision effectively mixes in a bit ofyl(x):
jlc =jl(x) +yl(x)
This is inevitable because bothjl andyl satisfy the same differential equation and, on that account, the same recurrence relation. The admixture ofyl becomes a problem when the numerical value ofyl(x) is much larger than that of jl(x) because even a minuscule amount of a very large number may be large (this is the case forx'2).
19 / 27
Errors in computations
An example: Errors in calculations spherical Bessel functions The simplest solution is so-calledMiller’s device: just use downward recursion, starting from some large valuel =L. This avoids subtractive cancellation by taking small values ofjl+1(x) andjl(x) and producing a larger jl−1(x) by addition. While the error may still behave like a Neumann function, the actual
magnitude of the error willdecrease quickly as we move downward to smallerl values.
20 / 27
Errors in computations
An example: Errors in calculations spherical Bessel functions If we start iterating downward with arbitrary values forjL+1(c) and jL(c), after a short while we will arrive at the correctl dependence for a given value ofx. Although the numerical values ofjl(c), l = 1, ...,L−1 so obtained will not be correct because they depend upon the arbitrary values assumed valuesjL+1(c) andjL(c), therelative values will be accurate. But because we know that
j0(x) = sin(x)/x, what is left is the normalization:
21 / 27
Errors in computations
An example: Errors in calculations spherical Bessel functions
Realization of Miller’s device (downward iteration)
22 / 27
Error analysis in numerical simulations
Algorithm performance and error scaling Accuracy of simulations depend on
I accuracy of an algorithm
I a machine precision (round-offs)
Algorithmic errortypically scales with the number of steps N like al 'αN−β
α andβ are algorithm-specific constants
Round-off errortypically scales with the number of stepsN like ro 'D·Nµ
whereµis the “diffusion exponent“ andD is the “diffusion constant”. Typicallyµ'1/2 (diffusion-like dynamics).
23 / 27
Error analysis in numerical simulations
Therefore,overall errorscales with the number of steps N as ov 'αN−β+D·Nµ
How to estimate the overall error? It is trivial when the exact solution is known (but then there is no need for simulations).
Otherwise use the following recipe:
I simulate algorithmA with N steps and get the result,X(N)
I simulate algorithmA with M =sN,s is an integer larger than 1, steps and get the result, X(M)
I calculate4(N) =kX(M)−X(N)kand plot the dependence 4(N) vsN
24 / 27
Error analysis in numerical simulations
An example: integration of nonlinear oscillator Equation of motion
¨
x=−γx˙−x−x3+asin(ωt) Two variables: coordinate x and velocity υ= ˙x.
Initial conditions: This is the state of the system at the initial instant of timet = 0,x(0) = 0 and υ(0) = 0.
Task: to propagate the system over the time t = 10T, where T = 2π/ω.
Overall error: E(N) =p
(x(sN)−x(N))2+ (υ(sN)−υ(N))2, s = 4.
Algorithms: Euler’s scheme and the 4th-order Runge-Kutta method(will discuss them latter on).
Number of stepsN is the number of integration steps per period T. It sets the time-step h=T/N.
25 / 27
Error analysis in numerical simulations
An example: integration of nonlinear oscillator
Performances of the two algorithms
26 / 27
Error analysis in numerical simulations
Good reading
David Goldberg, What Every Computer Scientist Should Know About Floating Point Arithmetic, ACM Computing Surveys23, 5 -
48 (1991).
27 / 27
