OCamlでニュートン法(2)

さて、前回微分の計算をするところまでで終わっていたんだけど今回は実際に Newton-Raphson で実装してみよう。
$\Large \mathbf{x_{n+1}=x_{n}-\frac{f(x_{n})}{f\prime(x_{n})}}$

(* small number *)
let dx = 1.0e-10;;

(* threshold *)
let threshold = 1.0e-15;;

(* derivative *)
let deriv f x dx = (f (x +. dx) -. f x) /. dx;;

(* X n+1 *)
let xn_plus_one f x = x -. (f x) /. (deriv f x dx);;

(* newton-raphson *)
let newton_raphson f i =
  let rec newton x =
    let next_num = xn_plus_one f x in
    Printf.printf "%.20f\n" next_num;
    if abs_float (next_num -. x) < threshold then ()
    else newton next_num
  in newton i;;

実行すると

# #use "newton.ml";;
val dx : float = 1e-10
val threshold : float = 1e-15
val deriv : (float -> float) -> float -> float -> float = <fun>
val xn_plus_one : (float -> float) -> float -> float = <fun>
val newton_raphson : (float -> float) -> float -> unit = <fun>
# let f x = x ** 2.0 -. 2.0;;
val f : float -> float = <fun>
# newton_raphson f 2.0;;
1.50000004137018194683
1.41666667356169995173
1.41421568647730544477
1.41421356237427708891
1.41421356237309514547
1.41421356237309492343
- : unit = ()
# newton_raphson f 10.0;;
5.10001781366008355434
2.74608272554679633259
1.73719668537143467901
1.44423846224395591165
1.41452565448231015743
1.41421359644891775353
1.41421356237308581960
1.41421356237309514547
1.41421356237309492343

このように２の平方根が求まるが、見ればわかるとおり二分法と比べて大幅に収束が速くなっている。ただし、誤差の関係上thresholdの値を1.0e-16以下にしてしまうと収束しなくなってしまうという問題もある。