『Maple における DK 法の実装』
(代数方程式の数値解法)

代数方程式の根を求める数値解法には様々な手法がありますが、もっとも有名な Newton 法をベースにした Duran-Kerner 法を数式処理システム Maple で実装するプログラミング・サンプルを紹介します。
Newton 法は与えられた初期値を元にして反復によりひとつの根を求めますが、DK 法はすべての根を一斉に求めることができます。

>   

 Newton 法と差積

関数 f(x) の根 α に対して、Newton 法は以下の定義で与えられます。

   x[n+1] = x[n]-f(x[n])/diff(f(x[n]),x)

ここで x[n]  は、ある根αの n 回反復により得られた近似値です。

Newton 法では、右辺2項目の分母で示しているように、微分値を計算する必要がありますが、この部分を以下のような 差積 で置き換えることができます。いま根の初期値を alpha[i]  とすると、上記微分値は

Product(alpha[i]-alpha[j],j = 1 .. d)

と表すことができます。ただし、j ≠ i です。

この差積を計算する関数を Maple で定義しておきます。

>    restart:

>    # compute difference product w.r.t i-th value
diffprod := proc(v::list, i::integer)
  local tmpv;

  tmpv := subsop(i=NULL, v); # drop i-th element
  product(v[i]-tmpv[j], j=1..nops(tmpv));
end proc:

試してみます;

>    diffprod([v1, v2, v3, v4, v5], 2); # w.r.t the 2nd value

(v2-v1)*(v2-v3)*(v2-v4)*(v2-v5)

DK 法とは、Newton 法の反復式の中で微分値の計算を行う部分を上記の差積で置き換えた算法です。微分値を差積で計算させることにより、どの根も同じ程度の誤差値を持つようになり、解法の暴走などが比較的防げます。

 DK 法の実装

それでは、DK 法を実装していきます。

まず、DK 法の1ステップあたりの計算を行うプロシージャを定義しておきましょう。引数には、関数、変数名、値のリスト、を取るようにします。

>    oneStepDK := proc(p::algebraic, var::name, v::list)
  [seq(
     v[k] - eval(p, var=v[k])/diffprod(v, k),
     k=1..nops(v)
   )];
end proc:

正しく動作するか試してみます。

>    oneStepDK((x-1.02)*(x-2.34)*(x-3.22), x, [0.2, 1.1, 5.2]);

[1.377665778, 1.156992954, 4.045341268]

計算が行えました。

それでは、最終的な DK 法反復を行うプロシージャを定義します。プロシージャの引数には、多項式、変数名、初期値リストの各情報に加えて、絶対許容誤差の最大値と最大反復回数を指定できるようにしておきます。

>    # define the DK-method for root-finding
solveDK := proc(p::algebraic, var::name, vals::list,
                err::numeric, maxite::posint)
           # err : upper limit of maxmimum absolute error
           # maxite : max number of iteration

  option `Copyright (c) 2003, Cybernet Systems,Co.,Ltd. All rights reserved.`;
  local dkPrev, dkCur, maxerr, n;

  # initial parameters
  dkPrev := oneStepDK(p, var, vals);
  n := 0:
  maxerr := err*2:

  while ((maxerr > err) and (n < maxite)) do
    dkCur := oneStepDK(p, var, dkPrev);
    maxerr := max(op(
                map(abs, dkCur - dkPrev)
              ));
    dkPrev := dkCur;
    n := n+1;
  end do:
  return(dkCur, n, maxerr);
end proc:

>   

 DK 法で計算してみよう

このワークシートで定義した DK 法を用いて、さっそく根の探索を行ってみましょう。

>   

 Case 1:

以下の5次多項式を定義します。

>    poly1 := expand(product(x-sin(i), i=1..5));

