History of Computational Physics

In [1]:
%matplotlib inline
%config InlineBackend.print_figure_kwargs={'bbox_inches':None}
import numpy as np
import matplotlib.pyplot as pl
import matplotlib.patches as pt
In [2]:
%%html
<style>
.annotate {color: tomato}
</style>

Introduction

Computing has been used to give answers to pressing questions for a very long time. As the tools and means to execute calculations improved, more and more calculations have been possible. Calculations have been done by humans on paper, by scribbling in sand, with tally sticks, or by using mechanical devices, such as the abacus, or the sliding rule. With such devices one can do maybe 2 or 3 additions or multiplications of moderate size numbers per minute. Only after the invention of modern computers in 1940's the computational science has gained enormous power, when thousands of operations per seconds become possible. Today's computing devices measure their number processing speed in GFlop (billion of floating point operations per second). As of 2021, it is believe that in few years supercomputers will break into ExaFlop regime.

The ever changing nature of computing technology makes defining computational physics hard, because new and more challenging problems start to become feasible with the newest improvements.

This notebook illustrates one of the earliest achievement in computational science: the exact calculation of the mathematical constant $\pi$. It is common knowledge for us that $\pi = 3.1415926\ldots$, but 2000 years ago we only knew that $\pi$ is about 3. Here is how Archimedes did it in his famous work Κύκλου μέτρησις (Kuklou metrēsis) ca. 250 BC. Although technology changed radically, the guiding principles, the methodology and the spirit behind this calculation remained the same to this day.

What is the question here?\ What are the digits that come after 3 in the number $\pi$?

What is the theoretical framework?\ Basic plane geometry of regular polygons is the only knowledge required to follow these calculations.

Mathematical Formulation\ The main idea is that the length of the circle is larger that the perimeter of any polygon inscribed of it. The length of circle is also smaller that the perimeter of any polygon within which the circle is inscribed.

Numerical tools and techniques\ The calculation of $\pi$ is impossible to any number of digits (the current record is few billion digits, please check in report back, google for example "record pi digits"), but with this procedure lower and higher bounds can be calculated by iteration, placing $\pi$ within tighter and tighter bounds, leading to convergent results.

What do we learn from this?\ Even when results are not known exactly, or impossible to know, approximate results are satisfactory for all practical purposes. The results of most calculations that require $\pi$ don't change sensibly if $\pi$ is considered with 12 or 15 digits of precision. However, in rare cases when it may matter, it is reassuring that we have a procedure that can produce extra digits, of course, at some computational extra cost.

How long is the circle?

The number $\pi$ is the ratio between the length of a circle and its diameter. But how to measure the length of circle? Archimedes' idea was to circumscribe and inscribe regular polygons with more and more sides. At the first step one uses square, the next step octagons and so on. Each step uses polygons with double the number of sides than the previous polygon. That means that the angle at the center of the circle subtended by one side is halved at each step.

In [17]:
fig = pl.figure(figsize=(6,6))
ax = fig.add_axes([0.05,0.05,0.9,0.9])
ax.axis([-6.7,6.7,-6.7,6.7])
ax.axis("off")
ax.set_aspect(1)

ax.add_patch(pt.Circle((0,0),5, fill = False, lw = 2))
ax.add_patch(pt.Circle((0,0),0.1, facecolor="black"))
ax.add_patch(pt.RegularPolygon((0,0), 6, radius=4.9,
        fill=False, edgecolor="tab:red", lw=3, orientation=30*np.pi/180))
ax.add_patch(pt.RegularPolygon((0,0), 12, radius=5,
        fill=False, edgecolor="tab:orange", lw=1, ls="--",orientation=30*np.pi/180))
ax.add_patch(pt.RegularPolygon((0,0), 4, radius=5.09*pow(2,0.5),
        fill=False, edgecolor="tab:blue", lw=3, orientation=-0.5*45*np.pi/180))
ax.add_patch(pt.RegularPolygon((0,0), 8, radius= 5.45,
        fill=False, edgecolor="tab:green", lw=1,ls="--", orientation=-0.5*45*np.pi/180))
ax.text(0.1,0.1,"O", fontsize=14)
ax.text(-4.5,4.2,"$P_4$", fontsize=14, color="tab:blue")
ax.text(-1.5,5.3,"$P_8$", fontsize=14, color="tab:green")
ax.text(-1,3.7,"$p_6$", fontsize=14, color="tab:red")
ax.text(0.,4.4,"$p_{12}$", fontsize=14, color="tab:orange");

This figure illustrates a hexagon $p_6$ and a dodecagon $p_{12}$ inscribed in a circle, and also a square $P_4$ and a octogon $P_8$ circumscribed to the same circle. The drawing itself is created with Python code, you can peak inside and see the commands for each graphic element, lines, text, etc.

