fetburner.core

コアダンプ

ニュートン法とDKA法の実装

大学でニュートン法DKA法で方程式の解を求める課題が出たので、今回は趣向を変えてStandard MLで実装しました。

実装は以下の通りです。

(* 添字も渡すmap *)
fun mapi f =
  rev
    o #1
    o foldl (fn (x, (result, i)) => (f (i, x) :: result, i + 1)) ([], 0)
 
(* 回数の上限に達するか終了条件を満たすまで計算を反復する *)
fun iterateWhile fin f =
  let 
    fun iterateWhileAux acc 0 x = acc
      | iterateWhileAux acc n x =
          if fin x then acc
          else iterateWhileAux (x :: acc) (n - 1) (f x)
  in
    iterateWhileAux []
  end
 
local
  open Complex
  (* 総乗 *)
  fun producti f = foldl op* (fromInt 1) o mapi f
 
  (* ニュートン法の収束判定 *)
  fun finNewton eps f x = magnitude (f x) < eps
 
  (* DK法の収束判定 *)
  fun finDK eps f = List.all (finNewton eps f)
 
  (* 多項式の微分係数の近似値 *)
  fun f' c0 j zs zj =
    producti (fn (i, zi) => if i = j then c0 else zj - zi) zs
 
  (* DK法の初期値を求める *)
  fun aberth (c0 :: c1 :: cs) r =
    let 
      val n = length (c1 :: cs)
      (* あまり美しくない *)
      open Real
      fun angle k =
        Math.pi * (2.0 * fromInt k + 0.5) / fromInt n
      open Complex
    in
      List.tabulate (n, fn k =>
        ~c1 / (c0 * fromInt n)
        + fromReal r * exp (i * fromReal (angle k)))
    end
in
  (* ホーナー法による多項式の計算 *)
  fun horner cs z =
    foldl (fn (c, sum) => sum * z + c) (fromInt 0) cs
 
  (* ニュートン法における計算の1ステップ *)
  fun newton1 f f' x =
    x - f x / f' x
 
  (* 指定された精度で収束するまでニュートン法を反復する *)
  fun newton eps f f' =
    iterateWhile (finNewton eps f) (newton1 f f')
 
  (* DK法の1ステップ *)
  fun durandKerner1 (cs as (c0 :: _)) zs =
    mapi (fn (j, zj) => newton1 (horner cs) (f' c0 j zs) zj) zs
  (* 指定された精度で収束するまでDK法を反復する *)
  fun durandKerner eps r cs n =
    iterateWhile (finDK eps (horner cs)) (durandKerner1 cs) n (aberth cs r)
end
 
local
  open Complex
in
  (* x^2-6+9=0の解をニュートン法で求める *)
  val q1Newton =
    newton 1E~6
      (horner (map fromInt [1, ~6, 9]))
      (horner (map fromInt [2, ~6])) 100 (fromReal 0.0)
  (* x^3-1=0の解をニュートン法で求める *)
  val q2Newton =
    map
      (newton 1E~6 
        (horner (map fromInt [1, 0, 0, ~1]))
        (horner (map fromInt [3, 0, 0])) 100)
      [fromReal 0.5, {real = ~1.0, imag = 1.0}, {real = ~1.0, imag = ~1.0}]
  (* x^6-4x^4+2x^2-1=0の解をニュートン法で求める *)
  val q3Newton =
    map 
      (newton 1E~6
        (horner (map fromInt [8, 0, ~4, 0, 2, 0, ~1]))
        (horner (map fromInt [48, 0, ~16, 0, 4, 0])) 1000)
      [fromInt 1,
       fromInt ~1,
       {real = 1.0, imag = 1.0},
       {real = ~1.0, imag = 1.0},
       {real = ~1.0, imag = ~1.0},
       {real = 1.0, imag = ~1.0}]
 
  (* x^2-6+9=0の解をDKA法で求める *)
  val q1DK = durandKerner 1E~6 1.0 (map fromInt [1, ~6, 9]) 100
  (* x^3-1=0の解をDKA法で求める *)
  val q2DK = durandKerner 1E~6 1.0 (map fromInt [1, 0, 0, ~1]) 100
  (* x^6-4x^4+2x^2-1=0の解をDKA法で求める *)
  val q3DK = durandKerner 1E~6 10.0 (map fromInt [8, 0, ~4, 0, 2, 0, ~1]) 100
end

ここでは複素数計算に適当にググって出てきたhttp://www.cs.unm.edu/~rlpm/code/のライブラリを使っています。 ただし、このライブラリーは除算がバグっていて以下のような修正をしないと正しい値が求められません。

diff -u Complex/Complex.sml Complex_after/Complex.sml
--- Complex/Complex.sml 2005-12-17 02:00:43.000000000 +0900
+++ Complex_after/Complex.sml   2014-11-09 21:31:34.000000000 +0900
@@ -88,7 +88,7 @@
           fun f x = Real./(x,Real.+(Real.*(a,a),Real.*(b,b)))
       in
           {real=f (Real.+(Real.*(x,a),Real.*(y,b))),
-           imag=f (Real.+(Real.*(y,a),Real.*(x,b)))}
+           imag=f (Real.-(Real.*(y,a),Real.*(x,b)))}
       end

このバグのせいで休日が1日潰れてしまいました。 数値計算のプログラムは同じ様な型の値を大量に扱うせいで型の庇護が受けられない上、値を見ただけではどこがおかしいのか分かりづらいのでデバッグは非常につらいものがあります。

Standard MLを選んだばかりに苦労した形になって後味が悪いですが、 OCamlStandard ML、同じML系の言語にも関わらず構文や命名規則等色々な点が違っていて面白いです。 研究室次第では来年から書く事になるので今から慣れておくのも良さそうですね。