Snap.svgでDurand-Kerner-Aberth(DKA)法を可視化

Durand-Kerner-Aberth(DKA)法とは

多項式の解と係数の関係を使い、複素数を考えることにより多項式の解のすべてを一度に得ようという方法である。その際、Aberthの提案した初期値を使うとうまく収束しやすい。アルゴリズムの提唱者であるDurand、Kernerと合わせてDKA法と言われている。

計算ルーチン

方程式

a_nx^n+a_{n-1}*x^{n-1}+\dots+a_1x+a_0=0

初期値の計算

z_j^{(0)}=-\frac{a_{n-1}}{na_n}+R\exp\left(i\frac{2\pi}{n}\left(j-\frac{3}{4}\right) \right)\\

ここで$i$は虚数単位で$R$はすべての解より大きな定数が良いらしい。$R$は

x^n+|a_{n-2}/a_n|x^n+\dots+|a_1/a_n|x+|a_0/a_n|=0

の解を使うこれはnewton法で解く。

実装

init.js

function init(){
//*****************//
//*** 係数の取得 ***//
//*****************//
    n=$('#degree').val(), arr=[];
    for( let i=0; i<=n; i++ ) arr.push(parseFloat($('#term'+i).val()));
    if( arr[n]===0 ) alert(n+'次の項が0です。');

    beta=-arr[n-1]/(arr[n]*n);
//*** c_i=a_i/a_n ***//
    const arr2=[];
    for( let i=0; i<n-1; i++ ) arr2.push(-mathjs.abs(arr[i]/arr[n]));
    const funcR=(x)=>{
        let val=mathjs.pow(x, n);
        for( let i=0; i<n-1; i++ ) val+=arr2[i]*Math.pow(x, i);
        return val;
    }
//*** newton法 ***//
    R=math.solver.newton(funcR, 100);
//*** 初期値の設定 ***//
    ans0=[], ans1=[];
    for( let i=0; i<n; i++ ){
        const tmp=2*Math.PI*i/n+3/(2*n);
        ans0.push(mathjs.add(beta, mathjs.multiply(R, mathjs.exp(mathjs.multiply(mathjs.sqrt(-1),tmp)))));
        ans1.push(0);
    }
    console.log(ans0);

    console.log('b=', beta, " R=", R);
//*** グラフの描画範囲設定 ***//
    const range=1.1*(Math.abs(beta)+R);
    graph.set_range(-range, range, -range, range);
//*** 初期値の円を描く ***//
    graph.draw_circle(beta, 0, R);
//*** 色の設定 ***//
    color=[];
    for( let i=0; i<n; i++ ) color.push(svg.color.deg(360*i/n));

    next();
}

うん、苦しい。パラメーターとかグラフはどこからでもいじれるようにグローバル変数にした。見た目の実装と計算アルゴリズムが絶妙に混じって見づらい...、それに輪をかけてmathjsの計算部分が...orz。

newton法

newton.js

import { diff } from '../diff.js'

const MAX_ITER=1000;

export default function(func, x0, err){
    if( err===void 0 || err<0 ) err=1.0e-8;

    let x1=x0-func(x0)/diff.central(func, x0);
    for( let ite=0; ite<MAX_ITER; ite++ ){
        x0=x1;
        x1=x0-func(x0)/diff.central(func, x0);
        if( Math.abs(x1-x0)<err ) return x1;
    }
    throw new Error('math.solver.newton '+MAX_ITER+'回で収束しませんでした');
}

diff.js

export const diff={
    central: function(func, x, delta){
        delta= delta || 1.0e-8;
        return (func(x+delta)-func(x-delta))/(2.0*delta);
    }
}

svgの描写用ユーティリティのアップデート

svg_update.js