Our first goal is to calculate the sizes of the sides of inscribed polygons. The sequence of polygons goes like $p_6$, $p_{12}, \dots p_{6\times 2^k}$, with $k = 0, 1, 2, \ldots$. At every step, the number of sides doubles and the size of each side $\ell$ becomes smaller and smaller, until eventually $\ell \rightarrow 0$. The figure below shows what happens during one transition. The blue line is $\ell_j = 2 AC$, while the red line is the side for the next iteration $\ell_{j+1} = AB$. Here $\ell_j$ is the size of a side of the inscribed regular polygon with $6\times 2^j$ sides, and $\ell_{j+1}$ is, of course the length of a side for the next polygon, with double number of sides.

The geometry of the triangle AOB requires that $$AC \cdot OB = AB \cdot OD$$ By assuming that the radius OB is unity in some units, we get $$\frac 12 \ell_j = \ell_{j+1} \sqrt{ 1 - (\ell_{j+1}/2)^2}$$ or $$\ell_{j+1}^4 - 4 \ell_{j+1}^2 + \ell_j^2 = 0$$ As a quadratic equation in $\ell_{j+1}^2$, it has solutions $$\ell_{j+1}^2 = 2 \pm \sqrt{4 - \ell_j^2}$$ Which solution should we keep here? Since we know that $\ell_j \rightarrow 0$, only the solution with "$-$" makes sense, so we get the recursion $$\ell_{j+1} = \sqrt{2 - \sqrt{4 - \ell_j^2}} \quad\quad (1)$$ which should start with $\ell_0 = 1$ for the size of the inscribed regular hexagon, then calculate $\ell_1$ for the size of the polygon with 12 sides, and so on.

The perimeter for each inscribed polygon that gradually approach a circle is $6\times 2^{j} \ell_j$, and therefore the number $\pi$ is greater than the half of the perimeter, $$\pi \ge 6\times 2^{j-1} \ell_j \quad\quad(2)$$

In [6]:
fig = pl.figure(figsize=(6,6))
ax = fig.add_axes([0.05,0.05,0.9,0.9])
ax.axis([-10,10,-2,18])
ax.axis("off")
ax.set_aspect(1)
ax.set_frame_on(True)

def draw_point(p, l, label_pos="SW"):
    pos = {"N":(0, 1), "NW":(-0.7, 0.7), "W":(-1,0), "SW":(-0.7,-0.7), 
           "S":(0,-1), "SE":( 0.7,-0.7), "E":( 1,0), "NE":( 0.7, 0.7)}
    w = 16
    ofst = (w*pos[label_pos][0], w*pos[label_pos][1])
    ax.annotate(l, p, xytext = ofst,ha="center",va="center",
                textcoords="offset points", fontsize=14)
    ax.plot(p[0],p[1], "o", color="tab:blue")

def draw_line(p1, p2, **kwargs):
    kwargs = dict(kwargs)
    if not 'c' in kwargs.keys() :  kwargs["c"]="black"
    ax.plot((p1[0],p2[0]),(p1[1],p2[1]), **kwargs)

R, ϕ, dϕ = 16., 30., 4.0
phi1 = (90 + ϕ)*np.pi/180
phi2 = (90 - ϕ)*np.pi/180

a1 = pt.Arc((0,0),2*R,2*R,0,90-(ϕ+dϕ),90+(ϕ+dϕ), edgecolor="gray", lw=2)
ax.add_patch(a1)

pO = (0,0)
pA = (R*np.cos(phi1), R*np.sin(phi2))
p1 = (R*np.cos(phi2), R*np.sin(phi2))
pB = (0, R)
pC = (0, R*np.sin(phi1))
pD = (0.5*R*np.cos(phi1), 0.5*R*(1 + np.sin(phi1)))

draw_point(pO,"O")
draw_point(pA,"A","N")
draw_point(pB,"B","NE")
draw_point(pC,"C")
draw_point(pD,"D","SE")
draw_point(p1,"")
draw_line(pO,pA)
draw_line(pO,pB)
draw_line(pO,p1)
draw_line(pA,p1,c="tab:blue", lw=3)
draw_line(pA,pB,c="tab:red", lw=3)
draw_line(pO,pD,c="black", ls="--")

#pl.savefig("fig.pdf")

We can carry the same procedure for circumscribed polygons $P_{4\times 2^k}$ with $k=0,1,2,\ldots$. The blue line is $L_j = 2 AB$, while the red line is the size for the next iteration $L_{j+1} = 2 CD$.

