DKA法
をテンプレートにして作成
[
トップ
|
新規
|
一覧
|
検索
|
最終更新
|
ヘルプ
]
開始行:
代数方程式は重根も別々に考えると&ref(rf_dka.eq1.gif,nolin...
ここまでの方法では適切な初期値や範囲を与えることで1つの解...
一方で代数方程式での複数解を求める場合は初期値や探索範囲...
ここでは&ref(rf_dka.eq1.gif,nolink,70%);個の複素数解を一...
まず，代数方程式が複素解を持つとして，未知数を&ref(rf_dka...
#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,n...
同様に&ref(rf_dka.eq10.gif,nolink,70%);のとき左辺第&ref(r...
#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,no...
#ref(rf_dka.eq15.gif,nolink,70%)
ここで，&ref(rf_dka.eq13.gif,nolink,70%);．
この式をワイヤストラス(Weierstrass)法，もしくは，デュラン...
((http://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_meth...
以下では上式をDK式と呼ぶことにする．
DK式により&ref(rf_dka.eq1.gif,nolink,70%);個の複素数解が...
ニュートン法を用いているため&ref(rf_dka.eq14.gif,nolink,7...
そこで，以下のアバース(Aberth)の初期値を用いる．
#ref(rf_dka.eq17.gif,nolink,70%)
ここで，&ref(rf_dka.eq18.gif,nolink,70%);，&ref(rf_dka.eq...
また，ここでの&ref(rf_dka.eq20.gif,nolink,70%);は&ref(rf_...
Aberthの初期値を用い，DK式を反復計算することで，
&ref(rf_dka.eq1.gif,nolink,70%);個の根を得る方法をDKA(Dur...
***DKA法の手順 [#t0cedca7]
+中心&ref(rf_dka.eq22.gif,nolink,70%);と半径&ref(rf_dka.e...
+DK式で&ref(rf_dka.eq14.gif,nolink,70%);を更新
収束するまで2を繰り返すことで，近似解を得る．
Aberthの初期値の半径を求める方法は様々に提案されているが...
Aberthの方法((O. Aberth, ``Iteration methods for finding ...
&ref(rf_dka.eq23.gif,nolink,70%);として以下の方程式を得る．
#ref(rf_dka.eq24.gif,nolink,70%)
ここで，&ref(rf_dka.eq25.gif,nolink,70%);はホーナーの方法...
&ref(rf_dka.eq26.gif,nolink,70%); でかつ，すべての係数&re...
#ref(rf_dka.eq28.gif,nolink,70%)
の代数方程式を解いた時の解&ref(rf_dka.eq29.gif,nolink,70%...
よって，&ref(rf_dka.eq30.gif,nolink,70%);の解は，&ref(rf_...
&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...
{
// 半径算出のための方程式の係数
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;
// 多項式S(w)の係数
vector<double> b(cd.size());
b[0] = abs(cd[0]);
for(int i = 1; i <= n; ++i){
b[i] = -abs(cd[i]);
}
// Aberthの初期値の半径をニュートン法で算出
int m = 0; // 係数cの中で0でないものの個数
for(int i = 0; i <= n; ++i) m += (fabs(c[i].real()) >...
double rmax = 0.0; // 半径の最大値
for(int i = 1; i <= n; ++i){
double ri = pow(m*fabs(c[i].real())/fabs(c[0].rea...
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%);の...
#code(C){{
/*!
* ホーナー法(組立除法)
- P(x) = a0 x^n + a1 x^(n-1) + ... + a_(n-1) x + a_n ...
- 商は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, vec...
{
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<complex...
{
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.eq47.gif,nolink,70%)
となる．オイラーの公式(&ref(rf_dka.eq48.gif,nolink,70%);)...
#ref(rf_dka.eq49.gif,nolink,70%)
が得られる．ただし，ここでの&ref(rf_dka.eq50.gif,nolink,7...
また，&ref(rf_dka.eq51.gif,nolink,70%);を角度に加えたのは...
終了行:
代数方程式は重根も別々に考えると&ref(rf_dka.eq1.gif,nolin...
ここまでの方法では適切な初期値や範囲を与えることで1つの解...
一方で代数方程式での複数解を求める場合は初期値や探索範囲...
ここでは&ref(rf_dka.eq1.gif,nolink,70%);個の複素数解を一...
まず，代数方程式が複素解を持つとして，未知数を&ref(rf_dka...
#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,n...
同様に&ref(rf_dka.eq10.gif,nolink,70%);のとき左辺第&ref(r...
#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,no...
#ref(rf_dka.eq15.gif,nolink,70%)
ここで，&ref(rf_dka.eq13.gif,nolink,70%);．
この式をワイヤストラス(Weierstrass)法，もしくは，デュラン...
((http://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_meth...
以下では上式をDK式と呼ぶことにする．
DK式により&ref(rf_dka.eq1.gif,nolink,70%);個の複素数解が...
ニュートン法を用いているため&ref(rf_dka.eq14.gif,nolink,7...
そこで，以下のアバース(Aberth)の初期値を用いる．
#ref(rf_dka.eq17.gif,nolink,70%)
ここで，&ref(rf_dka.eq18.gif,nolink,70%);，&ref(rf_dka.eq...
また，ここでの&ref(rf_dka.eq20.gif,nolink,70%);は&ref(rf_...
Aberthの初期値を用い，DK式を反復計算することで，
&ref(rf_dka.eq1.gif,nolink,70%);個の根を得る方法をDKA(Dur...
***DKA法の手順 [#t0cedca7]
+中心&ref(rf_dka.eq22.gif,nolink,70%);と半径&ref(rf_dka.e...
+DK式で&ref(rf_dka.eq14.gif,nolink,70%);を更新
収束するまで2を繰り返すことで，近似解を得る．
Aberthの初期値の半径を求める方法は様々に提案されているが...
Aberthの方法((O. Aberth, ``Iteration methods for finding ...
&ref(rf_dka.eq23.gif,nolink,70%);として以下の方程式を得る．
#ref(rf_dka.eq24.gif,nolink,70%)
ここで，&ref(rf_dka.eq25.gif,nolink,70%);はホーナーの方法...
&ref(rf_dka.eq26.gif,nolink,70%); でかつ，すべての係数&re...
#ref(rf_dka.eq28.gif,nolink,70%)
の代数方程式を解いた時の解&ref(rf_dka.eq29.gif,nolink,70%...
よって，&ref(rf_dka.eq30.gif,nolink,70%);の解は，&ref(rf_...
&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...
{
// 半径算出のための方程式の係数
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;
// 多項式S(w)の係数
vector<double> b(cd.size());
b[0] = abs(cd[0]);
for(int i = 1; i <= n; ++i){
b[i] = -abs(cd[i]);
}
// Aberthの初期値の半径をニュートン法で算出
int m = 0; // 係数cの中で0でないものの個数
for(int i = 0; i <= n; ++i) m += (fabs(c[i].real()) >...
double rmax = 0.0; // 半径の最大値
for(int i = 1; i <= n; ++i){
double ri = pow(m*fabs(c[i].real())/fabs(c[0].rea...
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%);の...
#code(C){{
/*!
* ホーナー法(組立除法)
- P(x) = a0 x^n + a1 x^(n-1) + ... + a_(n-1) x + a_n ...
- 商は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, vec...
{
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<complex...
{
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.eq47.gif,nolink,70%)
となる．オイラーの公式(&ref(rf_dka.eq48.gif,nolink,70%);)...
#ref(rf_dka.eq49.gif,nolink,70%)
が得られる．ただし，ここでの&ref(rf_dka.eq50.gif,nolink,7...
また，&ref(rf_dka.eq51.gif,nolink,70%);を角度に加えたのは...
ページ名:
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
