非線形方程式の根を数値的に求める方法について
線形方程式と非線型方程式 †
についての方程式
の解を求めたい.
もし,
が
のような線形方程式(1次方程式)であるならば,その解は,
と計算できる(ただし,
).
しかし,
が非線形方程式であったならどうだろうか.
たとえば,以下の
次代数方程式(algebraic equation)を考える.
2次方程式(
)ならば解の方程式,3次や4次方程式の場合はカルダノ(Cardano)やフェラリ(Ferrari)の公式で
解析的に解くことができるが,5次以上だと直接的な解法はない.
また,超越方程式(transcendental equation)と呼ばれる無限次代数方程式(
)もある.
たとえば,
は無限べき級数で表されるため,
を含む方程式は超越方程式となる.
超越方程式に関してもほとんどの場合直接的な解法がない.
以下では非線形方程式を数値的に解く方法として,
- 2分法
- ニュートン・ラフソン法
- 多次元のニュートン・ラフソン法
- ホーナー法
- DKA法
について述べる.
2分法 †
関数
を区間
を考えたとき,
と
の符号が異なる場合(
),区間
内に少なくとも1つの解が存在する.
区間
の中点での関数値
を求め,その値の符号が
と同じならば,
解は
内に存在する.
このとき,解の存在する区間は
になっている.
同様にして,
の中点での関数値を使うことで,解の存在区間は
となる.
この作業を繰り返すことで方程式の解を求める方法を2分法(bisection method)と呼ぶ.
 |
図1 2分法 |
図1は2分法で近似解を求める手順を示している.
現在の解の存在区間を
とする.
は
で初期化される.
中点を
とする.
図1では
であるので,解は
に存在するので,
とする.
次の中点を
であり,
図1では
であるので,解は
に存在する.
これを反復して,
と求めていく.
2分法の手順を以下に示す.
,
とする.
- 中点
における関数値
を求める.
の場合
,
の場合
とする.
- 2に戻る
反復を終了するための収束条件は以下となる.
2分法のコード例を以下に示す.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
| |
int bisection(double func(const double), double xl, double xr, double &x, int &max_iter, double &eps)
{
double f = func(xl);
double fmid = func(xr);
if(f*fmid >= 0.0) return 0.0;
int k;
double dx = fabs(xr-xl), xmid;
for(k = 0; k < max_iter; ++k){
xmid = 0.5*(xl+xr); dx *= 0.5;
fmid = func(xmid);
if(dx < eps || fmid == 0.0){
x = xmid;
max_iter = k; eps = dx;
return 1;
}
if(f*fmid < 0){
xr = xmid;
}
else{
xl = xmid;
f = fmid;
}
}
max_iter = k; eps = dx;
return 0;
}
|
ニュートン・ラフソン法 †
2分法は指定した区間内の解は必ず求めるものの,解の収束は遅い.
関数
が連続で微分可能関数であるならば,関数の勾配値を用いることで収束を早めることができる.
解の初期値を
とする.
での関数
の接線は,
となる.この接線が
軸と交差する(
)点は,
を次の近似値として同様の処理を繰り返すことで解に近づいていく.図2はこの手順を示したものである.
この方法をニュートン・ラフソン法(Newton-Raphson method),もしくは単にニュートン法(Newton method)と呼ぶ.
ニュートン法の式は以下となる.
反復回数を多くすると限りなく
の解に近づく.
反復を終了するための収束条件は,
となる.
ニュートン法は2分法と異なり解のある区間でなく初期近似値1つだけを指定するだけでよく,収束も速い.
しかし,
だけでなくその微分値
も必要で,初期値によっては解へ収束しない場合もある.
 |