Because triangles OCD and OEB are similar we have that $CD = OC\times EB/OE = EB/OE$, because we assume that the radius OC = 1. On the other hand, $AB \times OF = 2 EB \times OE$, and therefore $CD = AB/2 \times OF/OE^2$. For the next step in reasoning we look at triangle FAB and write $$FA^2 + AB^2 = (2 EB)^2 = (OF - 1)^2 + AB^2 = OF^2 -2 OF + 1 + AB^2$$ On the other hand, in triangle OAB: $$OB^2 = AB^2 + 1$$ and since $OF=OB$ we have that $$2 AB^2 - 2 OF + 2 = 4 EB^2 \rightarrow OF = AB^2 + 1 - 2 EB^2$$ But $AB^2 + 1 = OE^2 + EB^2 = OB^2$, and therefore $$OF = OE^2 - EB^2$$

Finally, this results in having that $$CD = \frac 12 AB \times \frac{OE^2 - EB^2}{OE^2} = \frac 12 AB \left( 1 - (EB/OE)^2\right) = \frac 12 AB ( 1 - CD^2)$$, which provides the recurence $$\frac 12 L_{j+1} = \frac 1 4 L_{j} [1 - \frac 14 L_{j+1}^2]$$ that can be solved as a quadratic equation $$L_j L_{j+1}^2 + 8 L_{j+1} - 4 L_j = 0$$ to give $$L_{j+1} = \frac{4}{L_{j}} \left( \sqrt{1 + \frac{L_j^2}{4}} - 1\right)\quad\quad(3)$$ If we start from a square with $L_0 = 2$, we can calculate the size of the octogon side $L_1$ and so on. The number $\pi$ in this case is smaller than the half of the perimeter of the circumscribed polygon with $4\times 2^j$ sides $$\pi \le 4 \times 2^{j-1} L_{j} \quad\quad (4)$$

In [7]:
fig = pl.figure(figsize=(6,6))
ax = fig.add_axes([0.05,0.05,0.9,0.9])
ax.axis([-12,8,-2,18])
ax.axis("off")
ax.set_aspect(1)
ax.set_frame_on(True)

R, ϕ, dϕ = 12., 45/2., 14.0
phi = ϕ*np.pi/180
a1 = pt.Arc((0,0),2*R,2*R,0,90-(ϕ+dϕ),120+(ϕ+dϕ), edgecolor="gray", lw=2)
ax.add_patch(a1)

pO = (0,0)
pA = (-R*np.sin(phi), R*np.cos(phi))
pD = ( R*np.sin(phi)/np.cos(phi), R)
pC = (0, R)
pB = (-R*np.sin(phi)+R*np.cos(phi), R*np.cos(phi) + R*np.sin(phi))
p1 = (-R*np.sin(phi)-R*np.cos(phi), R*np.cos(phi) - R*np.sin(phi))
p2 = ( -R*np.sin(phi)/np.cos(phi), R)
pE = (0, R*np.cos(phi) + R*np.sin(phi))
pF = (-R*1.41*np.sin(phi), R*1.41*np.cos(phi))
draw_point(pO,"O")
draw_point(p2,"")
draw_point(pA,"A","S")
draw_point(pC,"C")
draw_point(pD,"D","E")
draw_point(pB,"B","E")
draw_point(pE,"E")
draw_point(pF,"F")
draw_line(pO, pF)
draw_line(pO, pB)
draw_line(pO, pE, ls="--")
draw_line(pB, p1, c="tab:blue", lw=3)
draw_line(pD, p2, c="tab:red", lw=3)
draw_line(pB, pF)

Recurrence

Using equations (1), (2), (3) and (4) we can calculate successively the sizes $\ell_j$ and $L_j$ of polygons with doubling number of sides, starting with hexagons and squares, and "squeeze" the number $\pi$ between shrinking bounds because:

$$ 6\cdot 2^{j-1} \ell_j \le \pi \le 4\cdot 2^{j-1} L_j$$

The python code below does just that, working out the first 6 levels. At the first level (level 0) we can only say that $\pi$ is between 3 and 4, half of the perimeters of inscribed hexagon and circumscribed square. However, with each iteration, the boudn become tighter and tighter.

In [3]:
l = 1
L = 2
for j in range(6):
    p = 6*pow(2,j-1)*l
    P = 4*pow(2,j-1)*L
    print("level {0:d}: {1:.16f} < π < {2:.16f}".format(j,p,P))
    l = pow(2 - pow(4 -  l*l,0.5),0.5)
    L = 4*(pow(1 + L*L/4, 0.5) - 1)/L

level 0: 3.0000000000000000 < π < 4.0000000000000000
level 1: 3.1058285412302498 < π < 3.3137084989847612
level 2: 3.1326286132812369 < π < 3.1825978780745294
level 3: 3.1393502030468721 < π < 3.1517249074292581
level 4: 3.1410319508905298 < π < 3.1441183852458665
level 5: 3.1414524722853443 < π < 3.1422236299423449

Programming Assignment\ Extend the calculation to higher levels, say up to 20. What do you observe? Is this something that you would expect? Can you explain in mathematical terms what is happening?

In [ ]: