![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
#author("2019-10-17T16:28:45+09:00","default:pbcglab_user","pbcglab_user") #author("2019-10-17T18:10:12+09:00","default:pbcglab_user","pbcglab_user") 代数方程式は重根も別々に考えると&ref(rf_dka.eq1.gif,nolink,70%);個の複素数解を持つ(ここでは実数も複素数の一部としている)． ここまでの方法では適切な初期値や範囲を与えることで1つの解が得られた． 一方で代数方程式での複数解を求める場合は初期値や探索範囲を適切に設定し直さなければならない． ここでは&ref(rf_dka.eq1.gif,nolink,70%);個の複素数解を一度に計算する方法について述べる． まず，代数方程式が複素解を持つとして，未知数を&ref(rf_dka.eq2.gif,nolink,70%);で表すことにする． #ref(rf_dka.eq3.gif,nolink,70%) 代数方程式におけるニュートン法を考える． 上記の代数方程式の解を&ref(rf_dka.eq4.gif,nolink,70%);とすると，代数方程式は以下のように変形できる． #ref(rf_dka.eq5.gif,nolink,70%) これを微分すると， #ref(rf_dka.eq6.gif,nolink,70%) &ref(rf_dka.eq7.gif,nolink,70%);が解&ref(rf_dka.eq8.gif,nolink,70%);に等しいとき，&ref(rf_dka.eq9.gif,nolink,70%);であるので，上式の左辺は第1項のみ残る． 同様に&ref(rf_dka.eq10.gif,nolink,70%);のとき左辺第&ref(rf_dka.eq11.gif,nolink,70%);項のみ残る．これを式で表すと， #ref(rf_dka.eq12.gif,nolink,70%) ここで，&ref(rf_dka.eq13.gif,nolink,70%);である． &ref(rf_dka.eq14.gif,nolink,70%);を&ref(rf_dka.eq1.gif,nolink,70%);のそれぞれの解の近似値とすると，ニュートン法の公式から， #ref(rf_dka.eq15.gif,nolink,70%) ここで，&ref(rf_dka.eq13.gif,nolink,70%);． この式をワイヤストラス(Weierstrass)法，もしくは，デュラン・ケルナー(Durand-Kerner)法と呼ぶ ((http://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method))． 以下では上式をDK式と呼ぶことにする． DK式により&ref(rf_dka.eq1.gif,nolink,70%);個の複素数解が計算できるが， ニュートン法を用いているため&ref(rf_dka.eq14.gif,nolink,70%);の初期値&ref(rf_dka.eq16.gif,nolink,70%);を適切に選ぶ必要がある． そこで，以下のアバース(Aberth)の初期値を用いる． #ref(rf_dka.eq17.gif,nolink,70%) ここで，&ref(rf_dka.eq18.gif,nolink,70%);，&ref(rf_dka.eq19.gif,nolink,70%);は複素平面上で&ref(rf_dka.eq1.gif,nolink,70%);個の解を包むような円の半径である． また，ここでの&ref(rf_dka.eq20.gif,nolink,70%);は&ref(rf_dka.eq21.gif,nolink,70%);を示す． Aberthの初期値を用い，DK式を反復計算することで， &ref(rf_dka.eq1.gif,nolink,70%);個の根を得る方法をDKA(Durand-Kerner-Aberth)法と呼ぶ． ***DKA法の手順 [#t0cedca7] +中心&ref(rf_dka.eq22.gif,nolink,70%);と半径&ref(rf_dka.eq19.gif,nolink,70%);を計算し，初期値&ref(rf_dka.eq16.gif,nolink,70%);を計算 +DK式で&ref(rf_dka.eq14.gif,nolink,70%);を更新 収束するまで2を繰り返すことで，近似解を得る． Aberthの初期値の半径を求める方法は様々に提案されているが，ここでは Aberthの方法((O. Aberth, ``Iteration methods for finding all zeros of a polynomial simultaneously'', Math. Comput. 27, pp.339-344, 1973. ))を説明する． &ref(rf_dka.eq23.gif,nolink,70%);として以下の方程式を得る． #ref(rf_dka.eq24.gif,nolink,70%) ここで，&ref(rf_dka.eq25.gif,nolink,70%);はホーナーの方法で得られた係数である． ここで，&ref(rf_dka.eq25.gif,nolink,70%);はホーナーの方法を応用した組立除法で得られた係数である． &ref(rf_dka.eq26.gif,nolink,70%); でかつ，すべての係数&ref(rf_dka.eq25.gif,nolink,70%);が非ゼロとしたとき，&ref(rf_dka.eq27.gif,nolink,70%);の解は， #ref(rf_dka.eq28.gif,nolink,70%) の代数方程式を解いた時の解&ref(rf_dka.eq29.gif,nolink,70%);を半径とした円内に存在する． よって，&ref(rf_dka.eq30.gif,nolink,70%);の解は，&ref(rf_dka.eq31.gif,nolink,70%);内にある．この&ref(rf_dka.eq19.gif,nolink,70%);を初期値を算出するための半径として用いる． &ref(rf_dka.eq32.gif,nolink,70%);は単一の実数解を持つのでニュートン法を使って解ける． DKA法を使って複数の根を一度に計算するコード例を以下に示す． まず，Aberthの初期値を求める関数である． #code(C){{ /*! * Aberthの方法で初期値を算出 * @param[out] z 解探索のための初期値 * @param[in] c 多項式の係数 * @param[in] n 方程式の次数 * @param[in] max_iter 半径計算のための最大反復数 * @param[in] eps 半径計算のための許容誤差 * @return 1:成功,0:失敗 */ int aberth(vector<complexf> &z, const vector<complexf> &c, int n, int max_iter, double eps) { // 半径算出のための方程式の係数 vector<complexf> a(n+1); complexf zc = -c[1]/(c[0]*(double)n); horner(c, a, n, zc); vector<complexf> cd(n+1, 1.0); // 係数c' double c1n = -c[1].real()/n; cd[0] = c[0]; double rm; // 組立除法での余り vector<double> a(n+1), tmp(n); // 係数格納用の一時的な変数 for(int i = 0; i <= n; ++i) a[i] = c[i].real(); // zの多項式をz+c1/nで割っていくことでwの多項式の係数を求める for(int i = n; i > 1; --i){ horner(a, c1n, tmp, rm, i); cd[i] = rm; a = tmp; } cd[1] = a[1]+c1n; vector<double> b(a.size()); b[0] = abs(a[0]); // 多項式S(w)の係数 vector<double> b(cd.size()); b[0] = abs(cd[0]); for(int i = 1; i <= n; ++i){ b[i] = -abs(a[i]); b[i] = -abs(cd[i]); } // Aberthの初期値の半径をニュートン法で算出 double r = 100.0; int m = 0; // 係数cの中で0でないものの個数 for(int i = 0; i <= n; ++i) m += (fabs(c[i].real()) > 1e-6 ? 1 : 0); double rmax = 0.0; // 半径の最大値 for(int i = 1; i <= n; ++i){ double ri = pow(m*fabs(c[i].real())/fabs(c[0].real()), 1.0/(i+1.0)); if(ri > rmax) rmax = ri; } double r = rmax; newton(b, n, r, max_iter, eps); // Aberthの初期値 complexf zc = -c[1]/(c[0]*(double)n); for(int j = 0; j < n; ++j){ double theta = (2*RX_PI/(double)n)*j+RX_PI/(2.0*n); z[j] = zc+r*complexf(cos(theta), sin(theta)); } return 1; } }} ここでは複素数を扱うために，C++のcomplexを使っている． #include <complex> // 複素数 typedef complex<double> complexf; 半径を算出するために，ニュートン法のコードを変更し， 係数列を与えると関数値と導関数値をホーナー法で計算するようにした． また，半径算出のための&ref(rf_dka.eq33.gif,nolink,70%);の係数算出に使っているhorner関数は以下． #code(C){{ /*! * ホーナー法(組立除法)の係数の値を求める * @param[in] a 代数方程式の係数 * @param[out] b x=x0のときの組立除法の係数 * @param[in] x0 係数を計算するxの値 */ template<class T> inline void horner(const vector<T> &a, vector<T> &b, int n, T x0) /*! * ホーナー法(組立除法) - P(x) = a0 x^n + a1 x^(n-1) + ... + a_(n-1) x + a_n を (x-b) で割ったときの商と余りを返す - 商はn-1次の多項式の係数として返す * @param[in] a 代数方程式の係数 * @param[in] b 割る1次式の係数(x-b) * @param[out] c 商であるn-1次の多項式の係数 * @param[out] rm あまり * @param[in] n 方程式の次数(配列aの大きさはn+1,配列cはn) */ inline void horner(const vector<double> &a, double b, vector<double> &c, double &rm, int n) { if(n <= 2) return; b = a; for(int i = 0; i <= n; ++i){ for(int j = 1; j <= n-i; ++j){ b[j] += x0*b[j-1]; } if(n <= 1) return; rm = a[0]; // 最終的に余りになる c.resize(n); for(int i = 1; i < n+1; ++i){ c[i-1] = rm; rm *= b; rm += a[i]; } } }} 初期値を算出後，以下のワイヤストラス法の関数に多項式の係数とともに渡すことで解を得る． #code(C){{ /*! * ワイヤストラス法(DK公式) * @param[inout] z 初期値位置を受け取り，解を返す * @param[in] c 多項式の係数 * @param[in] n 方程式の次数 * @param[inout] max_iter 最大反復数(反復終了後,実際の反復数を返す) * @param[inout] eps 許容誤差(反復終了後,実際の誤差を返す) * @return 1:成功,0:失敗 */ int weierstrass(vector<complexf> &z, const vector<complexf> &c, int n, int &max_iter, double &eps) { double e = 0.0, ej; vector<complexf> zp; complexf f, df; int k; for(k = 0; k < max_iter; ++k){ zp = z; // DK式の計算 for(int j = 0; j < n; ++j){ f = func(z[j], c, n); df = c[0]; for(int i = 0; i < n; ++i){ if(i != j){ df *= zp[j]-zp[i]; } } z[j] = zp[j]-f/df; } // 誤差の算出 e = 0.0; for(int j = 0; j < n; ++j){ if((ej = abs(func(z[j], c, n))) > e){ e = ej; } } // 収束判定 if(e < eps){ max_iter = k; eps = e; return 1; } } eps = e; return 0; } }} ***Aberthの初期値について [#dccca43a] 代数方程式の解と係数の関係を導いてみる． まず，2次方程式の場合を考える． #ref(rf_dka.eq34.gif,nolink,70%) 解&ref(rf_dka.eq35.gif,nolink,70%);を使うと方程式は以下のように書き直せる． #ref(rf_dka.eq36.gif,nolink,70%) 元々の式と係数を比べることで，以下の関係式が得られる． #ref(rf_dka.eq37.gif,nolink,70%) 同様にしてn次方程式の場合， #ref(rf_dka.eq38.gif,nolink,70%) よって，n次方程式の解と係数の関係式は， #ref(rf_dka.eq39.gif,nolink,70%) ここで，&ref(rf_dka.eq40.gif,nolink,70%);である． たとえば，3次方程式(&ref(rf_dka.eq41.gif,nolink,70%);)の場合， #ref(rf_dka.eq42.gif,nolink,70%) &ref(rf_dka.eq43.gif,nolink,70%);のときの関係式から， #ref(rf_dka.eq44.gif,nolink,70%) Aberthの初期値では解&ref(rf_dka.eq45.gif,nolink,70%);の複素平面における重心 #ref(rf_dka.eq46.gif,nolink,70%) を中心とした半径&ref(rf_dka.eq19.gif,nolink,70%);の円周上に等間隔に&ref(rf_dka.eq16.gif,nolink,70%);を配置する．つまり， #ref(rf_dka.eq47.gif,nolink,70%) となる．オイラーの公式(&ref(rf_dka.eq48.gif,nolink,70%);)から， #ref(rf_dka.eq49.gif,nolink,70%) が得られる．ただし，ここでの&ref(rf_dka.eq50.gif,nolink,70%);は&ref(rf_dka.eq21.gif,nolink,70%);を示す． また，&ref(rf_dka.eq51.gif,nolink,70%);を角度に加えたのは，初期値&ref(rf_dka.eq16.gif,nolink,70%);が実軸に対称に分布することを防ぐためである(共役複素数になることを防ぐ)．