export function draw_circle(disp, x0, y0, r, N){
    N=N | 100;
    let path_data='', is_over_range=true;
    for( let i=0; i<=N; i++ ){
        const x=x0+r*Math.cos(2*Math.PI*i/N);
        const y=y0+r*Math.sin(2*Math.PI*i/N);

        if( disp.ymin<y && y<disp.ymax && disp.xmin<x && x<disp.xmax ){
            if( is_over_range ) path_data+='M'+disp.convert_x(x)+','+disp.convert_y(y)+' ';
            else path_data+='L'+disp.convert_x(x)+','+disp.convert_y(y)+' ';
            is_over_range=false;
        }
        else is_over_range=true;
    }
    return disp.svg().path(path_data).attr({ fill: 'none', strokeWidth: 1, stroke: 'black' });
}

export function draw_point(disp, x, y){
    return disp.svg().circle(disp.convert_x(x), disp.convert_y(y), 3);
}

export default draw_line(disp, x0, y0, x1, y1){
    return disp.svg().line(disp.convert_x(x0), disp.convert_y(y0), disp.convert_x(x1), disp.convert_y(y1))
        .attr({ stroke: 'black', strokeWidth: 2 });
}

export default svg.color.deg(degree){
    return 'hsb('+degree+'deg,1,1)';
}

それぞれ、displayクラスのメソッドにしてgraphクラスから呼ぶようにしてあるがその部分は略す。displayでアスペクト比が崩れているのでcircleを使わずpathオブジェクトで書いている。
色の設定のためHSB色指定を用いた。

解の改善

z^{(k+1)}_j=z_j^{(k)}-\frac{f(z_j^{(k)})}{a_n \left(\prod_{j\neq l}\left(z^{(k)}_j-z^{(k)}_l\right) \right)}

$\prod$は総乗記号、要は掛け算、$j\neq l$も自分自身の引き算で割ると無限大になるからプログラムに落としこむとforループとif( j==l ) continue分母は現在の値、それを使って最も変化量が大きそうな項で0点を近似する。

実装

next.js

function next(){
    for( let i=0; i<n; i++ ){
        ans1[i]=0;
        for( let j=0; j<=n; j++ ){ ans1[i]=mathjs.add(ans1[i], mathjs.multiply(arr[j], mathjs.pow(ans0[i], j))); }
        ans1[i]=mathjs.divide(ans1[i], arr[n]);
        for( let j=0; j<n; j++ ){
            if( j===i ) continue;
            ans1[i]=mathjs.divide(ans1[i], mathjs.subtract(ans0[i], ans0[j]));
        }
        ans1[i]=mathjs.subtract(ans0[i], ans1[i]);
    }

//*** 前のポイントと線の描写 ***//
    for( let i=0; i<ans1.length; i++ ){
        graph.add_line(ans0[i].re, ans0[i].im, ans1[i].re, ans1[i].im).attr({ stroke: color[i] });
    }
    for( let i=0; i<ans0.length; i++ ){ graph.add_point(ans0[i].re, ans0[i].im).attr({fill: color[i], stroke: color[i]}); }

    const err=[];
    for( let i=0; i<n; i++ ) err.push(mathjs.abs(mathjs.subtract(ans0[i], ans1[i])));
    const thre=get_cut();
    if( err.filter((e)=>{ return e>thre; }).length===0 ){
        for( let i=0; i<ans1.length; i++ ){ graph.add_point(ans1[i].re, ans1[i].im, 4).attr({ fill: color[i], stroke: 'black' }); }
        return true;
    }
//*** 次のポイントの描写 ***//
    for( let i=0; i<ans1.length; i++ ){ graph.add_point(ans1[i].re, ans1[i].im).attr({fill: color[i], stroke: color[i]}); }
    for( let i=0; i<n; i++ ) ans0[i]=ans1[i];
    return false;
}

function get_cut(){
    const val=parseFloat($('#err').val());
    const dig=parseInt($('#err_digit').val());
    return val*Math.pow(10, dig);
}

