9. Special numbers: \(\pi\)
[status: written, but incomplete]
9.1. Motivation, prerequisites, plan
9.1.1. Motivation
Our purpose here is a recreational tour through the various aspects of the number \(\pi\), including calculating it with the use of random numbers. This brings out various features of probability and statistics. This topic also gives us the opportunity to write several small programs that calculate \(\pi\) and do related things. It also allows us to go mention, in passing, various “higher math” aspects of \(\pi\) which get students comfortable with the terminology of those areas of math.
But this is not entirely recreational: the number \(\pi\) is deeply tied to the geometry of circles and spheres, and thus comes up in many physical laws: just think of anything that comes from a single point and spreads uniformly in all directions: the radiation from a (rather idealized) star, the electrical force field of an isolated electron, the gravitational field of an isolated mass. The geometrical aspect of these fluxes means that at a distance \(r\) from the origin you will see an intensity that is proportional to \(\frac{1}{4\pi r^2}\)
Among the physics formulae that involve \(\pi\) you will see:
Einstein’s equation for general relativity:
\[G_{\mu\nu} = \frac{8\pi G}{c^4} T_{\mu\nu}\]The Heisenberg uncertainty principle for position and momentum in quantum mechanics:
\[\Delta x \Delta p \geq \frac{h}{4\pi}\]Since \(\pi\) occurs so much in quantum mechanics, we have even coined a special constant and symbol called \(\hbar\) (h-bar):
\[\hbar = \frac{h}{2\pi}\]so that the uncertainty principle can be written as:
\[\Delta x \Delta p \geq \frac{\hbar}{2}\]
9.1.2. Prerequisites
The 10-hour “serious programming” course.
The “Data files and first plots” mini-course in Section 2
Random number basics Section 11
9.1.3. Plan
Our plan here is not a deeply rooted one: this is a playful romp through various aspects of \(\pi\), with the intention of following whatever tangents might come up.
9.2. A collection of factoids
https://www.livescience.com/34132-what-makes-pi-special.html
among a collection of random whole numbers, the probability that any two numbers have no common factor — that they are “relatively prime” — is equal to \(6/\pi^2\). Strange, no?
Finally, \(\pi\) emerges in the shapes of rivers. A river’s windiness is determined by its “meandering ratio,” or the ratio of the river’s actual length to the distance from its source to its mouth as the crow flies. Rivers that flow straight from source to mouth have small meandering ratios, while ones that lollygag along the way have high ones. Turns out, the average meandering ratio of rivers approaches — you guessed it — \(\pi\).
https://www.angio.net/pi/whypi.html
https://www.newyorker.com/tech/annals-of-technology/pi-day-why-pi-matters
\(\pi\) is an irrational number: it can never be written as a fraction. In fact it is even more elusive: it can also never be written as the solution to a polynomial equation, so we say that it is not just irrational: it is also transcendental.
But there are a lot of cute approximations to \(\pi\) with fractions. The two which you might remember are:
and
For those who know some trigonometry, Machin’s formula is exact:
Another favorite of mine is Stirling’s approximation which gives an approximation to the factorial function (see also Section 20):
Stirling’s formula works for large-ish values of \(n\), and approximates the exact value of \(n!\) as \(n \rightarrow \infty.\) If you try it you will find that \(10! = 3628800\), while Stirling’s formula gives \(3598695.61893\), which is off by about 0.8%.
9.3. Calculating \(\pi\): ancient history
(FIXME must improve sources on this.)
Bible: fountain which is 30 paces around and 10 paces across.
Archimedes: inscribe polygons; go up to 6; he went up to 96. Depending on the energy level of the classroom, try to do the actual calculation!
9.4. Calculating \(\pi\): monte carlo method
This is a good first introduction to monte carlo integration, which allows us to discuss monte carlo methods in general.
Introduce the European city of gambling, refer to James Bond, and then dive in to the method.
The method involves shooting darts into a square which has a circle inscribed in it. Draw the picture of a circle inside a square, and draw points of random darts hitting it.
The fraction of darts that fall in the circle is proportional to the fraction of areas:
The area of the circle is \(\pi r^2\), and that of the square is \(\pi l^2\). Since we have constructed this so that \(l = 2r\) we get:
This gives us:
Now what I usually do is write a live program which has a loop that
throws 1000 darts. It does so by calculating x = random.random() *
2 - 1 and the same for y. This gives us a dart in a square. Then
using the pythagoras theorm with if sqrt(x*x + y*y) < 1 we can
determine if the dart is in the circle. We add all that up and
estimate \(\pi\).
I usually write this program (12 lines) live while the students write it with me. I write the program so that it prints, for each dart, four things: the index of the loop, the x coordinate, the y coordinate, and the estimate of \(\pi\) so far.
After experimenting with 1000 darts, then 100000, then a million, we go back to 1000 and redirect the output into a file.
This file can be plotted with a line using columns 1 and 4 (estimate of \(\pi\) vs. n_darts), and with points using columns 2 and 3 (the locations of the darts).
9.5. Calculating \(\pi\): series that converge to \(\pi\)
Further reading: https://en.wikipedia.org/wiki/List_of_formulae_involving_%CF%80
Remember the terminology: a sequence \(\{x_i\}\) is an ordered list of numbers with a criterion that gives you the next one in the sequence. An infinite series \(\sum_{k=0}^{\infty} x_k\) is the sum of the sequence \(\{x_k\}\). (Although you can then also say that the sequence of partial sums (i.e. up to N instead of up to infinity) of a series is a sequence, so the terms interweave…)
Over the years people have discovered many infinite series that converge to \(\pi\).
9.5.1. Madhava-Leibniz series
Properties: this series is exact, but it converges very slowly (sublinear convergence). To get 10 digits you need five billion terms of the series. See https://en.wikipedia.org/wiki/Leibniz_formula_for_%CF%80
9.5.2. “Efficient” infinite series
The Bailey-Borwein-Plouffe formula has the interesting feature that it can be used to pick out any binary digit of \(\pi\) (although the digit extraction algorithm is lengthy):
Let us write a program which calculates this series up to a certain value of N:
#! /usr/bin/env python3
import sys
import math
def main():
"""Calculate a series up to a certain value of k which is given on the
command line
"""
# series_term = series_term_leibniz
series_term = series_term_zeta_2
verbose = True
try:
lower = int(sys.argv[1])
upper = int(sys.argv[2])
except:
print('usage: %s lower_bound upper_bound' % sys.argv[0])
sys.exit(1)
sum_to_here = 0
for k in range(lower, upper+1):
sum_to_here += series_term(k)
# pi_approximation = 4 * sum_to_here
pi_approximation = math.sqrt(6 * sum_to_here)
if verbose: # print the intermediate values
print(k, ' ', sum_to_here, ' ', pi_approximation,
' ', math.fabs(math.pi - pi_approximation), ' ',
math.fabs(math.pi - pi_approximation) / math.pi)
print('result after %d iterations: %g' % (upper, sum_to_here))
## FIXME: the formula below changes for different series
print('## pi_approximation_final:', pi_approximation)
def series_term_leibniz(k):
"""Returns the kth term of the leibniz series for pi. Sum starts at
1. Note that to print the resulting calculation of pi you will
need to print 4*sum_to_here.
"""
return (-1)**(k+1) / (2*k - 1)
def series_term_zeta_2(k):
"""Returns the kth term of the riemann function zeta(2). Sum starts
at 1. Note that to print the resulting calculation of pi you will
need to print sqrt(6.0*sum_to_here)."""
return 1.0/(k**2)
## you can add more series term functions here
if __name__ == '__main__':
main()
Note that the results of these series often need to then be
square-rooted, or squared, or divided by something, so make sure to
modify the main() function slightly each time so that you can see
\(\pi\) clearly in the results. By default I have set it to multiply the
summation by 4 for the Leibniz formula.
Try running this program and see how rapidly it converges to \(\pi\). The Leibniz formula seems to require a factor of 10 more iterations to get just one more digit, which is quite slow.
Right now series-to-N.py has (at least) two displeasing
qualities: it asks you for the lower limit of the sum, and it
requires that you change the code to print the resulting value of
\(\pi\) depending on which series you are summing. Modify the
program so that both those pieces of information are associated
with the function, instead of having to give them as input or
changing the code.
For each of the infinite series you have programmed, make a plot of the convergence to \(\pi\) versus how many terms you sum. Then research the theoretical formulae that tell you how rapidly these series converge.
9.5.3. Formulae based on the Riemann zeta function
In general the Riemann zeta function is:
Evaluated at \(s = 2\) we get the series:
and at \(s = 4\) we get:
Write a function that calculates the Riemann zeta function for complex values of \(s\) and reproduce the attractive domain coloring plots shown in the images in https://en.wikipedia.org/wiki/Riemann_zeta_function
You will need to research how to calculate the zeta function for complex values.
9.6. Relationships between special numbers
One of the most striking math identities is Euler’s identity:
It seems to combine the five most memorable numbers (0, 1, i, e, \(\pi\)) with the three basic arithmetic operations (addition, multiplication, exponentiation).
To see how mathematicians and philosophers have gone poetic on Euler’s identity, see the wikipedia page: https://en.wikipedia.org/wiki/Euler%27s_identity#Mathematical_beauty
In class we can introduce a brief discussion of the Taylor series and show how the formula comes about, possibly adding the trigonometric identity to it. This will usually depend on a reading of the math level (and fatigue) of the students.