数値コンピューティングのレポート進捗

numeric.sml

(* 添字も渡す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
      mapi
        (fn (k, _) =>
          ~c1 / (c0 * fromInt n)
          + fromReal r * exp (i * fromReal (angle k)))
        (c1 :: cs)
    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
  val q1Newton =
    newton 1E~6
      (horner (map fromInt [1, ~6, 9]))
      (horner (map fromInt [2, ~6])) 100 (fromReal 0.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}]
  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}]

  val q1DK = durandKerner 1E~6 1.0 (map fromInt [1, ~6, 9]) 100
  val q2DK = durandKerner 1E~6 1.0 (map fromInt [1, 0, 0, ~1]) 100
  val q3DK = durandKerner 1E~6 10.0 (map fromInt [8, 0, ~4, 0, 2, 0, ~1]) 100
end

patch

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