Ocamlで二分法(2) 平方根を求める＆二分法の汎用部品化

昨日に続いて二分法を使い、今度は平方根を求めてみる。

二分法を使う場合に
対象となる関数が連続関数であり、 $f(a) < 0$ かつ $f(b) > 0$ の場合、 $f(x) = 0$ となる根 $x$ が区間 $(a,b)$ に存在する

ということなので $\sqrt{a}$ を求めるにあたり、 $f(x) = x^2 - a$ という式において $f(x) = 0$ を満たす $x$ が解となる。（ちなみにこの式がピンとくるまでに結構時間がかかったｗ）
このように $y = 0$ となる点が解となるように関数を作るのが二分法のキモらしい
ちなみに、二分法における解の精度は $\frac{x2 - x1 }{2 ^ n}$ ということである。

では２の平方根を求めるプログラムを実装してみよう。

$f(x) = x^2 - 2$

type sign = Positive | Negative ;;

let f x = x ** 2.0 -. 2.0 ;;

let g x =
  if x >= 0.0 then
    Positive
  else
    Negative
;;

let rec bis ( x1, x2, x3, x4 ) =
  if x3 >= x4 then
    (x1, x2, x3, x4)
  else
    let xm = ( x1 +. x2 ) /. 2.0 in
        if g ( f xm ) = g ( f x1 ) then
          bis ( xm, x2, x3 + 1, x4 )
    else
      bis ( x1, xm, x3 + 1, x4 )
;;

実行してみる。

# bis (1.0, 2.0, 0, 30);;
- : float * float * int * int =
(1.41421356145292521, 1.41421356238424778, 30, 30)

今回はサクッと行った！
さてお気づきの通り、前のと比べると関数 f の定義しか変わっていない。
なので次は高階関数を使って関数 f を引数に指定できるようにしてみよう。

type sign = Positive | Negative ;;

let g x =
  if x >= 0.0 then
    Positive
  else
    Negative
;;

let rec bis ( f, x1, x2, x3, x4 ) =
  if x3 >= x4 then
    (x1, x2, x3, x4)
  else
    let xm = ( x1 +. x2 ) /. 2.0 in
        if g ( f xm ) = g ( f x1 ) then
          bis ( xm, x2, x3 + 1, x4 )
    else
      bis ( x1, xm, x3 + 1, x4 )
;;

と、これだけ。
実際にテストしてみる。

# let f x = x ** 2.0 -. 3.0;;
val f : float -> float = <fun>
# bis (f, 1.0, 2.0, 0, 30);;
- : float * float * int * int =
(1.73205080721527338, 1.73205080814659595, 30, 30)
# let f x = x ** 2.0 -. 4.0;;
val f : float -> float = <fun>
# bis (f, 1.0, 3.0, 0, 30);;
- : float * float * int * int = (1.99999999813735485, 2., 30, 30)

できた。
これまでは試行回数を指定してきたわけだけど、次に最小誤差に到達したらその時点でプログラムを終了させるようにしてみる。ここで最小誤差は 1.0e-15 とする。

type sign = Positive | Negative ;;
let min_f = 1.0e-15;;

let g x =
  if x >= 0.0 then
    Positive
  else
    Negative
;;

let rec bis ( f, x1, x2, x3 ) =
  let xm = ( x1 +. x2 ) /. 2.0 in
  if ( x2 -. x1) < min_f then
    (x1, x2, x3 + 1)
  else
    if g ( f xm ) = g ( f x1 ) then
      bis ( f, xm, x2, x3 + 1 )
    else
      bis ( f, x1, xm, x3 + 1 )
;;

さらにグラフを書けるように中点の変遷をテキストに出力するようにする。

type sign = Positive | Negative ;;
let min_f = 1.0e-15;;
let out_file = "out.txt";;

let g x =
  if x >= 0.0 then
    Positive
  else
    Negative
;;

let rec bis ( oc, f, x1, x2, x3 ) =
  let xm = ( x1 +. x2 ) /. 2.0 in
  Printf.fprintf oc "%.17f\n" xm;
  if ( x2 -. x1) < min_f then
    Printf.printf "x1: %.17f x2: %.17f xm: %.17f try: %d\n" x1 x2 xm ( x3 + 1 )
  else
    if g ( f xm ) = g ( f x1 ) then
      bis ( oc, f, xm, x2, x3 + 1 )
    else
      bis ( oc, f, x1, xm, x3 + 1 )
;;

let fbis ( f, x1, x2, x3 ) =
  let oc = open_out out_file in
    bis ( oc, f, x1, x2, x3 );
  close_out oc;
;;

で、結果が

# let f x = cos ( x /. 2.0 );;
val f : float -> float = <fun>
# fbis ( f, 0.0, 6.0, 0);;
x1: 3.14159265358979312 x2: 3.14159265358979356 xm: 3.14159265358979312 try: 54
- : unit = ()

out.txtにも中点の値が書きこまれている。ちなみにこのプログラムを書いている途中、またもや盛大にハマっていた。Printf.printf の行の最後にあるセミコロンの必要性の有無が未だにさっぱりわかっていないのだ。まだまだヘタレである。
で、Excelでグラフおこしながら思ったんだけど、収束が速いってどのくらいからを速いというのだろうか？