図2 ニュートン法 |
ニュートン法のコード例を以下に示す.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
| |
int newton(double func(const double), double dfunc(const double),
double &x, int &max_iter, double &eps)
{
double f, df, dx;
int k;
for(k = 0; k < max_iter; ++k){
f = func(x);
df = dfunc(x);
x = xn-f/df;
dx = fabs(f/df);
if(dx < eps || fabs(f) < eps){
max_iter = k; eps = dx;
return 1;
}
}
max_iter = k; eps = dx;
return 0;
}
|
多次元のニュートン・ラフソン法 †
ニュートン法を連立非線形方程式に一般化する.
ベクトル表記では,
ここで,
である.
ステップ目の近似値を
とし,
の周りで上式をテイラー展開する.
ここで,
は
を
の要素とするヤコビ行列である.
2次以上の項を無視すると,連立非線形方程式は以下となる.
とすると,
が得られる.
この式は
を未知数とした線形連立方程式であり,
LU分解などで解くことで,
が得られる.
そして,以下の式で
を計算する.
例)2元連立非線形方程式の場合 †
そして,
である.よって,
に関する式は以下となる.
について解くと,
これらを用いて
を更新する.
ホーナー法(組立除法) †
2分法やニュートン法を使って以下の
次の代数方程式を解くときの演算回数を考える.
の値を計算する際の乗算回数は
回で
加算が
回である.
乗算回数が
の2乗で増えるため,たとえば,
では55回だが,
だと5050回にもなる.
乗算回数を少なくするために代数方程式を以下のように変形する.
この場合,一番内側の括弧内から計算していくと,
となり,最終的に
となる.
このときの乗算回数は
回,加算も
回である.
特に
が大きな時に計算回数を大幅に減らすことができる.
たとえば,乗算回数は
で10回,
でも100回で済む.
この計算方法はホーナー(Horner)法と呼ばれる.
ホーナー法で代数方程式の値とその導関数を求めるコード例を以下に示す.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
| |
template<class T>
inline T func_h(double x, const vector<T> &b, int n)
{
T f = b[0];
for(int i = 1; i <= n; ++i){
f = b[i]+f*x;
}
return f;
}
template<class T>
inline T dfunc_h(double x, const vector<T> &b, int n)
{
T df = n*b[0];
for(int i = 1; i <= n-1; ++i){
df = (n-i)*b[i]+df*x;
}
return df;
}
|
テンプレート関数にしているのはDKA法で複素数を扱う関係上,様々な型に対応できるようにしたいためである.
DKA法 †
代数方程式は重根も別々に考えると
個の複素数解を持つ(ここでは実数も複素数の一部としている).
ここまでの方法では適切な初期値や範囲を与えることで1つの解が得られた.
一方で代数方程式での複数解を求める場合は初期値や探索範囲を適切に設定し直さなければならない.
ここでは
個の複素数解を一度に計算する方法について述べる.
まず,代数方程式が複素解を持つとして,未知数を
で表すことにする.
代数方程式におけるニュートン法を考える.
上記の代数方程式の解を
とすると,代数方程式は以下のように変形できる.
これを微分すると,
が解
に等しいとき,
であるので,上式の左辺は第1項のみ残る.
同様に
のとき左辺第
項のみ残る.これを式で表すと,
ここで,
である.
を
のそれぞれの解の近似値とすると,ニュートン法の公式から,
ここで,
.
この式をワイヤストラス(Weierstrass)法,もしくは,デュラン・ケルナー(Durand-Kerner)法と呼ぶ
*1.
以下では上式をDK式と呼ぶことにする.
DK式により
個の複素数解が計算できるが,
ニュートン法を用いているため
の初期値
を適切に選ぶ必要がある.
そこで,以下のアバース(Aberth)の初期値を用いる.
ここで,
,
は複素平面上で
個の解を包むような円の半径である.
また,ここでの
は
を示す.
Aberthの初期値を用い,DK式を反復計算することで,
個の根を得る方法をDKA(Durand-Kerner-Aberth)法と呼ぶ.
DKA法の手順 †
- 中心
と半径
を計算し,初期値
を計算
- DK式で
を更新
収束するまで2を繰り返すことで,近似解を得る.
Aberthの初期値の半径を求める方法は様々に提案されているが,ここでは
Aberthの方法*2を説明する.
として以下の方程式を得る.
ここで,
はホーナーの方法を応用した組立除法で得られた係数である.
でかつ,すべての係数
が非ゼロとしたとき,
の解は,
の代数方程式を解いた時の解
を半径とした円内に存在する.
よって,
の解は,
内にある.この
を初期値を算出するための半径として用いる.
は単一の実数解を持つのでニュートン法を使って解ける.
DKA法を使って複数の根を一度に計算するコード例を以下に示す.
まず,Aberthの初期値を求める関数である.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
| |
int aberth(vector<complexf> &z, const vector<complexf> &c, int n, int max_iter, double eps)
{
vector<complexf> cd(n+1, 1.0); 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();
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(cd.size());
b[0] = abs(cd[0]);
for(int i = 1; i <= n; ++i){
b[i] = -abs(cd[i]);
}
int m = 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);
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;
半径を算出するために,ニュートン法のコードを変更し,
係数列を与えると関数値と導関数値をホーナー法で計算するようにした.
また,半径算出のための
の係数算出に使っているhorner関数は以下.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
| |
inline void horner(const vector<double> &a, double b, vector<double> &c, double &rm, int n)
{
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];
}
}
|
初期値を算出後,以下のワイヤストラス法の関数に多項式の係数とともに渡すことで解を得る.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
| |
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;
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の初期値について †
代数方程式の解と係数の関係を導いてみる.
まず,2次方程式の場合を考える.
解
を使うと方程式は以下のように書き直せる.
元々の式と係数を比べることで,以下の関係式が得られる.
同様にしてn次方程式の場合,
よって,n次方程式の解と係数の関係式は,
ここで,
である.
たとえば,3次方程式(
)の場合,
のときの関係式から,
Aberthの初期値では解
の複素平面における重心
を中心とした半径
の円周上に等間隔に
を配置する.つまり,
となる.オイラーの公式(
)から,
が得られる.ただし,ここでの
は
を示す.
また,
を角度に加えたのは,初期値
が実軸に対称に分布することを防ぐためである(共役複素数になることを防ぐ).
ソースコード †
上記のコードを含むVisual Studio 2010用のソースコードを以下に置く.
参考文献 †