ページ

2011年7月4日月曜日

C言語で高次方程式を解く

Cのプログラム中で高次方程式の数値解を求めようとして、躓いた。 4次方程式までなら解の公式があるのでそれを使えば良いのだが、どうやら5次以上になると近似解法を使うしかないらしい。いろいろと調べていると、この問題は「コンパニオン行列の固有値を求める」問題に帰着できることが分かった。つまり「コンパニオン行列」さえ作ってやれば、任意の固有値問題解法アルゴリズムで解くことが出来る。MATLABやScilabで多項式の根を求める関数も、この方法を採用しているらしい。今回は固有値の計算にGNU Scientific Library(GSL)の行列演算を利用して、実数係数高次方程式の複素解を求めるラッパー関数を作ってみたので、メモしておく。
※とりあえず手っ取り早く使える物を作ったという話なので、数学的な原理や計算精度についての考察は一切していない。悪しからず。

コンパニオン行列
山川栄樹博士(関西大学)の講義資料(PDF)などによると、
方程式
 について、行列

 あるいはの転置行列をコンパニオン行列とよび、この行列の固有値は方程式の解に等しいとのこと。原理はよくわからない(そこまで厳密に数学をやらなくても許されるのが工学屋の良いところ?)。

 GNU Scientific Library (GSL)
科学技術計算用のC言語ライブラリで、行列や複素数の扱いはもちろんのこと、微分方程式の数値解法なんかも入っているらしい。公式サイトからダウンロードしてインストールする。一般的なUnix系OS(mac含む)+gccなら、シンプルにconfigure, make, make installで使えるようになるはず。

プログラム
GSLを用いた固有値の解法については、九州大学金子邦彦研究室のwebページが大変参考になった。関数の仕様は以下の通りとした。

int roots(double *a, int len_a, int order, double *r, int len_r);
*a方程式の係数を低次から順に入力。
len_a配列aの長さ。方程式の次数+1以上でなければならない。
order方程式の次数。
*r解の出力先配列。解は複素数になり、各解の実数部と虚数部を交互に出力。
len_r配列rの長さ。解は複素数になるので、len_aの2倍以上なければならない。

ソースコードはこちら。GSLライブラリを使うので、リンクオプション -lgsl -lgslcblas が必要。

使用例
プログラム実行時の引数として、ゼロ次の係数から順に与える。
方程式
 を解くならば(解は

$ ./roots -5 7 -3 1
roots
1.000000 + 2.000000i
1.000000 + -2.000000i
1.000000 + 0.000000i

といった具合。この例は3次方程式なのであまりご利益はないけれども、もっと高次の方程式も解けるはず。

0 件のコメント: