OCamlでニュートン法(1) 微分の計算

さて、ニュートン法とは二分法のとこでもやった $y=f(x)$ と $x$ 軸との交点を求める場合に関数の接線を用いるという方法。
$\Large f\prime(x_{n})= \frac{f(x_{n})-0}{x_{n}-x_{n+1}}$
で、 $x_{n+1}$ について解いて、
$\Large \mathbf{x_{n+1}=x_{n}-\frac{f(x_{n})}{f\prime(x_{n})}}$
これを使う。
このうち、今回は $f\prime(x_{n})$ の計算の仕方。

$f\prime(x_{n})$ を求めるのに一番簡単なのは
$\lim_{dx\to0} \left(\frac{f(x+dx) - f(x)}{dx}\right)$
において $dx$ を十分小さな値にして計算する方法。

差分法というらしい。

ということで $f(x) = x^2$ を例に実際に試してみる。解析的手法だと $f\prime(x)=2x$ であり、xが1の場合は $f\prime(1)=2$ となる。では数値計算でも同様の結果となるか？

let f x = x ** 2.0;;

let rec recderiv x dx t =
    if t < 500 then
      let deriv x dx = ( f ( x +. dx ) -. f x ) /. dx in
      Printf.printf "dx: %.17f  dy/dx: %.17f  error: %.17f\n" dx (deriv x dx) (abs_float (2.0 -. (deriv x dx)));
      recderiv x (dx /. 10.0) (t + 1)
    else
      ()
;;

としてこれを deriv.ml というファイル名で保存。
で、実行。

# #use "deriv.ml";;
val f : float -> float = <fun>
val recderiv : float -> float -> int -> unit = <fun>
# recderiv 1.0 1.0e-5 0;;
dx: 0.00001000000000000  dy/dx: 2.00001000001392981  error: 0.00001000001392981
dx: 0.00000100000000000  dy/dx: 2.00000099992436686  error: 0.00000099992436686
dx: 0.00000010000000000  dy/dx: 2.00000010108780613  error: 0.00000010108780613
dx: 0.00000001000000000  dy/dx: 1.99999998784505761  error: 0.00000001215494239
dx: 0.00000000100000000  dy/dx: 2.00000016548074155  error: 0.00000016548074155
dx: 0.00000000010000000  dy/dx: 2.00000016548074155  error: 0.00000016548074155
dx: 0.00000000001000000  dy/dx: 2.00000016548074155  error: 0.00000016548074155
dx: 0.00000000000100000  dy/dx: 2.00017780116468158  error: 0.00017780116468158
dx: 0.00000000000010000  dy/dx: 1.99840144432528155  error: 0.00159855567471845
dx: 0.00000000000001000  dy/dx: 1.99840144432528155  error: 0.00159855567471845
dx: 0.00000000000000100  dy/dx: 2.22044604925031308  error: 0.22044604925031308
dx: 0.00000000000000010  dy/dx: 0.00000000000000000  error: 2.00000000000000000
dx: 0.00000000000000001  dy/dx: 0.00000000000000000  error: 2.00000000000000000
dx: 0.00000000000000000  dy/dx: 0.00000000000000000  error: 2.00000000000000000

結果を見るに dx が小さすぎると誤差が大きくなる傾向にある。これは打ち切り誤差や丸め誤差といった問題のためのようだ。次に任意精度を使って解いてみる。OCamlではNumというモジュールが相当するようなのでこのモジュールを使って試してみよう。

#load "nums.cma";;
open Num;;

let f x = x **/ Int 2;;

let rec recnumderiv a b t =
  if t < 500 then
    let deriv x dx = ( f ( x +/ dx ) -/ f x ) // dx in
    let ans = deriv a b in
        Printf.printf "%s\n" (approx_num_fix 180 ans);
    recnumderiv a ( b // (Int 10 ) ) (t + 1)
  else
    ()
;;

これをnumderiv.mlというファイル名として実行。

# #use "numderiv.ml";;
val f : Num.num -> Num.num = <fun>
val recnumderiv : Num.num -> Num.num -> int -> unit = <fun>
# recnumderiv (Int 1) (Int 1 // Int (-5)) 0;;

結果は実行してもらえばわかるが見事に 2 に収束していく。
ただこの Num モジュール、任意精度を扱えるのはいいがすごく扱いにくいのだった。