連立方程式計算プログラム、その4。
2007/09/24
JavaScriptで連立方程式を計算する試み第4弾。
今回は行列的な方法をつかわず、人間が手計算する際の計算方法を用いる。行列を使う場合、計算量は指数関数的に上がっていき、わずか9次で手一杯の状況だった(第3弾)。今回の方式では次数が増えても多項式的にしか計算量が増さないため、これまでより遥かに高速に計算できる。連立方程式を解く目的では、恐らくこれが最後になるだろう。
2007/09/24-2009/12/24
人間が手計算する際に行うのと同様に、入力された方程式から未知数を1つずつ消去して最初の未知数を求め、今度は逆に計算を遡って全ての未知数を求める。
計算量を配列の操作回数で求めてみる。まず、n元連立方程式単独では、n(n+1)個の配列要素を消費する。このn個の方程式を操作して、未知数を1つ減らしたn-1元連立方程式を作る。この作業を繰り返し、最終的に1元まで還元して未知数を求める。ここまでの配列操作回数は1〜n元までの和となり、その総計はn(n+1)(n+2)/3となる。逆行列を求める前段階で消費する配列要素数n!に比べれば、圧倒的に低コストである。
(2009/12/24追記)詳しい計算方法を具体的に書いておく。一重に未知数を1つずつ削除すると言ってもわかりづらいようだ。
3元1次連立方程式を例に計算過程を示す。先ず、解く3元連立方程式が次のようなものであるとする。
式を簡単に表現するために行列を用いて表す。
この行列を次のように左側を単位行列となるように式変形することができれば、右側のA,B,Cが求めたいa,b,cになる。
算法の都合上、さらに式を簡単にするために右側1列を合わせた4行3列の行列として書き表しておく。
この行列のi行j列目の成分を、Ci,jと表記する。例えばC0,0=2、C1,3=19である。またi行全体はCiとする。
前述のように、左側3×3を単位行列に書き換える。C0,0が1になるように、C0全体をC0,0で割る。即ちC0=C0/C0,0を実行する。
次にC1,0とC2,0を0にしたい。C0,0は1なので、CiをC0にCi,0を掛けたもので引く。即ち、Ci=Ci-C0×Ci,0とする。もちろん、0行自身は除く。
以上の操作を、全ての行で行えばいい。次はC1,1=1となるようにC1=C1/C1,1とする。
全く同様に次にC0,1とC2,1を0にする。Ci=Ci-C1×Ci,1を実行。
3行目に対しても全く同じ操作を行う。
よって、解(a,b,c) = (4,-1,3)が得られた。
これをn行n列に一般化したものをJavaScriptとして記述すると次のようになる。ただし入力は前述のCi,jに対応する2次配列c
、出力は配列で解が不定解のときisBad
がtrueになる。解自体はans
に格納されている。
cal = function( c ) { var n = c.length; var ans = new Array( n ); var isBad = false; var i,j,k,div1,div2; for( i=0; i<n; i++ ) { div1 = c[i][i]; if(!div1) { isBad = true; break; } for( j=i; j<=n; j++ )c[i][j] /= div1; for( j=0; j<n; j++ ) if(i!=j) { div2 = c[j][i]; for( k=i; k<=n; k++ )c[j][k] -= div2*c[i][k]; } } for( i=0; i<n; i++ )ans[i] = c[i][n]; return [ans,isBad]; }
2007/09/24-2009/12/24
百規模のフォームは余りに巨大すぎるため、表示は省スペース化されている。行(横)が係数、列(縦)が式である。値は2〜300まで指定可能。環境によっては、20〜30を超えた辺りから画面外にはみ出してしまう。また、百規模のテーブルはセル数が1万に達し、計算量は少なくても出力されるHTMLがかなり大きなものになる。
問題点として、手元のIE6が強制終了することがあったが、計算途中にわざとalertを挟むことで、(手元では)なんとか回避できている。他の環境や、そもそも本当に回避できているかは不明なので、予め警告しておく。
(2009/12/24)解の計算方法を効率化し、解が定まらないときはNaN出力ではなく--となり不定解を検出できるようになった。