練習問題 9.1 ニュートン法

(任意精度の)有理数計算を行う Num モジュールを使って、3.6.5 節で取り上げた Newton-Raphson 法のプログラムを書き直しなさい。いくつか使うであろう関数を下に挙げておきます。

この問題を解くまえに、Newton-Raphson法について考えなくてはならない。

最初この本を読んだときには、ややこしそうなので、読み飛ばしたのだった。
しかし、この練習問題を解くために、読み直さなくてはならない。
ところが、やはり書かれていることがわかりにくかったので、ニュートン法について説明してあるページで 勉強し、その知識をもとにプログラムを書いた。

Newton-Raphson 法のプログラム

x * x – 2 = 0 となるような関数 f(x) = x*x – 2 について考える。

x(n+1) = x(n) - f(x(n)) / f'(x(n))
       = x(n) - (x(n) * x(n) - 2) / 2 * x(n)

       n=1 のときの x(1) の値を、仮に 5 とすると、

       x(2) = x(1) - (x(1) * x(1) - 2) / 2 * x(1)
            = 5    - ( 5   *  5   - 2) / 2 * 5
            = 5    -       23          / 10
            = 5    -       2.3
            = 2.7

        x(3) = x(2) - (x(2) * x(2) - 2) / 2 * x(2)
             = 2.7  - (2.7  * 2.7  - 2) / 2 * 2.7
             = 2.7  - (7.29        - 2) / 5.4
             = 2.7  -  5.29             / 5.4
             = 2.7  -  0.97962963
             = 1.72037037

f(x) の x のときの微分係数を求める
dx を 0.1 * 10 のマイナス10乗 として考える

let deriv f =
    let dx = 0.1e-10 in
    fun x -> (f(x +. dx) -. f(x)) /. dx;;
(* 実行例 *)
deriv (fun x -> x *. x) 1.0;;
- : float = 2.000000165480742 

f(x) = 0 となるような関数 f について
x が s のときの値を求める。
すなわち、x1 が s のときの x2 をもとめる。

let nr f s =
    s -. (f s) /. (deriv f s);;
(* 実行例 *)
nr (fun x -> x *. x) 1.0;;
- : float = 0.500000041370182058

f(x) = 0 となる x をもとめる
f — f(x) = 0 となる f(x)
s — newton式に与える初期値
(f(x)の形によって適当に与えることとする)
dx — どこまで近似値をもとめるか。
0.1e-6 まで求めることとする。
x2 – x1 の値が、dx より小さくなったら解とする。

let rec newton f s =
    let dx = 0.1e-6 in
    let x1 = nr f s in
    let x2 = nr f x1 in
    if abs_float(x2 -. x1) > dx
    then newton f x1
    else x1;;

ルート5 を求める

newton (fun x -> x *. x -. 5.0) 1.0;;
(* - : float = 2.23606797750364272 *)

ルートを求める関数 heihokon
(ただし、引数は実数で)

let heihokon = newton (fun x -> x *. x -. n) 1.0;;

参考

ニュートン法について(メモ)

ニュートン法とは何か??ニュートン法で解く方程式の近似解