例えば、図のように1辺が1の正方形の中に半径1、
中心角90°の扇形が接している状況を考えます。

この扇形の面積と正方形の面積の比は、いうまでもなく
1×1×π×90/360 : 1×1　=　π/4 : 1になるはずです。

そこで、図の黒い点のように正方形内部に任意の座標(x,y)をとることを考えます。
この任意の座標(x,y)が円の中に存在する確率はどうなると思いますか？

実は、(x,y)が完全にランダムであれば、π/4の確率で円内部に存在するでしょう。

そこで、任意の点(x,y)を統計的にたくさんとり、
円内部の点の数を数えれば、π/4を計算することが可能です。

このような方法をモンテカルロ法といいます。
実際のモンテカルロ法のプログラムを下記に記します。

Sub main()

Dim X, Y, D As Single
Dim i, counter As Long
Dim ob1 As Object

Set ob1 = Application.ThisWorkbook.Worksheets("Sheet1")

Randomize

counter = 0

For i = 1 To 10000000

X = Rnd
Y = Rnd
D = X ^ 2 + Y ^ 2

If D < 1 Then
counter = counter + 1
End If

Next i

ob1.Cells(1, 1).Value = "円周率"
ob1.Cells(1, 2).Value = counter / i * 4

End Sub

計算結果⇒3.14132568586743

10000000点レベルで非常に長い時間計算してもこんな精度ですが、
モンテカルロ法は円の面積以外にも応用できますので、
積分の方法として知っておくには有意義でしょう。

[1回]

奇数を次のルールで並べてください。

1番目：1
2番目：1,3
3番目：1,3,5
4番目：1,3,5,7
:
:
n番目：1,3,5,7,9…2n-1

このように並べた奇数を各グループごとに足してみましょう。

1番目：1=1
2番目：1+3=4
3番目：1+3+5=9
4番目：1+3+5+7=16
:
:

皆様、お気づきですか？実は、各奇数の和はグループ番号の二乗になっています。
確認してみましょう。

1番目：1=1=1×1
2番目：1+3=4=2×2
3番目：1+3+5=9=3×3
4番目：1+3+5+7=16=4×4
:
:

さて、n番目の奇数の和というのは、1から数えてn個の奇数を足した和と同じであり、
このルールをみると、奇数の総和なんてすぐ出せます。

例えば、1+3+5+…+11は、 2n-1 = 11 ⇔ n =6より6個の奇数をたしています。
すなわち、1+3+5+…+11　= 6×6 = 36となります。

また、1+3+5+…+999までの奇数の和は、2n-1 = 999 ⇔ n = 500より、
500×500　= 250000だとすぐ計算できるでしょう。

[2回]

前回の日記では、円周率に収束する漸化式を考えました。

詳細は、円周率について-アルキメデスの方法を考えるを参考にしてください

この漸化式をアルゴリズムにしてみますと、非常に少ない量のプログラミングソースコードで高い精度の円周率を得ることが可能です。正25165824角形を用いて円周率を計算するプログラムを下記に示します。

ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー
Sub 円周率計算()

Dim P1, Q1 As Double
Dim P2, Q2 As Double
Dim i, j, k As Long

P1 = 3
Q1 = 2 * 3 ^ 0.5
j = 6

Cells(1, 1).Value = "i"
Cells(1, 2).Value = "P(n)"
Cells(1, 3).Value = "Q(n)"

Cells(2, 1).Value = j
Cells(2, 2).Value = P1
Cells(2, 3).Value = Q1

For i = 1 To 22

k = j * 2
Q2 = 2 * P1 * Q1 / (P1 + Q1)
P2 = (P1 * Q2) ^ 0.5

Cells(2 + i, 1).Value = k
Cells(2 + i, 2).Value = P2
Cells(2 + i, 3).Value = Q2

P1 = P2
Q1 = Q2
j = k

Next i


End Sub
ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー

この精度で計算しますと、

P(25165824)＜π＜Q（25165824）

3.14159265358978<π<3.14159265358981
を導くことができます。

まさにアルキメデスの考え方が、いかに偉大なものか思い知らされますよ。

※余談ですが、アルキメデスの96角形の計算結果の場合ですと

3.14103195089051<π<3.14271459964537

