[Date Index] [Thread Index] [Author Index]

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}}
]