Lab assignment #2: Numerical Errors and Integration

Lab assignment #2: Numerical Errors and Integration
General Advice
The topics of this lab are to put into practice knowledge about numerical errors, and the
trapezoid and Simpson’s rules for integration.
• Ask questions if you don’t understand something in this background material: maybe
we can explain things better, or maybe there are typos..
• Test your code as you go, not just when it is finished. The easiest way to test code
is with print(’’) statements. Print out values that you set or calculate to make sure
they are what you think they are.
• Recall, one way to practice modularity is to define external functions for repetitive
tasks. For example, you may want to write generic functions for Trapezoidal and
Simpson’s rules (or use the piece of code provided by the textbook online resources),
place them in a separate file, and call and use them in your answer files.
Computational Background
Solving integrals numerically Lecture notes, as well as Sections 5.1 – 5.3 of the text,
introduce the Trapezoidal Rule and Simpson’s Rule. The online resource for the textbook
provides the Python program trapezoidal.py, which you are free to use.
Reading data from a textfile The command
a = numpy.loadtxt('filename.txt')
will read data in the file ’filename.txt’ into the array a.
1
Standard deviation calculations To calculate the sample mean x and standard deviation σ of a sequence {x1, . . . , xn} using the standard formulas
x ≡
1
n
Xn
i=1
xi and σ ≡
vuut
1
n − 1
Xn
i=1
(xi − x)
2
(1)
requires two passes through the data: the first to calculate the mean and the second to
compute the standard deviation (by subtracting the mean in the term (xi−x) before squaring
it). An alternative, mathematically equivalent, formula for the standard deviation,
σ ≡
vuut
1
n − 1
 Xn
