Infinite sum of Python for loop

Question

I'm a geometry student. Recently we were doing all kinds of crazy circle stuff, and it occurred to me that I don't know why $\pi r^2$ is the area of a circle. I mean, how do I really know that's true, aside from just taking my teachers + books at their word?

So I tried to derive the formula myself. My strategy was to fill a circle with little squares. But I couldn't figure out how to generate successively smaller squares in the right spots. So instead I decided to graph just one quadrant of the circle (since all four quadrants are identical, I can get the area of the easy +x, +y quadrant and multiply the result by 4 at the end) and put little rectangles along the curve of the circle. The more rectangles I put, the closer I get to the correct area. If you graph it out, my idea looks like this:

Okay, so to try this in practice I used a Python script (less tedious):

from math import sqrt, pi

# explain algo of finding top right quadrant area
# thing with graphics via a notebook

# Based on Pythagorean circle function (based on r = x**2 + y**2)
def circle_y(radius, x):
  return sqrt(radius**2 - x**2)

def circleAreaApprox(radius, rectangles):
  area_approx = 0
  little_rectangles_width = 1 / rectangles * radius

  for i in range(rectangles):
    x = radius / rectangles * i
    little_rectangle_height = circle_y(radius, x)
    area_approx += little_rectangle_height * little_rectangles_width

  return area_approx * 4

This works. The more rectangles I put, the wrongness of my estimate goes down and down:

for i in range(3):
    rectangles = 6 * 10 ** i
    delta = circleAreaApprox(1, rectangles) - pi # For a unit circle area: pi * 1 ** 2 == pi
         print(delta)

Output

0.25372370203838557
0.030804314363409357
0.0032533219749364406

Even if you test with big numbers, it just gets closer and closer forever. Infinitely small rectangles circleAreaApprox(1, infinity) is presumably the true area. But I can't calculate that, because I'd have to loop forever, and that's too much time. How do I calculate the 'limit' of a for loop?

Ideally, in an intuitive way. I want to reduce the magic and really understand this, not 'solve' this by piling on more magic techniques (like the $\pi \times radius^2$ formula that made me curious in the first place).

Thanks!

Answer 1

This is an excellent question. You are following in Archimedes' footsteps and starting to invent integral calculus and the idea of a limit.

I will try to address (briefly!) the mathematical and philosophical issues here, not the programming question.

