有理関数型 Brent 法の提案
日大生産工 (院)
日大生産工
1.
○茂木 渉
篠原 正明
逆 2 次補間の公式は,式の次数が 2 と低いことなど
はじめに
非線形方程式の数値解法は数多く存在しているが,根
の理由から,Lagrange 補間 (Lagrange interpolation)
の存在する保証がない,存在するとしても根の個数がわ
によって導出するのが簡単である.Lagrange の補間多
からない,どんな問題にも適用できる万能な解法が存
項式は (1),(2) 式で与えられる.
在しないといった理由から,非線形方程式の根を求める
y=
ことは困難であるとされている.従って,手法の改良,
N
X
組み合わせ,使い分けなど,問題に応じて試行錯誤して
いかなければならない.
N (j6=k)
Lk (x) =
例えば,一般的に非線形方程式を解くには Newton束性がないなどの弱点もある.そこで序盤は二分法な
どで荒い解を求めておき,最終段階で根の精度に磨き
本論文では,導関数や関数形が利用できない方程式
有効な手法といわれている「Brent の方法」に「有理関
+
数補間」を組み合わせた解法アルゴリズムを提案する.
少なくとも 1 次収束させるというものであり,初期区間
が根を含み,その区間内で関数が計算できる限り,必ず
収束することが保証されており,確実性と速さを兼ね備
える優れた手法であるといわれている.
2.1.
逆 2 次補間
x =
+
+
(3)
[y − f (a)][y − f (b)]c
[f (c) − f (a)][f (c) − f (b)]
[y − f (b)][y − f (c)]a
[f (a) − f (b)][f (a) − f (c)]
[y − f (c)][y − f (a)]b
[f (b) − f (c)][f (b) − f (a)]
(4)
この逆 2 次関数と x 軸との交点は,y = 0 と置くこと
により,(5) 式のように書くことができる.このときの
x の値を根の推定値とする.
逆 2 次補間 (inverse quadratic interpolation) とは,
3 点を通る逆 2 次関数 (2 次関数の逆関数) を使い,根の
x =
近くに急速に収束させることをねらう手法である.根
の存在範囲を超えて発散してしまう危険性や二分法よ
りも収束が遅い場合もあるが,妥当な場合には 1.839 次
収束する.
(x − a)(x − b)
f (c)
(c − a)(c − b)
(x − b)(x − c)
f (a)
(a − b)(a − c)
(x − c)(x − a)
f (b)
(b − c)(b − a)
る.[x ↔ y, a ↔ f (a), b ↔ f (b), c ↔ f (c)]
超 1 次収束するはずの方法を監視しながら使用し,超
1 次収束しないようならば一時的に二分法に切り替え,
(2)
(3) 式の逆関数,つまり x と y を入れ替えた関数を作
は,二分法の確実性と,逆 2 次補間や割線法といった高
次の収束をする手法を組み合わせた解法である [1,3,4].
=
+
を問題として想定し,根の囲い込みを前提とする中で
Van Wjingaarden-Dekker-Brent 法 (Brent の方法)
x − xj
xk − xj
に代入すると,(3) 式のようになる.
y
Van Wjingaarden-Dekker-Brent 法
(1)
3 点を [a, f (a)], [b, f (b)], [c, f (c)] として補間式 (1),(2)
をかけるために Newton-Raphson 法を使う,といった
混成アルゴリズムがある.
Y
j=0
Raphson 法を使用するのが良いとされるが,大域的収
2.
Lk (x)yk
k=0
+
+
f (a)f (b)c
[f (c) − f (a)][f (c) − f (b)]
f (b)f (c)a
[f (a) − f (b)][f (a) − f (c)]
f (c)f (a)b
[f (b) − f (c)][f (b) − f (a)]
Proposal of Rational Function Type Brent’s Method
Wataru MOGI† and Masaaki SHINOHARA
(5)
3.
有理関数逆補間
補間に多項式ではなく,有理関数を用いる方法が有理
関数補間 (rational f unction interpolation) である.
多項式補間ではうまく近似できない場合も,有理関
数補間では近似できることがあり,収束次数の面でも有
利であることが示されている [2].
3 点 (xk , yk ), (xk−1 , yk−1 ), (xk−2 , yk−2 ) を基にして有
理関数補間をするには,最も簡単な有理 1 次式
y=
図 1. 逆 2 次補間
αx + β
γx + δ
(7)
を使用する.ここで,δ は任意で決められるため,δ = 1
とするのが一般的である.
2.2.
割線法
3 点を (7) 式に代入し,それらを「α, β, γ を未知数と
割線法 (secant method) は 2 点を通る直線を利用し
て根を推定する手法であり,線形関数で補間をすること
から,線形逆補間法とも呼ばれる.また,割線法の収束
する 3 元連立 1 次方程式」とみなして解き,(7) 式にお
いて y = 0 を満たす x の値を根の推定値とする.
但し,連立方程式をこのまま解くのは手間がかかる
次数は,黄金分割比 (=1.618) であることが知られてい
ので,以下のような工夫をする.
る.但し,確実に根を囲い込む保証はなく,状況によっ
[原点移動] 変数変換


 ξ
