円周率

円周率


円周率、π(パイ)で表される3.141592654...は無理数であり、その中でも方程式の根となりえない超越数と呼ばれる数となっている
ここでは、unityの標準機能を利用せずに純粋にコードにより円周率を求めてみる

まず円の中にある正六角形を「最初のドミノ」として考え漸化式を作る事から始める

円周率を求める起点となる考え方


基本的な考え方は以下になる。円周率そのものは円周と直径との"比"を表しているので、その直径はきりの良い数字にしておくと扱いやすくなる
ここではその半径を1とした
pi1.png
ここで、この三角形をより細かくしていく事を考える。たとえばこの三角形を2分割して正12角形になった状態を考える
pi2.png
最初は6角形だった。次に12角形を考えた。このように粗い状態から細かい状態へと砕いて、よりなめらかな24角形、48角形...と角数を増やしていくと
やがて真円に近い円周の近似値が得られるだろうと考える(これは積分の考え方に近い。つまり細かく砕いてなめらかにして統合する)

三角形の円周を求めるための肝心の一辺は?ですぐには分らないが、この?に12角形なら12掛ければ円周の長さが求められる
この辺の長さを知るには色々な方法がある。三平方の定理やsin、cos等の三角形の余弦定理等の利用が思いつく

円周率を見積もる


漸化式の関数を作る前に円周率を先に見積もってみる。考え方として上記の六角形の考えを延長し三角比を利用する
この円に内接する6角形と外接する6角形の間の数が円周率の値になる筈だと考える。図で示すと以下になる
pi3.png
A を外接する6角形の一辺の 1/2。B を内接する6角形の一辺。半径 r を 1 として、鋭角30度の直角三角形の比率を利用しながら、この図を数式化すると以下になる
求めた辺の長さに円周を求めるような倍率を最後に掛算すると円周率が出る
daum_equation_1422335494722.png
これにより円周率は 3 ~ 3.4641...の間の値になる事が見積もられた

この考えを応用すればθを細かくした状態でcosθを使って辺の長さを出し細分数を掛ければ円周率を簡単に求められる
という方法も考え付くが、これはちょっと置いておいて三平方の定理を利用して漸化式で求める方法を考える

円周率を求める漸化式


<漸化式の特徴>
まず、漸化式を作る前に、漸化式自体の特徴を整理しておく
前の計算結果が続く計算の根拠となる相互の関係を漸化関係と言う。この漸化関係を表す式を漸化式と呼ぶ

  • 同じ関数を前の計算結果で変数を更新しながら繰り返し実行する事で目的の未知数を手に入れる関数式 = 漸化式
  • 数式内には実行ごとに前の計算結果で更新される変数がある
  • 数式内には更新されない定数部がある
  • その数式の「目的を果たすための動作」を決定するのは「固定された式そのものの形と定数部」

この「更新される部分」と「更新されない部分」を強く意識して先ほどの三角形のパーツを見てみる
anが前回の計算結果、an+1が今回の計算結果となっている。図形内に直角三角形がふたつ二組あるのが確認できる。この組みになった直角三角形同士の関係を利用する
pi4.png
この関係を直角三角形と三平方の定理の関係を利用しながら数式化していく
三平方の定理はお馴染みの 斜辺c=対辺a+隣辺b の関係を表す式
daum_equation_1422520037481.png
これを利用して式を作ると以下になる。この 3 つの式を連立させ変形して漸化式にする
daum_equation_1422462096347.png

<式の変形に関するプランニング>
  1. an+1 の根を求める事が目的
  2. x と y は漸化式(関数)として式内には残らないように消滅させる必要がある
    (数学的関数は入力と引数が 1:1 になる事が前提。つまり入力 an、出力 an+1 となる。入力が an , x , y と3つになるようなことは無い)
  3. 従って式は最終的に an+1=○an となる
  4. an+1 を左辺とした時、式の変形の最中に右辺に an+1 があるのは望ましくない
    式の変形計算中に右辺に y を置いた時、この y には an+1 が関わってくるので右辺から y を消滅させることが第一目標となる

daum_equation_1422496123623.png

これでan+1の根を得る式が手に入った。早速unityで計算させてみる。値は扱える桁数の多いdouble型を利用したいので.NET側のMathライブラリを利用する

<PIsolve1 ver0.1>

using UnityEngine;
using System.Collections;
using System;

public class PIsolve1 : MonoBehaviour
{
      int calcLoop;
       double an, pi;

      void Start ()
       {
               calcLoop = 20;
               an = 1d;
               for (int i = 0; i < calcLoop; i++) {
                       an = Math.Sqrt (2d - Math.Sqrt (4d - an * an));
               }
               pi = an * 6d * Math.Pow (2d, calcLoop);    //○角形の辺の数だけ掛算する
               print (pi / 2d);
       }
}


出力
3.14167426502176

本来出力される筈の正確な値 WolframAlpha - pi
3.1415926535897932384626433832795028841971693993751058

漸化式は20回程ループさせる。値はdouble型なので浮動小数点有効15桁まで正常な値が得られると期待できる。計算させてみると・・・
どうも値の様子がおかしい。有効数字、5桁目から値が違う。計算回数からして精度的に、もう少し先まで合っている筈だ
何故こうなってしまうのか?原因を探ってみる為、一回づつPocketCasで手計算をしてみる。それが以下のもの
2015-01-30 00.45.52.png
この計算式の各段階のan+1の値を見ていると、その値が回数を追う毎に小さな値になっていく事がわかる

ここで計算機誤差に対して考察する。右記のページを参照 「計算機:誤差の問題
この考えを今回に適用すると「大きさの近いふたつの数字の引算の操作」が根号内で発生している事に気が付く。これが誤差の原因となっているようだ
この誤差を無くすため式に有理化を行うと以下の式になる。これをunityのコードにする
pi9.png

using UnityEngine;
using System.Collections;
using System;

public class PIsolve1 : MonoBehaviour
{
      int calcLoop;
       double an, pi;

      void Start ()
       {
               calcLoop = 25;
               an = 1d;
               for (int i = 0; i < calcLoop; i++) {
                       an = an / (Math.Sqrt (2 + Math.Sqrt (4 - an * an)));
               }
               pi = an * 6d * Math.Pow (2d, calcLoop);    //○角形の辺の数だけ掛算する
               print (pi / 2d);
       }
}

出力
3.14159265358979
本来出力される真値 WolframAlpha - pi
3.1415926535897932384626433832795028841971693993751058

漸化式のループ回数を25回にして全ての桁が真値となった。無事、円周率を求める事が出来た事になる

近似値


355/113 ≒ pi

  • 最終更新:2015-02-17 19:51:03

このWIKIを編集するにはパスワード入力が必要です

認証パスワード