Newton-Raphson法と単位円内接正多角形を使った円周率近似をScalaでやってみた

やってみた

C言語縛りは面白くないので例によってScalaでやってみた

こんなやつをコードにした

まずNewton-Raphson法は次の式におけるxの値を求める

このあと登場するけど、単位円内接正多角形を使う近似では漸化式をつかっていくわけだけども
その中で平方根を計算しないといけないのでCを任意実数定数として関数g(x)を次のようにする

すると次の式におけるxを近似的に求めることになる

Newton-Raphson法ではこれを次の式を満たす様に微分を使いながら近似していく

ここまでNewton-Raphson法

次に単位円を用いて円周率を求める。
しかし、円周がわからない。
これは円周率が不明であるためである。

そこで単位円内に内接正多角形を作って円周とみなす。
これは漸化式で表現でき正n角形のnが大きくなればなるほどその周の長さが円周の長さに近づく。

その内接正多角形の一辺の長さが次の式となる

というわけでコードにしてみた

実際に書いたコード

object calcPI{

    val max = 20

    def main(args:Array[String]):Unit = {
        var Ln = newtonRaphson(2)//L_1

        var n = 0
        var pow = 2
        while( n < max ){
            println( pow * Ln )
            Ln = newtonRaphson(2 - newtonRaphson(4-Ln*Ln))
            n += 1
            pow *= 2
        }
    }

    def newtonRaphson(c:BigDecimal):BigDecimal = {
        var x:BigDecimal = c
        var dx:BigDecimal = -func(x,c) / dfunc(x)

        while( 1.0e-16 < dx.abs  ){
            x += dx
            dx = -func(x,c) / dfunc(x)
        }

        return x
    }

    def func(x:BigDecimal, c:BigDecimal):BigDecimal = x*x - c
    def dfunc(x:BigDecimal):BigDecimal = 2.0*x
}

これを実行して、出力結果をdatファイルにおとしてgnuplotでグラフにするとこんな感じになる

BigDecimal使ってはいるけど精度的にはまだまだかなと思うけどある程度できた。