ξk−1


ξk−2
ては発散する可能性もある.
2 点を [a, f (a)], [b, f (b)] とすると,補間式は (6) 式の
ようになる.
x = b − f (b)
b−a
f (b) − f (a)
(6)
=
=
=
x − xk
xk−1 − xk
xk−2 − xk
(8)
αξ + β − γξy = y
(9)
と置き,有理 1 次式を
とすれば「x = xk → y = yk 」の条件から,ただちに
β = yk
(10)
が得られるので,あとは残りの 2 元連立 1 次方程式
(
αξk−1 − γξk−1 yk−1 = yk−1 − yk
(11)
αξk−2 − γξk−2 yk−2 = yk−2 − yk
から α を求める.連立 1 次方程式を Cramer の公式
図 2. 割線法
2.3.
二分法
二分法 (bisection method) は古典的で単純であるが,
失敗のない頑健 (robust) な方法である.
ある区間 [xl , xu ] において,f (xl ) と f (xu ) の符号が
(Cramer0 s rule) によって解くと,
y
k−1 − yk −ξk−1 yk−1 yk−2 − yk −ξk−2 yk−2 α =
ξ
k−1 −ξk−1 yk−1 ξk−2 −ξk−2 yk−2 =
ξk−2 yk−2 (yk − yk−1 ) − ξk−1 yk−1 (yk − yk−2 )
ξk−1 ξk−2 (yk−1 − yk−2 )
(12)
異なれば,区間内に根が存在することが中間値の定理に
よって保証される.ここで,区間の中点 xm = (xl +xu )/2
となる.これより (9) 式の有理関数補間式を y = 0 と置
における関数の値 f (xm ) の符号を調べれば,根が中点
いた 1 次方程式
のどちらにあるか絞り込むことができる.
0 = αξ + β
(13)
を解き,ξ を求め,
s := b −
x = xk + ξ
(14)
割線法で根の推定を試み,STEP9 へ.
b−a
s := b − f (b)
f (b) − f (a)
有理関数型 Brent 法
ここで,Brent の方法における逆 2 次補間の代わり
に,有理関数逆補間を用いれば,確実かつ高速なアルゴ
Brent の方法の基礎アルゴリズム [4] を基に作成したア
ルゴリズムを以下に示す.本アルゴリズム中に a, b, c, d, s
及び mf lag という変数を使用しているが,これらはそ
(i) a は f (a) と f (b) の符号が反対になるような
前反復時もしくはそれ以前の値.
(ii) b は最新の値で,f (x) = 0 に最も近い値.
(iii) c は前反復時の値.
(iv) d は前々反復時の値.
(v) s は補間または二分法によって得られた根の
推定値.
(vi) mf lag は二分法を使用したか否かを示すフ
ラグ.
[ST EP 1]
f (a) · f (b) < 0 となるような初期区間 [a, b]
を決定し,STEP2 へ.
[ST EP 2]
もし |f (a)| < |f (b)| ならば STEP3 へ,そ
うでなければ STEP4 へ.
[ST EP 3]
a と b を入れ替え,STEP4 へ.
[ST EP 4]
c := a とし,STEP5 へ.
[ST EP 5]
mf lag = 1 とし,STEP6 へ.
[ST EP 6]
もし,3 点 a, b, c が互いに異なる,つまり
f (a) 6= f (b) かつ f (b) 6= f (c) ならば STEP7(有
理関数逆補間) へ.そうでなければ STEP8(割
線法) へ.
[ST EP 7]
有理関数逆補間で根の推定を試み STEP9 へ.
(16)
[ST EP 9]
リズムができるだろうと考えた.
れぞれ次のような意味を持つ.
(15)
[ST EP 8]
を根の近似値とする.
4.
f (b)[c − b][a − b][f (c) − f (a)]
f (a)[a − b][f (b) − f (c)] − f (c)[c − b][f (b) − f (a)]
Brent の追加条件 1,2,3 のいずれかを満た
してしまうならば,補間は失敗であるため
STEP10(二分法) へ.満たさないのならば
補間値を採用するので STEP12 へ.
[Brent の追加条件]
1. s が (3a + b)/4 と b の間にない
2. mf lag = 1 であり,b 6= c かつ |s−b| ≥
|b − c|/2 である
3. mf lag = 0 であり,c 6= d かつ |s−b| ≥
|c − d|/2 である
[ST EP 10]
補間失敗.二分法で根を推定し,STEP11 へ.
a+b
s :=
(17)
2
[ST EP 11]
mf lag = 1 として,STEP13 へ.
[ST EP 12]
mf lag = 0 として,STEP13 へ.
[ST EP 13]
d := c として,STEP14 へ.
[ST EP 14]
c := b として,STEP15 へ.
[ST EP 15]
もし f (a) · f (b) < 0 ならば STEP16 へ.そ
うでなければ STEP17 へ.
[ST EP 16]
b := s と置き,STEP18 へ.
[ST EP 17]
a := s と置き,STEP18 へ.
[ST EP 18]
もし |f (a)| < |f (b)| ならば STEP19 へ.そ
うでなければ STEP20 へ.
[ST EP 19]
a と b を入れ替え,STEP20 へ.
[ST EP 20]
停止条件を満たしていれば b を返して終了
する.満たしていなければ STEP6 へ.
5.
数値計算例
提案した手法と Brent の方法とで,以下の非線形方
程式の問題を解き,両手法の比較を行った.各例題に記
された括弧内の値は初期区間を意味している.
また,計算の停止条件には,一般的によく使用される
ε 法と δ 法を用い,2 つを併用 (OR の形で使用) した場
合で,停止条件を満たすまでのそれぞれの反復回数及
びそのときの関数値と根の存在区間を比較した.
[ε 法]
関数の絶対値で収束判定をする方法.
|f (b)| < ε
(18)
(19)
但し,ε, δ は収束判定用定数であり,本論文では共に
10−8 としている.
例 1.
2e
x−1
−1 = 0
tanh x + 0.2x + 0.3 = 0
例 3.
x − sin x − cos x = 0
[0, 2] (22)
例 4.
ln x − x + 2 = 0
[2, 4] (23)
(x + 3)(x − 1)2 = 0
7回
2.55 × 10−11
9.28 × 10−5
例2
6回
5.83 × 10−10
7.49 × 10−4
例3
6回
1.95 × 10−9
2.07 × 10−6
例4
4回
1.54 × 10−9
6.83 × 10−6
例5
10 回
4.97 × 10−14
6.04 × 10−8
−10
例6
6回
5.50 × 10
3.43 × 10−4
例7
6回
2.82 × 10−12
7.88 × 10−9
表 2. 数値実験結果 (有理関数型 Brent 法)
回数
|f (b)|
|b − a|
6回
2.21 × 10−9
1.14 × 10−4
例2
6回
−10
3.39 × 10
5.82 × 10−4
例3
6回
2.43 × 10−9
2.04 × 10−6
例4
4回
1.27 × 10−9
5.62 × 10−6
例5
10 回
7.10 × 10−15
9.04 × 10−7
例6
6回
2.82 × 10−11
1.29 × 10−7
例7
5回
5.06 × 10−10
1.32 × 10−5
[−3, 3] (20)
例 2.
例 5.
例1
例1
[δ 法]
根の存在区間の幅で収束判定をする方法.
|b − a| < δ
表 1. 数値実験結果 (Brent の方法)
回数
|f (b)|
|b − a|
[−3, 3] (21)
[−4, 4/3] (24)
例 6.
tan x − 3x + 1 = 0
[0, 1] (25)
例 7.
x3 − 6x2 + 12x − 11 = 0
[3, 4] (26)
さらに,補間式内の演算回数を比較しても,(5) 式は
20 回,(15) 式は 17 回であり,これらのことから有理関
数補間の有効性や Brent の方法への組み込みの有用性
などは十分に示せたのではないかと考えられる.
具体的に実装をする上では,手法の切り替え時の無
駄な計算量の減少,丸め誤差の全体の戦略への影響を
考慮する必要があるため,それらを念頭に入れた有理
関数逆補間公式の導出やアルゴリズムの修正などは今
後の課題とする.
参考文献
6.
おわりに
関数の値のみ利用可能な非線形方程式の解法として
高い評価を得ている Brent の方法の改良案として,有
理関数補間を活用した有理関数型 Brent 法を提案した.
数値計算の比較により,反復回数は Brent の方法より
1 回少なくなるか同等程度であり,多くなるということ
はほぼない.反復回数が減る理由は,逆 2 次補間では
Brent の追加条件を満たしてしまうが,有理関数補間で
は満たさず,補間がうまくいくことによる.
また,同じ条件で逆 2 次補間と有理関数補間を使用
した場合,有理関数補間の方が根に近い値に推定され
ることが多く,同じ反復回数での最終的な数値解を比較
しても,(関数の絶対値や根の存在区間が 0 に近いとい
う意味で) 厳密解に近いことが多い.
[1] 玄光男・井田憲一:
「パソコン数値計算ライブラリ」,
株式会社 HBJ(1986)
[2] 戸川隼人:「UNIX ワークステーションによる科学
技術計算ハンドブック [基礎篇 C 言語版]」,サイエ
ンス社 (1992)
[3] William H.Press・Brian P.Flannery・
Saul A.Teokolsky・William T.Veterling(著),
丹 慶 勝 市・佐 藤 俊 郎・奥 村 晴 彦・小 林 誠 (訳):
「NUMERICAL RECIPES in C[日本語版]」 ,技術
評論社 (1993)
[4] Wikipedia-the free encyclopedia-:
「Brent’s method」,
http://en.wikipedia.org/wiki/Brent’s_method,
最終閲覧日 (2008/6/14)
ダウンロード

有理関数型Brent法の提案