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 djmp at earthlink.net http://home.earthlink.net/~djmp/ > 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) ? > > thanks in advance. > > -- > Luc Wastiaux - luc-w at usa.net > Notebook[{ Cell[CellGroupData[{ Cell["Calculating Pi", "Title"], Cell["\<\ David Park djmp at earthlink.net http://home.earthlink.net/~djmp/\ \>", "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}} ]