ブレント法のソースを表示
←
ブレント法
ナビゲーションに移動
検索に移動
あなたには「このページの編集」を行う権限がありません。理由は以下の通りです:
この操作は、次のグループに属する利用者のみが実行できます:
登録利用者
。
このページのソースの閲覧やコピーができます。
'''ブレント法'''（{{lang-en-short|Brent's method}}）は、[[二分法]]、[[割線法]]、[[逆2次補間]]を組み合わせた、複雑ではあるが広く用いられている、[[数値解析]]における[[求根アルゴリズム]]の1つである。二分法の安定さを有し、かつ安定でない他の手法と同程度に高速に解が求められる場合もある。可能な限り、より収束の早い割線法や逆2次補間を用い、必要に応じてより安定な二分法に切り替えて解を求めるという手法である。ブレント法は、[[1969年]]の{{ill|セオドラス・デッカー|en|Theodorus Dekker}}による初期のアルゴリズムを元にして、[[1973年]]に{{ill|リチャード・ブレント|en|Richard P. Brent}}により考案されたものである<ref name="Brent1973">Brent (1973). ''Algorithms for Minimization without Derivatives.'' Prentice-Hall, Englewood Cliffs, NJ. Reprinted by Dover Publications, Mineola, New York, January 2002. ISBN 0-486-41998-3. [http://maths-people.anu.edu.au/~brent/pd/rpb011i.pdf オリジナル版]が[[オーストラリア国立大学]]の著者のホームページより入手可能。</ref> 。 == デッカー法 == 二分法と割線法を組み合わせるという考えは、デッカーによるものである。 ここでは、方程式 {{math|''f''(''x'') {{=}} 0}} の解を求めたいとする。まず、二分法と同様に、{{math|''f''(''a''<sub>0</sub>)}} と {{math|''f''(''b''<sub>0</sub>)}} が互いに逆符号を持つような {{math|''a''<sub>0</sub>}} と {{math|''b''<sub>0</sub>}} の2点を初期値とする。{{math|''f''}} が区間 {{math|[''a''<sub>0</sub>, ''b''<sub>0</sub>]}} で連続であるとき、[[中間値の定理]]により、{{math|''a''<sub>0</sub>}} と {{math|''b''<sub>0</sub>}} の間に解が存在する。 反復計算の {{math|''k''}} ステップ目において、以下の3点が用いられる： * {{math|''b''<sub>''k''</sub>}} は現在の反復での値、つまりその時点で推定される {{math|''f''}} の解。 * {{math|''a''<sub>''k''</sub>}} は "反対点"、つまり {{math|''f''(''a''<sub>''k''</sub>)}} と {{math|''f''(''b''<sub>''k''</sub>)}} が逆符号を持つような点であり、したがって区間 {{math|[''a''<sub>''k''</sub>, ''b''<sub>''k''</sub>]}} に解が含まれる。また、{{math|{{abs|''f''(''b''<sub>''k''</sub>)}}}} は {{math|{{abs|''f''(''a''<sub>''k''</sub>)}}}} と等しいか、またはより小さい値とする。したがって {{math|''b''<sub>''k''</sub>}} は {{math|''a''<sub>''k''</sub>}} よりも求める解に近い。 * {{math|''b''<sub>''k''−1</sub>}} は1つ前の反復での値（最初の反復計算 ({{math|''k'' {{=}} 0}}) では {{math|''b''<sub>−1</sub> {{=}} ''a''<sub>0</sub>}} とする）。 次の反復値を求めるため、2つの値が計算される。1つは割線法により以下の式で求められる。 :<math> s = b_k - \frac{b_k-b_{k-1}}{f(b_k)-f(b_{k-1})} f(b_k), </math> 2つめは二分法により求められる。 :<math> m = \frac{a_k+b_k}{2}</math> 割線法を用いた場合の {{math|''s''}} が {{math|''b''<sub>''k''</sub>}} と {{math|''m''}} の間にある場合は次の反復で {{math|''b''<sub>''k''+1</sub> {{=}} ''s''}} となり、そうでない場合は中間点が使用され {{math|''b''<sub>''k''+1</sub> {{=}} ''m''}} となる。 そして、{{math|''f''(''a''<sub>''k''+1</sub>)}} と {{math|''f''(''b''<sub>''k''+1</sub>)}} とが逆符号を持つような、新しい反対点 {{math|''a''<sub>''k''+1</sub>}} が選ばれる。{{math|''f''(''a''<sub>''k''</sub>)}} と {{math|''f''(''b''<sub>''k''+1</sub>)}} が逆符号である場合、反対点は移動せず {{math|''a''<sub>''k''+1</sub> {{=}} ''a''<sub>''k''</sub>}} となる。両者が同符号であれば、{{math|''f''(''b''<sub>''k''+1</sub>)}} と {{math|''f''(''b''<sub>''k''</sub>)}} が逆符号となり、新たな反対点は {{math|''a''<sub>''k''+1</sub> {{=}} ''b''<sub>''k''</sub>}} となる。 最終的に、{{math|{{abs|''f''(''a''<sub>''k''+1</sub>)}} < {{abs|''f''(''b''<sub>''k''+1</sub>)}}}} となった場合は、{{math|''a''<sub>''k''+1</sub>}} は {{math|''b''<sub>''k''+1</sub>}} よりよい解の候補値となり、その結果 {{math|''a''<sub>''k''+1</sub>}} と {{math|''b''<sub>''k''+1</sub>}} の値が交換される。 以上がデッカー法による反復計算の1回分である。 == ブレント法 == デッカー法では、関数 ''f'' が滑らかな形状であれば、効率的に解が求められる。しかし、毎回の収束計算で割線法が選択される場合もあるが、このときは ''b''<sub>''k''</sub> の収束が非常に遅い（特に ''b''<sub>''k''</sub> − ''b''<sub>''k''−1</sub> が任意に小さい）この場合、デッカー法を用いると、二分法よりも非常に多くの収束計算が必要となる。 ブレントは、この問題を解決する小さな修正を提案した。割線法で得られた値が次の収束値として採用される前に、満たすべき判定式を追加した。具体的には、以下の2つの不等式である： * 与えられた許容値 <math>\delta</math> に対し、もし1つ前のステップで二分法を用いた場合は、不等式 :: <math> |\delta| < |b_k - b_{k-1}| </math> : による判定が行われ、満たしていなければ二分法が選択され、得た値を次の収束値とする。もし1つ前のステップが補間である場合は、不等式 :: <math> |\delta| < |b_{k-1} - b_{k-2}| </math> : が使用される。 * また、1つ前のステップが二分法の場合は、同じく :: <math>|s-b_k| < \begin{matrix} \frac12 \end{matrix} |b_k - b_{k-1}|</math> : を満たす必要がある。偽であれば二分法が使用され、得た値が次の収束値となる。1つ前が補間であった場合は、 :: <math>|s-b_k| < \begin{matrix} \frac12 \end{matrix} |b_{k-1} - b_{k-2}|</math> : が使用される。この修正により、k番目の収束計算で二分法が使用される場合、追加で行われる収束計算は多くとも <math>2log_2(|b_{k-1}-b_{k-2}|/\delta)</math> 回に収まる、という点が保証される。なぜなら、上記の条件によって、2回連続して行われる補間計算の回数を半分にでき、その後、多くとも <math>2log_2(|b_{k-1}-b_{k-2}|/\delta)</math> 回の収束計算により、値の変化は <math>\delta</math> よりも小さくでき、二分法に切り替えられるためである。ブレントは、二分法を使用する収束計算回数を N 回とするとき、必要な収束計算回数が多くとも合計 ''N''<sup>2</sup> 回であることを証明した。関数 ''f'' が滑らかな形状であれば、ブレント法は逆2次補間または[[線形補間]]を使用し、 [[収束速度]]は[[超線形]]となる。 さらにブレント法では、もし ''f''(''b''<sub>''k''</sub>)、''f''(''a''<sub>''k''</sub>)、''f''(''b''<sub>''k''−1</sub>) が明確な差があれば、（割線法で使用される）線形補間の代わりに逆2次補間を使用する。これにより収束性がわずかに上がる。結果として、''s''（線形補間または逆2次補間で得られる値）が採用される条件は、以下のように変更される：''s'' が (3''a''<sub>''k''</sub> + ''b''<sub>''k''</sub>) / 4 と ''b''<sub>''k''</sub>の間に存在するかどうか == 使用例 == 関数 ''f''(''x'') = (''x'' + 3)(''x'' − 1)<sup>2</sup> の解（''f''(''x'') = 0 となる x の値）を求める場合を例に挙げる。初期値として [''a''<sub>0</sub>, ''b''<sub>0</sub>] = [−4, 4/3] を考える。''f''(''a''<sub>0</sub>) = −25、''f''(''b''<sub>0</sub>) = 0.48148 である（数値はすべて概数とする）したがって、条件 ''f''(''a''<sub>0</sub>)''f''(''b''<sub>0</sub>) < 0 と |''f''(''b''<sub>0</sub>)| ≤ |''f''(''a''<sub>0</sub>)| を満たしている。 [[Image:Brent method example.png|thumb|関数 ''f''(''x'') = (''x'' + 3)(''x'' − 1)<sup>2</sup> のグラフ]] # 最初の収束計算では、以下の2点 (''b''<sub>−1</sub>, ''f''(''b''<sub>−1</sub>)) = (''a''<sub>0</sub>, ''f''(''a''<sub>0</sub>)) = (−4, −25) と (''b''<sub>0</sub>, ''f''(''b''<sub>0</sub>)) = (1.33333, 0.48148) の間で線形補間を行い、''s'' = 1.23256 を得る。(3''a''<sub>0</sub> + ''b''<sub>0</sub>) / 4 と ''b''<sub>0</sub> の間にあるため、この値は採用される。さらに ''f''(1.23256) = 0.22891 であるため、''a''<sub>1</sub> = ''a''<sub>0</sub>、''b''<sub>1</sub> = ''s'' = 1.23256 となる。 # 2回目の収束計算：(''a''<sub>1</sub>, ''f''(''a''<sub>1</sub>)) = (−4, −25) と (''b''<sub>0</sub>, ''f''(''b''<sub>0</sub>)) = (1.33333, 0.48148) と (''b''<sub>1</sub> を用いて逆2次補間を行い、''f''(''b''<sub>1</sub>)) = (1.23256, 0.22891) を得る。これにより、(3''a''<sub>1</sub> + ''b''<sub>1</sub>) / 4 と ''b''<sub>1</sub> の間にある 1.14205 が求められる。不等式 1.14205 − ''b''<sub>1</sub> ≤ ''b''<sub>0</sub> − ''b''<sub>−1</sub> / 2 を満たすため、この値を採用する。さらに ''f''(1.14205) = 0.083582 であるため、''a''<sub>2</sub> = ''a''<sub>1</sub>、''b''<sub>2</sub> = 1.14205 となる。 # 3回目の収束計算：(''a''<sub>2</sub>, ''f''(''a''<sub>2</sub>)) = (−4, −25) と (''b''<sub>1</sub>, ''f''(''b''<sub>1</sub>)) = (1.23256, 0.22891) と (''b''<sub>2</sub>, ''f''(''b''<sub>2</sub>)) = (1.14205, 0.083582) で逆2次補間を行う。収束値として (3''a''<sub>2</sub> + ''b''<sub>2</sub>) / 4 と ''b''<sub>2</sub> の間にある 1.09032 が得られる。しかし、ここでブレントによる追加条件が加わる。不等式 1.09032 − ''b''<sub>2</sub> ≤ ''b''<sub>1</sub> − ''b''<sub>0</sub> / 2 を満たさないため、この値は採用されない。代わりに、区間 [''a''<sub>2</sub>, ''b''<sub>2</sub>] の中間点 ''m'' = −1.42897 を計算しこれを採用する。''f''(''m'') = 9.26891 であり、''a''<sub>3</sub> = ''a''<sub>2</sub>、''b''<sub>3</sub> = −1.42897 となる。 # 4回目の収束計算：(''a''<sub>3</sub>, ''f''(''a''<sub>3</sub>)) = (−4, −25) と (''b''<sub>2</sub>, ''f''(''b''<sub>2</sub>)) = (1.14205, 0.083582) と (''b''<sub>3</sub>, ''f''(''b''<sub>3</sub>)) = (−1.42897, 9.26891) で逆2次補間を行う。収束値として 1.15448 が得られるが、この値は区間 (3''a''<sub>3</sub> + ''b''<sub>3</sub>) / 4 と ''b''<sub>3</sub>) には存在しない。したがって、中間点 ''m'' = −2.71449 に置き換わる。''f''(''m'') = 3.93934 であり、したがって ''a''<sub>4</sub> = ''a''<sub>3</sub>、''b''<sub>4</sub> = −2.71449 となる。 # 5回目の収束計算：逆2次補間により、区間内にある −3.45500 を得る。しかし、1つ前の計算に二分法を用いたため、不等式 −3.45500 − ''b''<sub>4</sub> ≤ ''b''<sub>4</sub> − ''b''<sub>3</sub> / 2 を満たす必要がある。この不等式は偽であるため、中間点 ''m'' = −3.35724 を採用する。''f''(''m'') = −6.78239 であるので、''m'' は新たな反対点 (''a''<sub>5</sub> = −3.35724) となり、収束値は (''b''<sub>5</sub> = ''b''<sub>4</sub>) のままである。 # 6回目の収束計算：''b''<sub>5</sub> = ''b''<sub>4</sub> であるため逆2次補間を使用できない。したがって、(''a''<sub>5</sub>, ''f''(''a''<sub>5</sub>)) = (−3.35724, −6.78239) と (''b''<sub>5</sub>, ''f''(''b''<sub>5</sub>)) = (−2.71449, 3.93934) の2点で線形補間を行う。''s'' = −2.95064 が得られ、この値はすべての条件を満たしている。しかし1つ前のステップから収束値が変わっていないので、この結果は採用せず、再度、二分法を用いる。''s'' = -3.03587、''f''(''s'') = -0.58418 となる。 # 7回目の収束計算：再度、逆2次補間を行う。''s'' = −2.99436 であり、すべての条件を満たしている。''f''(''s'') = 0.089961 であるので、''a''<sub>7</sub> = ''b''<sub>6</sub>、 ''b''<sub>7</sub> = −3.00219 となる（条件 |''f''(''b''<sub>7</sub>)| ≤ |''f''(''a''<sub>7</sub>)| を満たすよう、''a''<sub>7</sub> と ''b''<sub>7</sub> の値が交換される。） # 8回目の収束計算：''a''<sub>7</sub> = ''b''<sub>6</sub> であるため逆2次補間は使用できない。線形補間により ''s'' = −2.9999、f(s) = 0.0016 を得る。条件は満たしている。 # 以降の収束計算で、''b''<sub>9</sub> = −3 + 6·10<sup>−8</sup> と ''b''<sub>10</sub> = −3 − 3·10<sup>−15</sup> となり、得られる解が ''x'' = −3 に急速に近づく。 ''（訂正：9回目の計算 : f(s) = -1.4E-07、10回目の計算 : f(s) = 6.96E-12）'' == アルゴリズムの詳細 == '''input''' ''a'', ''b'', and a pointer to a subroutine for ''f'' calculate ''f''(''a'') calculate ''f''(''b'') '''if''' ''f''(''a'') ''f''(''b'') >= 0 '''then''' error-exit '''end if''' '''if''' |''f''(''a'')| < |''f''(''b'')| '''then swap''' (''a'',''b'') '''end if''' ''c'' := ''a'' '''set''' mflag '''repeat until''' ''f''(''b or s'') = 0 '''or''' |''b'' − ''a''| is small enough ''(convergence)'' '''if''' ''f''(''a'') ≠ ''f''(''c'') '''and''' ''f''(''b'') ≠ ''f''(''c'') '''then''' <math> s := \frac{af(b)f(c)}{(f(a)-f(b))(f(a)-f(c))} + \frac{bf(a)f(c)}{(f(b)-f(a))(f(b)-f(c))} + \frac{cf(a)f(b)}{(f(c)-f(a))(f(c)-f(b))} </math> ''(inverse quadratic interpolation)'' '''else''' <math> s := b - f(b) \frac{b-a}{f(b)-f(a)} </math> ''(secant rule)'' '''end if''' '''if''' ''(condition 1)'' ''s'' is not between (3''a'' + ''b'')/4 and ''b'' '''or''' ''(condition 2)'' (mflag is set '''and''' |''s''−''b''| ≥ |''b''−''c''| / 2) '''or''' ''(condition 3)'' (mflag is cleared '''and''' |''s''−''b''| ≥ |''c''−''d''| / 2) '''or''' ''(condition 4)'' (mflag is set '''and''' |''b''−''c''| < |<math>\delta</math>|) '''or''' ''(condition 5)'' (mflag is cleared '''and''' |''c''−''d''| < |<math>\delta</math>|) '''then''' <math> s := \frac{a+b}{2} </math> ''(bisection method)'' '''set''' mflag '''else''' '''clear''' mflag '''end if''' calculate ''f''(''s'') ''d'' := ''c'' ''(d is assigned for the first time here; it won't be used above on the first iteration because mflag is set) ''c'' := ''b'' '''if''' ''f''(''a'') ''f''(''s'') < 0 '''then''' ''b'' := ''s'' '''else''' ''a'' := ''s'' '''end if''' '''if''' |''f''(''a'')| < |''f''(''b'')| '''then swap''' (''a'',''b'') '''end if''' '''end repeat''' '''output''' ''b'' '' or s (return the root)'' == アルゴリズムの実装 == ブレントは[[1973年]]に[[ALGOL|ALGOL 60]]による記述でこのアルゴリズムを公開した<ref name="Brent1973" />。[[Netlib]]にはわずかな修正を行った[[Fortran]]版が公開されている<ref>[http://www.netlib.org/go/zeroin.f zeroin.f - netlib]</ref>。[[MATLAB]]関数 <tt>fzero</tt> にもブレント法が組み込まれている<ref>[https://jp.mathworks.com/help/matlab/ref/fzero.html 非線形関数の根 - MATLAB fzero - MathWorks 日本]</ref>。[[PARI/GP]]の<tt>solve</tt>についても同様である。その他の形で実装されたもの（[[C++]]、[[C言語]]、Fortran）は、書籍 Numerical Recipes (ISBN 4874085601) で参照できる。[[Apache Commons Math]] には、[[Java]]実装版が含まれている<ref>[https://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math4/analysis/solvers/BrentSolver.html BrentSolver (Apache Commons Math 4.0 API)]</ref>。[[GNU Scientific Library]] では gsl_root_fsolver_brent にて実装されている<ref>[https://www.gnu.org/software/gsl/doc/html/roots.html#c.gsl_root_fsolver_brent One Dimensional Root-Finding — GSL 2.5 documentation]</ref>。 == 参考文献 == {{Reflist}}<!--added under references heading by script-assisted edit--> * Atkinson (1989). ''An Introduction to Numerical Analysis'' (2nd ed.), Section 2.8. John Wiley and Sons. ISBN 0-471-50023-2. * R.P. Brent (1973). ''Algorithms for Minimization without Derivatives,'' Chapter 4. Prentice-Hall, Englewood Cliffs, NJ. ISBN 0-13-022335-2. * T.J. Dekker (1969). "Finding a zero by means of successive linear interpolation." In B. Dejon and P. Henrici (eds), ''Constructive Aspects of the Fundamental Theorem of Algebra'', Wiley-Interscience, London. ISBN 0471203009 == 外部リンク == * [http://users.hol.gr/~vpapanik/BrentWiki.dpr Wikipedia Brent algorithm - Implementation in Borland Delphi (console application) by vpapanik - 24 Aug 2010] * [http://users.hol.gr/~vpapanik/BrentWikiSolution.xls Wikipedia Brent example - Solution Table by vpapanik - 8 Sep 2009] * [http://www.netlib.org/go/zeroin.f zeroin.f] at [[Netlib]]. * [https://github.com/duncanmcnae/cpp_brent cpp_brent] for C++. *[http://math.fullerton.edu/mathews/n2003/BrentMethodMod.html Module for Brent's Method] by John H. Mathews * [http://newtonexcelbach.wordpress.com/2010/04/13/the-inverse-quadratic-method-3-brents-method/ Open source Excel VBA functions] and examples, including step by step evaluation of the example problem {{デフォルトソート:ふれんとほう}} [[Category:求根アルゴリズム]] [[Category:エポニム]]
このページで使用されているテンプレート:
テンプレート:Ill
(
ソースを閲覧
)
テンプレート:Lang-en-short
(
ソースを閲覧
)
テンプレート:Math
(
ソースを閲覧
)
テンプレート:Reflist
(
ソースを閲覧
)
ブレント法
に戻る。
ナビゲーション メニュー
個人用ツール
ログイン
名前空間
ページ
議論
日本語
表示
閲覧
ソースを閲覧
履歴表示
その他
検索
案内
メインページ
最近の更新
おまかせ表示
MediaWiki についてのヘルプ
特別ページ
ツール
リンク元
関連ページの更新状況
ページ情報