となります。紀元前の人間が円周率を3.14と定めたのがすごいと思うのです。

[1回]

円周率とは直径と円の周の比です。すなわち円周率をπ、直径をd、円周をLとして

L/d = π　…(1)

となるようなπの値が円周率です。このようなπを中学では3とおしえてた時代もありますけど、この数値は3.141592…程度です。

これはなぜでしょうか。そこで直径1の円に対して、円の中に内接する正六角形の周の長さP(6)と円の外に外接する正六角形Q(6)を考えてみよう。この場合いうまでもなくd = 1なので

P(6）　< L = πd　< Q(6)
P(6) < π < 　Q(6)　…(2)

という式が導けます。ここで、

P(6) =　0.5×6 = 3
Q(6) = 0.5*2/√3*6　= 6*√3/3 = 2√3

となるので、すなわち、

3 < π < 2√3 = 3.46　…(3)

が成り立ちます。円周率πは3より大きく3.46より小さい値であることが示されました。
円周率が3ということは、円ではなくて正六角形の周長にすぎないのです。

さて、当然内接および外接する正n角形のnが増大するほど、円周率計算の精度はあがります。
そこで、正n角系に拡張して考えてみましょう。
今、直径1の円に内接するn角形の周長をP(n),外接するn角形をQ(n)とすると

P(n) <　π　<　Q(n)　…(4)

が成り立ちます。三角関数で一般化すると、
P(n) = n*sin(180/n) 　…(5)
Q(n) = n*tan(180/n) …(6)
であり、したがって、
P(2n) = 2n * sin(180/2n) …(7)
Q(2n) = 2n * tan(180/2n)　　…(8)
とあらわせます。この(5)～(8)式が満たすべき漸化式を考えてみよう。

180/n = 2θとしてP(n) + Q(n)を計算すると、
P(n) + Q(n) = 2n*sin(180/n)*(1 + 1/(cos(180/n))
P(n) + Q(n) = 2n*sin(2θ)*(cos(2θ)+1)/cos(2θ)
P(n) + Q(n) = 8n*sin(θ)*cos(θ)*cos^2(θ)/cos(2θ)
P(n) + Q(n) = 8n*sin(θ)*cos^3(θ)/cos(2θ)…(9)
のようになります。さらに、P(n)*Q(n)を計算すると、

P(n)*Q(n) = 4n^2 * sin(180/n)*tan(180/n) = 4n^2 * sin^2(2θ)/cos(2θ)
P(n)*Q(n) = 4n^2 * sin^2(2θ)/cos(2θ) = 8n^2 * sin^2(θ)cos^2(θ)/cos(2θ)…(10)

のようになります。これらを用いて、次のような計算をすると、(8)、(9)、(10)より
(P(n)+Q(n))/{P(n)*Q(n)} = (2/(2n))*(cos(θ)/sin(θ)) = (2/(2n*tan(180/2n)) = 2/Q(2n)

1/Q(2n) = 0.5{1/P(n) + 1/Q(n)}
⇔Q(2n) = 2*(P(n)*Q(n))/(P(n)+Q(n))(11)

が導けます。次に1/{Q(2n)*P(n)}を計算してみましょう。すると(7)を考えると

1/{Q(2n)*P(n)}　=　1/{2n*tan(180/2n)} *1/ {n*sin(180/n)}
1/{Q(2n)*P(n)}　=　{cos(180/2n)}/{2n*sin(180/2n)} *1/ {2n*sin(180/2n)*cos(180/2n)}
1/{Q(2n)*P(n)}　=　1/{2n*sin(180/2n)} ^2 = 1/{P(2n)}^2
1/P(2n) = √{1/(Q(2n)*P(n))}⇔P(2n) =√{Q(2n)*P(n)} …(12)
　
が導けます。これらの漸化式を解析できれば円周率を導くことができるでしょう。
余談ですが、アルキメデスはP(96)およびQ(96)を計算して紀元前のころ
正確な円周率の計算に成功しています。

アルキメデスの方法を利用するなら、
P(96)＜π＜Q(96)
を解析する必要があります。これはなかなか大変です。手計算でも計算できますが、次回の日記ではExcelのプログラミングで円周率の計算を行ってみたいですね。ご期待ください。

[1回]