poly1 := sin(1)*x*sin(3)*sin(4)*sin(5)+sin(1)*sin(2)*x*sin(4)*sin(5)+sin(1)*sin(2)*sin(3)*x*sin(5)+sin(1)*sin(2)*sin(3)*sin(4)*x-sin(1)*sin(2)*sin(3)*sin(4)*sin(5)+x*sin(2)*sin(3)*sin(4)*sin(5)+sin(1)*...
poly1 := sin(1)*x*sin(3)*sin(4)*sin(5)+sin(1)*sin(2)*x*sin(4)*sin(5)+sin(1)*sin(2)*sin(3)*x*sin(5)+sin(1)*sin(2)*sin(3)*sin(4)*x-sin(1)*sin(2)*sin(3)*sin(4)*sin(5)+x*sin(2)*sin(3)*sin(4)*sin(5)+sin(1)*...
poly1 := sin(1)*x*sin(3)*sin(4)*sin(5)+sin(1)*sin(2)*x*sin(4)*sin(5)+sin(1)*sin(2)*sin(3)*x*sin(5)+sin(1)*sin(2)*sin(3)*sin(4)*x-sin(1)*sin(2)*sin(3)*sin(4)*sin(5)+x*sin(2)*sin(3)*sin(4)*sin(5)+sin(1)*...
poly1 := sin(1)*x*sin(3)*sin(4)*sin(5)+sin(1)*sin(2)*x*sin(4)*sin(5)+sin(1)*sin(2)*sin(3)*x*sin(5)+sin(1)*sin(2)*sin(3)*sin(4)*x-sin(1)*sin(2)*sin(3)*sin(4)*sin(5)+x*sin(2)*sin(3)*sin(4)*sin(5)+sin(1)*...
poly1 := sin(1)*x*sin(3)*sin(4)*sin(5)+sin(1)*sin(2)*x*sin(4)*sin(5)+sin(1)*sin(2)*sin(3)*x*sin(5)+sin(1)*sin(2)*sin(3)*sin(4)*x-sin(1)*sin(2)*sin(3)*sin(4)*sin(5)+x*sin(2)*sin(3)*sin(4)*sin(5)+sin(1)*...

この関数の根を DK 法で計算してみましょう。計算を行う際には方程式を近似化(数値化)しています。

>    solveDK(evalf(poly1), x, [1, 2, 3, 4, 5], 2e-8, 200);

[.9092974329, -.9589242756, -.7568024942, .1411200082, .8414709790], 59, .22e-8

戻り値の最初のリストが得られた近似根、2番目と3番目の戻り値は、与えた許容誤差内で反復した回数と反復を抜けたときの最大誤差値です。

DK 法で計算した上記の戻り値がどの程度の誤差を持っているのか、厳密に計算した根を数値化して比較してみましょう。

>    evalf( solve(poly1=0, x) );

.8414709848, .141120008, -.95892427, .9092974254, -.7568024940

初期値を適当に与えているわりには、最大許容誤差値が示しているとおりの近似値を得ていることがわかります。

>   

 Case 2:

より高精度で計算させることも可能です。

精度(桁数)をあげるには、Maple 内で定義されている精度の制御値を変更します。

>    Digits := 30:

以下の多項式に関する零点を求めてみましょう。見ればわかるとおり、この多項式は3つの正の近接根(絶対的に近い距離にある根)と2つの負の根を持っています。

>    poly2 := expand((x-0.256000)*(x-0.256002)*(x-0.255998)*(x+3.200000)*(x+5.128200));

poly2 := x^5+7.560200*x^4-.275318141075035914240000000000+10.210790399996*x^3+3.086660455571687116800000*x-10.982450790432288800*x^2

この多項式の根を DK 法で計算してみましょう。絶対許容誤差を 1e-40, 最大反復回数は 1000 回とします。

>    dksol := solveDK(poly2, x, [1/5, 1/4, 1/3, -2, -3], 1e-40, 1000);

dksol := [.256002000000000000006830689284, .256000000000000000009773303388, .255997999999999999990139401901, -5.12819999999999999999999999998, -3.19999999999999999999999999994], 1000, .13435229554e-19
dksol := [.256002000000000000006830689284, .256000000000000000009773303388, .255997999999999999990139401901, -5.12819999999999999999999999998, -3.19999999999999999999999999994], 1000, .13435229554e-19

今回実装した DK 法では、反復計算中に指定反復回数まで達するか、または指定された許容誤差を満たすようになれば反復を停止し答えを返すようになっています。

上記の計算結果を見ると、反復回数は 1000 (dksol[2])、最大絶対誤差値は 0.13*10^(-19) となっています。つまり 1000 回の反復計算を行ったが絶対誤差は 0.13*10^(-19) 程度までしかあげられなかった、ということになります。

>   

 参考文献

  • 新数学講座13『数値解析』 一松 信 著,朝倉書店,ISBN4-254-11443-5

閉じる