get_cutは終了閾値を得るための関数。
終了条件の判定は満たすものはfilterで除去してすべて除去されて要素数が0になったら終了とした。
やっぱり、mathjs.multiply/add/subtract/divideあたりはかなりキツイ。この辺は手動デバッグ全開でもともとほぼ一行だったものを展開もした。
また答えの置き換えは配列の置き換えでは死んだので要素ごとの置き換えにした。

html

index.html

<!DOCTYPE html>
<html>
  <head>
    <meta charset="utf-8">
    <link rel="stylesheet" type="text/css" href="css/graph.css">
    <title>Math by Javascript</title>
    <script src="lib/Snap.svg-0.5.1/dist/snap.svg-min.js"></script>
    <script src="node_modules/mathjax/MathJax.js?config=TeX-AMS_CHTML"></script>
    <script type="text/x-mathjax-config">
      MathJax.Hub.Config({
      displayAlign: "left",
      displayIndent: "2em"
      });
    </script>
    <script src="dist/DKA_sample.js"></script>
  </head>
  <body>
    <h1>Durand-Kerner-Aberth(DKA)法で多項方程式を解く</h1>
    <div id="canvas">
    </div>

    <form>
      n : <input id="degree" type="Number" min="0" max="10" step="1" value="5">
      <div id="term">
      </div>
      <div id="equation">
      </div>
      <div id="action">
        <button id="solve" type="button">解く</button>
        <input id="err" type="Number" min="0.1" max="9.9" step="0.1", value="1" size="1">
        e
        <input id="err_digit" type="Number" min="-12" max="0" step="1", value="-8" size="1">
      </div>
    </form>
  </body>
</html>

util.js

$(()=> { //ロード時実行関数
    make_input();
    $('#degree').change(()=>{ make_input() });
    $('#solve').click(()=>{
        if( $('#solve').html()==='解く' ){
            init();
            $('#solve').html('次');
        }
        else if( $('#solve').html()==='次' ){
            if( next() ){
                $('#solve').html('更新');
                confirm('終了しました');
            }
        }
        else if( $('#solve').html()==='更新' ){ window.location.reload(); }
    });

    canvas=svg.make.top('#canvas');
    graph=canvas.graph(0.1, 0.1, 0.95, 0.95);
});

//*** 方程式の入力部分を作る ***//
function make_input(){
    const n=$('#degree').val();
    const term=$('#term').empty();

    for( let i=n; i>=0; i-- ){
        const span=$('<span></span>)').html('\\(a_{'+i+'}\\):');
        span.append($('<input>').attr({ id: 'term'+i, size: '1', value: 1 }));
    term.append(span);
    }
    draw_eq();
    MathJax.Hub.Typeset();
}

//*** 方程式をmathjaxを使って書く ***//
function draw_eq(){
    const n=$('#degree').val(), arr=[];
    for( let i=0; i<=n; i++ ) arr.push(parseFloat($('#term'+i).val()));
    let tex='\\[';
    for( let i=n; i>=1; i-- ) tex+=arr[i].toString()+'x^{'+i+'}+';
    tex+=arr[0].toString()+'=0 \\]';
    $('#equation').html(tex);
}

感想

DKA法が可視化すると面白そうなアルゴリズムだから作ってみた。JavaScriptでの複素数の実装が思ったよりめんどくさかった、一応多少の改善策はある。
数値計算のアルゴリズムの可視化に挑戦したがアルゴリズムに可視化を混ぜるとどっちつかずの実装になってコードが混沌としてくる。そもそも逐次近似をブツブツ切って他の処理を挟むとか実用的ではない。
SVGのユーティリティは今回、pointとlineとcircleを足した。後、簡単にだが色についても使ってみた。Snap.svgにはあまりユーティリティはなく、直接、css要素を作っている感じである。RGB、以外にHSB,HSLなどがある。色による分類はこっちを使ったほうが良さそうである。もう少しユーティリティ化は進めたい。