RE: Calculating PI

• To: mathgroup at smc.vnet.net
• Subject: [mg28634] RE: [mg28592] Calculating PI
• From: "David Park" <djmp at earthlink.net>
• Date: Thu, 3 May 2001 04:28:25 -0400 (EDT)
• Sender: owner-wri-mathgroup at wolfram.com

```Luc,

It just happened that I worked a problem along this line several days ago.
It is from Harold M. Edwards' Advanced Calculus.

I have pasted the exercise in the form of a notebook expression at the end
of this posting. Just copy it and paste it into a blank notebook.

David Park

> From: Luc Wastiaux [mailto:luc-w at usa.net]
To: mathgroup at smc.vnet.net
> I would like to get started in mathematica programming..
>
> where can I find an algorithm that calculates PI (graphical method) ?
>
>
> --
> Luc Wastiaux - luc-w at usa.net
>

Notebook[{

Cell[CellGroupData[{
Cell["Calculating Pi", "Title"],

Cell["\<\
David Park
\>", "Subtitle"],

Cell[BoxData[
\(Needs["\<Graphics`Colors`\>"]\)], "Input"],

Cell[TextData[{
"This is Exercise 1 in Section 2.3 of Harold M. Edwards'  ",
StyleBox["Advanced Calculus: A Differential Forms Approach",
FontSlant->"Italic"],
"."
}], "Text"],

Cell[CellGroupData[{

Cell["Exercise 1 - Computation of  \[Pi]", "Subsubsection"],

Cell[TextData[{
"Calculate an approximation to ",
Cell[BoxData[
\(\[Pi] \[Equal] \[Integral]\_D \[DoubleStruckD]x\[Wedge]\
\[DoubleStruckD]y\)]],
" where ",
Cell[BoxData[
\(D\)]],
" is the domain defined by ",
Cell[BoxData[
\({{x, y} : x\^2 + y\^2 \[LessEqual] 1}\)]],
". We do this by dividing the unit square containing a quarter of \
the circle into an ",
Cell[BoxData[
\(n\[Times]n\)]],
" mesh and counting the squares definitely in the circle, the \
squares and the squares which are uncertain in that they straddle the \
boundary of the circle."
}], "Text"],

Cell["\<\
I tried using DensityPlot, but it can't be controlled precisely \
enough for this exercise. So, we will generate our own graphics. The \
following routine will generate the colored, outlined square for our \
graphics.\
\>", "Text"],

Cell[BoxData[
\(\(square[n_]\)[p_,
q_] := \[IndentingNewLine]Module[{d = 1/n, color, r1, r2,
poly}, \[IndentingNewLine]poly =
Polygon[
d {{p, q}, {p + 1, q}, {p + 1, q + 1}, {p, q + 1}, {p,
q}}]; \[IndentingNewLine]r1 =
d \@\( p\^2 + q\^2\); \[IndentingNewLine]r2 =
d \@\(\((p + 1)\)\^2 + \((q + 1)\)\^2\); \
\[IndentingNewLine]color = \[IndentingNewLine]Which[\
\[IndentingNewLine]r2 \[LessEqual] 1,
LightBlue, \[IndentingNewLine]r1 < 1 \[And] r2 > 1,
Tomato, \[IndentingNewLine]True,
White]; \[IndentingNewLine]{color, poly, Black,
Line @@ poly}\[IndentingNewLine]]\)], "Input"],

Cell["\<\
This makes our plot. If we go to too high a value of n, the grid \
lines dominate the picture, and also it takes a long time.\
\>", "Text"],

Cell[BoxData[
\(\(Module[\[IndentingNewLine]{n = 30}, \[IndentingNewLine]Show[
Graphics[\n\t{Table[\(square[n]\)[p, q], {p, 0, n - 1}, {q,
0, n - 1}], \[IndentingNewLine]Circle[{0, 0},
1, {0, \[Pi]/2}]}], \n\n
AspectRatio \[Rule] Automatic, PlotRange \[Rule] Automatic,
Background \[Rule] Linen, \n{Frame \[Rule] True,
FrameLabel \[Rule] {x, y},
PlotLabel \[Rule] "\<Approximating \[Pi]\>",
ImageSize \[Rule] 500}]];\)\)], "Input",
GeneratedCell->False],

Cell["\<\
We notice that at the end of each row (or the top of each column) we \
can have almost any number of tomato squares. So, how do we count the \
certain squares in a given row, say row q? Find the first p such that \
the top right hand corner is outside the circle.\
\>", "Text"],

Cell[CellGroupData[{

Cell[BoxData[{
\(\(1\/n\) \@\(\((p + 1)\)\^2 + \((q + 1)\)\^2\) >
1\), "\[IndentingNewLine]",
\(\(n  # &\) /@ %\), "\[IndentingNewLine]",
\(\(#\^2 &\) /@ %\), "\[IndentingNewLine]",
\(\(# - \((1 + q)\)\^2 &\) /@ %\), "\[IndentingNewLine]",
\(\(\@# &\) /@ % // PowerExpand\), "\[IndentingNewLine]",
\(\(# - 1 &\) /@ %\)}], "Input"],

Cell[BoxData[
\(\@\(\((1 + p)\)\^2 + \((1 + q)\)\^2\)\/n > 1\)], "Output"],

Cell[BoxData[
\(\@\(\((1 + p)\)\^2 + \((1 + q)\)\^2\) > n\)], "Output"],

Cell[BoxData[
\(\((1 + p)\)\^2 + \((1 + q)\)\^2 > n\^2\)], "Output"],

Cell[BoxData[
\(\((1 + p)\)\^2 > n\^2 - \((1 + q)\)\^2\)], "Output"],

Cell[BoxData[
\(1 + p > \@\(n\^2 - \((1 + q)\)\^2\)\)], "Output"],

Cell[BoxData[
\(p > \(-1\) + \@\(n\^2 - \((1 + q)\)\^2\)\)], "Output"]
}, Open  ]],

Cell["\<\
Taking the floor of the above number gives the last p column that is \
completely in the circle. Since the columns are numbered starting \
with zero we have to add 1 and the number of blue squares in row q \
for a partition n is...\
\>", "Text"],

Cell[BoxData[
\(numC[q_, n_] :=
Floor[\@\(n\^2 - \((1 + q)\)\^2\)]\)], "Input"],

Cell["\<\
To obtain the number of uncertain squares in a row we find the first \
p such that the lower left hand corner is on or outside the circle.\
\>", "Text"],

Cell[CellGroupData[{

Cell[BoxData[{
\(\(1\/n\) \@\(p\^2 + q\^2\) \[GreaterEqual]
1\), "\[IndentingNewLine]",
\(\(#  n &\) /@ %\), "\[IndentingNewLine]",
\(\(#\^2 &\) /@ %\), "\[IndentingNewLine]",
\(\(# - q\^2 &\) /@ %\), "\[IndentingNewLine]",
\(\(\@# &\) /@ % // PowerExpand\)}], "Input"],

Cell[BoxData[
\(\@\(p\^2 + q\^2\)\/n \[GreaterEqual] 1\)], "Output"],

Cell[BoxData[
\(\@\(p\^2 + q\^2\) \[GreaterEqual] n\)], "Output"],

Cell[BoxData[
\(p\^2 + q\^2 \[GreaterEqual] n\^2\)], "Output"],

Cell[BoxData[
\(p\^2 \[GreaterEqual] n\^2 - q\^2\)], "Output"],

Cell[BoxData[
\(p \[GreaterEqual] \@\(n\^2 - q\^2\)\)], "Output"]
}, Open  ]],

Cell["\<\
The p column number for the first white square is the ceiling of this \
number. This is also the total number of certain and uncertain \
squares. Subtracting the number of blue squares leaves the number of \
uncertain squares (tomato).\
\>", "Text"],

Cell[BoxData[
\(numU[q_, n_] :=
Ceiling[\@\(n\^2 - q\^2\)] - numC[q, n]\)], "Input"],

Cell["Let's check this out for some cases. Fill in the value of n.", \
"Text"],

Cell[CellGroupData[{

Cell[BoxData[
\(With[\[IndentingNewLine]{n = 10}, \[IndentingNewLine]Show[
Graphics[\n\t{Table[\(square[n]\)[p, q], {p, 0, n - 1}, {q,
0, n - 1}], \[IndentingNewLine]Circle[{0, 0},
1, {0, \[Pi]/2}]}], \n\nAspectRatio \[Rule] Automatic,
PlotRange \[Rule] Automatic,
Background \[Rule] Linen, \n{Frame \[Rule] True,
FrameLabel \[Rule] {x, y},
PlotLabel \[Rule] "\<Approximating \[Pi]\>",
ImageSize \[Rule] 500}]; \[IndentingNewLine]{certain @@
Table[numC[q, n], {q, 0, n - 1}],
uncertain @@ Table[numU[q, n], {q, 0, n - 1}]} //
TableForm]\)], "Input",
GeneratedCell->False],

Cell[BoxData[
InterpretationBox[GridBox[{
{\(certain[9, 9, 9, 9, 8, 8, 7, 6, 4, 0]\)},
{\(uncertain[1, 1, 1, 1, 2, 1, 1, 2, 2, 5]\)}
},
RowSpacings->1,
ColumnSpacings->3,
RowAlignments->Baseline,
ColumnAlignments->{Left}],
TableForm[ {
certain[ 9, 9, 9, 9, 8, 8, 7, 6, 4, 0],
uncertain[ 1, 1, 1, 1, 2, 1, 1, 2, 2, 5]}]]], "Output"]
}, Open  ]],

Cell[TextData[{
"We can now calculate the value of our approximation.  We give the \
uncertain squares a weighting on ",
Cell[BoxData[
\(1/2\)]],
"."
}], "Text"],

Cell[BoxData[
\(\[Pi]approximation[n_] :=
4/n\^2 \((Plus @@ Table[numC[q, n], {q, 0, n - 1}] + \(1\/2\)
Plus @@ Table[
numU[q, n], {q, 0, n - 1}])\)\)], "Input"],

Cell[CellGroupData[{

Cell[BoxData[
\(\[Pi]approximation /@ {10, 20, 100, 500, 1000} //
N\)], "Input"],

Cell[BoxData[
\({3.1`, 3.125`, 3.1406`, 3.141336`, 3.141534`}\)], "Output"]
}, Open  ]],