i=1
x
2
i − nx
2
!
, (2)
might seem preferable because, if you think about it, you can calculate σ in eqn. (2) with
only a single pass through the data. This means, for example, that you could calculate these
statistics on a live incoming data stream as it updates. However, the one-pass method has
numerical issues that you will investigate in Q1.
The canned routine numpy.std(array, ddof=1) also implements eqn. (1). You can
think of as a “correct” calculation.
Numerical (roundoff) error: Read Section 4.2 of Newman’s textbook, which discusses
characteristics of machine error. One important point is that you can treat errors on numerical calculations as random and independent. (This is a little confusing because you will find
that errors on a given computer are often reproducible: they’ll come out the same if you do
the calculation the same way multiple times on your computer. But there is typically no way
to predict what this error value is, i.e. it could be different on different computers, or even if
you do a software update on your computer.) As a result of this, you can use standard error
statistics as in experimental physics to figure out how error propagates through numerical
calculations. This results in expressions like (4.7) in the text, describing the error on the
sum (series) of N terms x = x1 + . . . + xN :
σ = C
√
N
p
x
2
, (3)
where C ≈ 10−16 is the error constant defined on p.130 and the overbar indicates a mean
value. This means that the more numbers we include in a given series, the larger the error
(by O(
√
N)). Even if the mean error is small compared to the individual terms the fractional
error
σ
P
i
xi
=
C
√
N
p
x
2
x
(4)
can be really large if the mean value x is small, as is the case when we sum over large
numbers of opposite sign.
Page 2
Scipy constants Thanks to scipy.constants, you don’t have to look up fundamental
constants anymore! Follow the link below to see how: https://docs.scipy.org/doc/
scipy/reference/constants.html
Physics Background
Black body radiation The black body function can be written as a function of wavenumber ν and temperature T, using the Planck constant h, the speed of light c and the Boltzmann
constant k:
B =
2hc2
ν
3
e
hcν
kT − 1
(5)
The total energy per unit area emitted by a black body follows Stefan’s law:
W = σT4
, (6)
where σ is the Stefan-Boltzmann constant. (Of course, you remember all this if I was your
prof for PHY252.)
Questions
1. [25%] Exploring numerical issues with standard deviation calculations
(a) Write pseudocode to test the relative error found when you estimate the standard
deviation using the two formulas (1) and (2), treating the numpy method
numpy.std(array, ddof=1)
as the correct answer. The input for this calculation will be a supplied dataset consisting of a one-dimensional array of values that is read in using numpy.loadtxt().
Note: The relative error of a value x compared to some true value y is (x − y)/y.
In implementing (2), you will need to account for the possibility that this approach
could result in taking the square root of a negative number. (Stop and think about
why this check is necessary.) You can implement a stopping condition if this is
the case, or print a warning.
Submit the pseudo-code.
(b) Now, using this pseudocode, write a program that uses (1) to calculate the standard deviation of Michelsen’s speed of light data (in 103 km s−1
), which is stored in
the file cdata.txt, and which was taken from John Baez’s (U.C. Riverside) website:
http://math.ucr.edu/home/baez/physics/Relativity/SpeedOfLight/measure_
c.html
Calculate the relative error with respect to numpy.std(array, ddof=1). Now do
the same for eqn. (2). Which relative error is larger in magnitude?
Submit the code, printed output, and written answer.
Page 3
(c) To explore this question further, we will evaluate the standard deviation of a
sequence with a predetermined sample variance. The function
numpy.random.normal(mean, sigma, n)
returns a sequence of length n of values drawn randomly from a normal distribution with mean mean and standard deviation sigma. Now generate two normally
distributed sequences, one with
mean, sigma, n = (0., 1., 2000)
and another one with
mean, sigma, n = (1.e7, 1., 2000),
with the same standard deviation but a larger mean. Then evaluate the relative
error of (1) and (2), compared to the numpy.std call. How does the relative error
behave for the two sequences?
Now that you have investigated a few cases, can you explain the difference in the
errors in the two methods, both for these distributions and for the data in Q1b?
Submit printed output and written answer.
(d) Can you think of a simple workaround for the problems with the one-pass method
encountered here? Try this workaround and see if it fixes the problem.
Submit pseudo-code, code, printed output and written answer.
2. [25%] Trapezoidal and Simpson’s rules for integration
We seek to evaluate the integral
Z 1
0
4
1 + x
2
dx. (7)
(a) What is the exact value of this integral?
(b) For N = 4 slices, compare the value you obtain when using the Trapezoidal vs.
Simpson’s rule, and with the exact value.
(c) For each method (Trapezoidal and Simpson), how many slices do you need to
approximate the integral with an error of O(10−9
)? We are only looking for a
rough estimate for the number of slices for each method (i.e., N = 2n
, with
n = 2, 3, 4 . . . , and N or n being the answer). For this question, do it in a
“dumb” way: increase the number of slices until you hit the mark. How long does
it takes to compute the integral with an error of O(10−9
) for each method?
Note: the integration is fast in both methods. In order to get an accurate timing,
you may want to repeat the same integration many times over. Recall that we
described a timing method in last week’s lab, but you are free to choose functions
with better performance.
(d) Adapt the “practical estimation of errors” of the textbook (§ 5.2.1, p. 153) to
the trapezoidal method only to obtain the error estimation for N2 = 32 (using
N1 = 16).
Page 4
(e) Why wouldn’t the practical estimation method work with the Simpson’s rule in
our particular case? How (in a few words; no need to implement it) would we
need to adapt the practical estimation of errors method for it to work with the
Simpson’s rule in our particular case?
For the whole exercise on this integral, submit your pseudo-code, code,
printed outputs, and written answers.
3. [25%] Stefan-Boltzmann constant
(a) Show that in eqn (5), the integration of B over ν can be written as
W = π
Z ∞
0
Bdν = C1
Z ∞
0
x
3
e
x − 1
dx. (8)
What is in C1?
Submit your written answer.
(b) Write a program to calculate the value for W given the temperature T. Explain
the method used to integrate over the infinite range, and give an estimate for the
accuracy of the method.
Submit your pseudocode, code, and written answer.
(c) Use your answer above to derive a value for the Stefan-Boltzmann constant (see
eqn 6) in SI units, to three significant figures or more. Check your result against
the value given in scipy.constants.
Submit your code, and written answer.
4. [25%] Exploring roundoff error
The idea of this exercise is to explore the effects of roundoff error for a polynomial
calculated a couple of ways. Consider p(u) = (1 − u)
8
in the vicinity of u = 1.
Algebraically, this is equivalent to the following expansion:
q(u) = 1 − 8u + 28u
2 − 56u
3 + 70u
4 − 56u
5 + 28u
6 − 8u
7 + u
8
. (9)
But numerically, p and q are not exactly the same.
(a) On the same graph, plot p(u) and q(u) very close to u = 1, for example picking
500 points in the range 0.98 < u < 1.02. Which plot appears noisier? Can you
explain why?
Submit your graph
(b) Now plot p(u)−q(u), and the histogram of this quantity p(u)−q(u) (for u near 1).
Do you think there is a connection between the distribution in this histogram and
the statistical quantity expressed in eqn (3) above? To check, first calculate the
standard deviation of this distribution (you can use the std function in numpy).
Then calculate the estimate obtained by using equation (3), with C = 10−16
.
Page 5
State how you calculated the other terms in (3). Hint: We are looking for order
of magnitude consistency here, do not worry about O(30 − 50%) differences.
Submit your plots and written answer
(c) From equation (4) above, show that for values of u somewhat greater than
0.980 but less than 1.0, the error is around 100%. Verify this by plotting
or printing out abs(p-q)/abs(p) for u starting at u = 0.980 and increasing
slowly up to about u = 0.984 (it might be different on different computers).
This fractional error quantity is noisy and diverges quickly as u approaches 1.0,
so you might need to plot or print several values to get a good estimate of the
values of u at which the error approaches 100%.
Submit your plots and written answer
(d) Roundoff error doesn’t just apply to series, it also comes up in products and
quotients. For the same u near 1.0 as in parts a-b, calculate the standard deviation (error) of the numerically calculated quantity f = u**8/((u**4)*(u**4)).
This quantity will show a range of values around 1.0, with roundoff error. You
can get a sense of the error by plotting f-1 versus u. Compare this error to the
estimate in equation (4.5) on p.131 of the textbook (don’t worry about the factor
of √
2).
Submit your plot(s) and written answer
Page 6