You are right to worry about a process that has to go on forever. The way mathematicians deal with that question is to replace the infinitely many operations it would take to "reach a limit" by infinitely many inequalities any one of which can be justified in a predictable finite number of steps. If in your picture you calculate the total area of the inscribed slices just as you have calculated the area of the circumscribed ones you can show (with logic, not a program) that the difference between those two areas is as small as you please as long as you are willing to use thin enough rectangles. Then you can argue (though it's not easy) that there is just one number less than all the overestimates and greater than all the underestimates. For a circle of radius $1$ we call that number $\pi$.

The next job is to over and underestimate the circumference of the unit circle with the same kind of argument, using circumscribed and inscribed polygons. There too you can show that they tell you a number for the circumference.

The final step is to show that number is exactly twice the $\pi$ that you found for the area.

For a circle of radius $r$ the circumference will be $r$ times as large, so $2\pi r$, and the area will be $r^2$ times as large, so $\pi r^2$. (Carefully proving those proportionalities for curved shapes like circles requires estimations and limits.)

Answer 2

In the limit as your mini squares shrink to no size, the area goes from being a sum of their areas to an integral. Consider the circle $x^2+y^2=r^2$. In the positive quadrant, one quarter of its area is$$\int_0^r\sqrt{r^2-x^2}dx=\int_0^{\pi/2}r^2\cos^2tdt$$(by substituting $x=r\sin t$). The mean values of $\cos^2t$ and $\sin^2t$ sum to $1=\cos^2t+\sin^2t$, and are equal as the functions differ only in a phase shift, so the average of each function is $\tfrac12$, and we're integrating over a half-period of $\cos^2t$, but the full period just reflects this. So the circle's area is $4\tfrac{\pi}{2}\tfrac{r^2}{2}=\pi r^2$.

This isn't how people originally worked it out, it's just what your method points us to. The simplest solution is to cut the circle into thinner and thinner sectors instead, and arrange these alternately into something resembling a rectangle, of height $r$. Its width is a half-circumference $\pi r$.

This argument can be formalised with calculus too; it boils down to the infinitesimal area element being $\rho d\rho d\theta$, with $\rho$ an arbitrary internal point's distance from the centre:$$\int_0^r\rho d\rho\int_0^{2\pi}d\theta=\tfrac12r^22\pi=\pi r^2.$$

But your approach is interesting from the perspective of more recent mathematics. A Monte Carlo method computes probabilities from integrals or vice versa. That the circle comprises a proportion $\tfrac{\pi}{4}$ of the smallest enclosing square's area means a randomly chosen point in that square has probability $\tfrac{\pi}{4}$ of being in the circle.

Answer 3

Typically this kind of issue is handled by seeing a stable value emerge, with tolerance under some desired value, from successively finer approximation.

---

In this case if you double the number of rectangles, you approximately halve the error - this is what you would expect from treating each sub-section of the circle as a straight line section, with the excess of area having a rectangle taken out of it, leaving two triangles of the same area. Anyway, this gives you a way to get extra precision by combining successive estimates.

It's actually quite interesting to look at the difference between this extrapolated estimate and a trapezoidal approach; effectively by this picture of section of the curve:

showing that here the two-rectangles-for-one loses the upper right rectangle of area so extrapolating that forward can be thought of as losing the additional wedges off the top of the two rectangles. For comparison the green trapezoid gives an underestimate of the area under the curve, in this case.

---

Using your base routines for the "extrapolated area" gives the following :

rec = 1
est = [circleAreaApprox(1, 1)]  # all estimates on unit circle
for i in range(20):
    rec *= 2
    est.append(circleAreaApprox(1, rec)) 
    print("rectangles= ",rec, "basic error=",est[-1]-pi,
          "extrapolated error=", 2*est[-1]-est[-2]-pi)

which is calculating both $(e_k-\pi)$ and $(2e_k-e_{k-1}-\pi)$, gives

rectangles=  2 basic error= 0.5904581539790841 extrapolated error= 0.32250896154796127
rectangles=  4 basic error= 0.35411641451264764 extrapolated error= 0.1177746750462112
rectangles=  8 basic error= 0.19822649076738008 extrapolated error= 0.042336567022112526
rectangles=  16 basic error= 0.10666038423794832 extrapolated error= 0.015094277708516568
rectangles=  32 basic error= 0.05600976928733781 extrapolated error= 0.00535915433672729
rectangles=  64 basic error= 0.028954259189892806 extrapolated error= 0.0018987490924478045
rectangles=  128 basic error= 0.014813138806821335 extrapolated error= 0.000672018423749865
rectangles=  256 basic error= 0.007525429367436942 extrapolated error= 0.00023771992805254882
rectangles=  512 basic error= 0.0038047491298613956 extrapolated error= 8.406889228584902e-05
rectangles=  1024 basic error= 0.0019172379492267133 extrapolated error= 2.972676859203105e-05
rectangles=  2048 basic error= 0.0009638743216493495 extrapolated error= 1.0510694071985682e-05
rectangles=  4096 basic error= 0.00048379526796571426 extrapolated error= 3.7162142820790223e-06
rectangles=  8192 basic error= 0.00024255458489719217 extrapolated error= 1.3139018286700832e-06
rectangles=  16384 basic error= 0.00012150956160006388 extrapolated error= 4.6453830293557985e-07
rectangles=  32768 basic error= 6.083690070024517e-05 extrapolated error= 1.642398004264578e-07
rectangles=  65536 basic error= 3.0447484171247652e-05 extrapolated error= 5.806764225013694e-08
rectangles=  131072 basic error= 1.5234007097131297e-05 extrapolated error= 2.053002301494189e-08
rectangles=  262144 basic error= 7.62063280745906e-06 extrapolated error= 7.258517786823404e-09
rectangles=  524288 basic error= 3.8115995590892737e-06 extrapolated error= 2.566310719487319e-09
rectangles=  1048576 basic error= 1.9062534803993003e-06 extrapolated error= 9.074017093269049e-10

and by the time you have a lot of rectangles, you are gaining $1000\times$ precision from extrapolation.