Cell[TextData[{
"The area of each uncertain square is ",
Cell[BoxData[
\(1/n\^2\)]],
". It's actual weighting should be in the range 0 to 1. In our \
calculation we give it a weighting of ",
Cell[BoxData[
\(1/2\)]],
". Therefore, the maximum error from each uncertain square is ",
Cell[BoxData[
\(1/\((2  n\^2)\)\)]],
" and the total error is ",
Cell[BoxData[
\(1/\((2  n\^2)\) U\_n\)]],
". This is different than Edwards figure of ",
Cell[BoxData[
\(2/n\^2  U\_n\)]],
". "
}], "Text"],

Cell["\<\
We might as well define the function for the total number of \
uncertain squares as a function of n.\
\>", "Text"],

Cell[BoxData[
\(numU[n_] :=
Plus @@ Table[numU[q, n], {q, 0, n - 1}]\)], "Input"],

Cell[TextData[{
"Then the uncertainty for ",
Cell[BoxData[
\(n \[Equal] 10\)]],
" and ",
Cell[BoxData[
\(n \[Equal] 20\)]],
" and ",
Cell[BoxData[
\(n \[Equal] 500\)]],
" is..."
}], "Text"],

Cell[CellGroupData[{

Cell[BoxData[
\(\(\(1\/\(2  #\^2\)\) numU[#] &\) /@ {10, 20, 500} //
N\)], "Input"],

Cell[BoxData[
\({0.085`, 0.04625`, 0.001986`}\)], "Output"]
}, Open  ]],

Cell["The actual error is less than these estimates.", "Text"],

Cell[CellGroupData[{

Cell[BoxData[
\(\[Pi] - \[Pi]approximation /@ {10, 20, 500} // N\)], "Input"],

Cell[BoxData[
\({0.04159265358979303`, 0.016592653589793116`,
0.00025665358979320985`}\)], "Output"]
}, Open  ]],

Cell["\<\
Edwards suggests that we can obtain a simplified formula for the \
number of uncertain squares. Try a linear fit.\
\>", "Text"],

Cell[BoxData[
\(\(data = Table[{n, numU[n]}, {n, 2, 40, 2}];\)\)], "Input"],

Cell[CellGroupData[{

Cell[BoxData[
\(numUfit[n_] = Fit[data, \ {1, n}, n]\)], "Input"],

Cell[BoxData[
\(\(-1.0631578947368316`\) + 1.9744360902255647`\ n\)], "Output"]
}, Open  ]],

Cell[BoxData[
\(\(Plot[{numUfit[n], 2  n}, {n, 0, 40},
PlotStyle \[Rule] {Blue, Red}, Background \[Rule] Linen,
Frame \[Rule] True, FrameLabel \[Rule] {n, U\_n},
PlotLabel \[Rule] "\<Uncertain Squares and Linear Fit\>",
ImageSize \[Rule]
500, \[IndentingNewLine]Epilog \[Rule] {AbsolutePointSize[
4], Point /@ data, Blue,
Text[U\_n \[Equal] numUfit[n],
Scaled[{0.05, 0.9}], {\(-1\), 0}], Red,
Text[U\_n \[Equal] 2  n,
Scaled[{0.05, 0.8}], {\(-1\), 0}]}];\)\)], "Input"],

Cell[TextData[{
Cell[BoxData[
\(2  n\)]],
" is easier to use and appears to be an upper bound, certainly for \
large ",
Cell[BoxData[
\(n\)]],
". Therefore the total error will be less than..."
}], "Text"],

Cell[CellGroupData[{

Cell[BoxData[
\(\(1\/\(2  n\^2\)\) 2  n\)], "Input"],

Cell[BoxData[
\(1\/n\)], "Output"]
}, Open  ]],

Cell["This checks with our previous calculation of error.", "Text"],

Cell[CellGroupData[{

Cell[BoxData[
\(1/{10, 20, 500} // N\)], "Input"],

Cell[BoxData[
\({0.1`, 0.05`, 0.002`}\)], "Output"]
}, Open  ]],

Cell[TextData[{
"Using ",
Cell[BoxData[
\(n \[GreaterEqual] 100\)]],
" will guarantee 2 place accuracy. But the result is actually much \
better."
}], "Text"],

Cell[CellGroupData[{

Cell[BoxData[
\(\[Pi] - \[Pi]approximation[100] // N\)], "Input"],

Cell[BoxData[
\(0.0009926535897930577`\)], "Output"]
}, Open  ]],

Cell["\<\
The error estimate, actually the error bound, assumes that each \
uncertain square is being used in the worse case and that the errors \
for all the squares are in the same direction. To the extent that the \
uncertain squares are randomly placed with respect to the \
circumference of the circle, their errors will tend to cancel.\
\>", "Text"]
}, Open  ]]
}, Open  ]]
},
FrontEndVersion->"4.1 for Microsoft Windows",
ScreenRectangle->{{0, 1280}, {0, 943}},
WindowSize->{639, 745},
WindowMargins->{{4, Automatic}, {Automatic, 1}}
]

```

• Prev by Date: Re: ViewPoint and RealTime3D
• Next by Date: Simple indefinite integral disagrees with table
• Previous by thread: Re: Calculating PI
• Next by thread: Re: Defining a function