計算物理学２第 9 回：非線形方程式の解法
ver. 2018/6/21
五次以上の高次の多項式や、複雑な関数では代数解が存在しない場合が多くあります。今回は、一般的な非線形方程式 f (x) = 0 の解 x を数値的に求める方法について説明します。

1 二分法
まずは一次元の場合 f (x) = 0 を考えます。一次元の問題の場合はまずは関数のグラフを書いてみて解のおおよその位置を検討することが大事です。f (x) は図 1 のような関数であったとします。求めたいのは x 軸との交点での x の値です。二分法 (bisection method) ではまず解の存在する区間を見つけ、その区間のサイズを狭めて近似解に収束させます。関数 f (x) の連続な領域で、x = x1 , x2 の二点で f (x1 ) と f (x2 ) の符号が異符号であるものが見つかったとします。図 1 では f (x1 ) < 0, f (x2 ) > 0 としています。このとき、f (x) が連続であることから、解は x1 < x < x2 にあります。続いてこの中点 x3 = (x1 + x2 )/2 を考えます。f (x3 ) の符号は正負のどちらかですが、f (x1 ) と f (x3 ) の符号が同じであれば、x3 のほうが x1 よりも解に近いため解は x3 < x < x2 の範囲にあります。プログラムの中で x1 を x3 で置き換えます。もしも f (x2 ) と f (x3 ) の符号が同じであれば、x3 のほうが x2 よりも解に近く、解は x1 < x < x3 の範囲にあります。x2 を x3 と置き換えます。図 1 では f (x3 ) > 0 であるため、プログラムの中で x2 を x3 で置き換えます。この段階で解はx1 < x < x3 の範囲にあることになり、区間の大きさが半分となりました。再び中点 x4 = (x1 + x3 )/2 でのf (x4 ) の符号を考え、これと同符号の点を置き換えます。図 1 では f (x4 ) < 0 なので x1 を x4 で置き換えます。続いて x5 = (x4 + x3 )/2 での f (x5 ) > 0 であるため、x3 を x5 で置き換えます。このように中点を新しい境界の値とすることで解の存在する領域を半分ずつ狭めてゆき、反復法によって解を求めます。二分法は他の方法に比べて収束に時間がかかりますが (収束の精度は k 回行って (x2 − x1 )/2k )、確実に収束するため信頼性が高いです。

2 はさみうち法・割線法 (セカント法)
二分法では x1 ,x2 の中点を用いて解の存在する領域を狭めていきましたが、中点ではなく、[x1 , x2 ]内 の 任 意 の 点 で も 構 い ま せ ん 。は さ み う ち 法 (false position method/regular falsi method) で は(x1 , f (x1 )), (x2 , f (x2 )) を通る直線が x 軸と交わる点を中点の代わりに用います。このような直線は
y=f (x2 ) − f (x1 )(x − x2 ) + f (x2 )x2 − x1(1)と書けますので、x 軸との交点 x3 は

x3 = x2 −x2 − x1f (x2 )f (x2 ) − f (x1 )(2)となります。f (x3 ) の符号をみて、x1 , x2 のどちらかを二分法と同様に置き換えることで解を求めます。はさみうち法は二分法より速く収束することもありますが、遅いときもあります。
割 線 法 (secant method) は 、最 新 の 2 点 を 使 っ て 直 線 を 引 く 方 法 で す 。こ れ ま で に 反 復 に よ っ て

x1 , x2 , · · · xn が求まっているとき、(xn−1 , f (xn−1 )), (xn , f (xn )) を通る直線が x 軸と交わる点を xn+1 と

1

f(x)

f(x3)>0f(x2)>0x4=(x1+x3)/2x1

x3=(x1+x2)/2x2xx5=(x3+x4)/2
f(x4)<0f(x1)<0

図1

二分法

します。

xn+1 = xn −xn − xn−1f (xn )f (xn ) − f (xn−1 )(3)

はさみうち法では f (xn+1 ) の符号を見て xn−1 , xn のどちらを xn+1 で置き換えるかを判定しますが、割線法では必ず xn−1 で置き換えます。割線法は二分法やはさみうち法よりも早く収束しますが、割線法では新しい
点 xn+1 は xn−1 と xn の間に必ずあるわけではありません。そのため、特に引いた直線が x 軸に平行に近い
場合は xn+1 がもともと考えていた空間 [x1 , x2 ] から大きくはずれてしまう可能性があります。そのような場合は反復を続けても収束解が得られないこともあります。例えば [x1 , x2 ] から外れた点が反復で得られた場合は二分法に戻すなどの対策が必要です。

3 Newton 法
Newton 法 (Newton-Raphson 法) は微分を用いて次の反復点を決定します。適当な値 x1 から計算を開始し、反復法によって真の値に近づけていきます。いま xn での値を使って xn+1 での値を求めます。f (x) をx = xn で Taylor 展開すると
f (xn + ∆x) = f (xn ) + f ′ (xn )∆x + O((∆x)2 )(4)となります。∆x の一次までで近似し、f (x) = 0 とすると

∆x = −2f (xn )f ′ (xn )

(5)

f(x)

f(x)
f(x2)>0x3x2x4x1x1x3x2x4x5

xxx5
f(x1)<0

図2

はさみうち法 (左図) と割線法 (右図)となります。つまり、xn+1 として

xn+1 = xn + ∆x = xn −f (xn )f ′ (xn )(6)として x の値を更新します。これは図 3 のように xn での f (x) の接線が x 軸と交わる点を xn+1 とすることに対応します。Newton 法では xn での微係数の値を計算する必要があります。f (x) が多項式などの場合で解析的に f ′ (x) が計算できる場合はそれを用いればよいですが、f (x) が複雑で直接の微分が不可能である場合(例えば行列要素が x の関数になっている行列の固有値や固有ベクトルを用いて f (x) が求められる場合、など)
は数値微分を用いて f ′ (xn ) を評価します。このとき複雑な公式を用いる必要はなく、[f (xn + δ) − f (xn )]/δ 、あるいは [f (xn + δ) − f (xn − δ)]/(2δ) などで近似するので十分です。後者のほうがよい近似ですが、一点あたり f (x) を 3 回計算する必要があります。Newton 法は割線法よりも早く収束するため非常によく使われますが、初期値の選び方を間違えると失敗します。特に微係数が 0 に近くなる場合 (f (x) が x 軸に対して平行に近い場合) は x 軸との交点が離れた場所に飛んでしまいます。階段関数に近いような振る舞いを示す関数の場合は Newton 法では収束させることができず、二分法を用いる必要があります。

3f(x)

x4

x3x2x1x

図 3 Newton 法

4 多次元の Newton 法
多変数の場合 (x = {x1 , x2 , · · · xn }) を考えます。f (x) = {f1 (x), f2 (x), · · · fn (x)} を n 個の関数として

f1 (x1 , x2 , · · · , xn ) =0,(7)f2 (x1 , x2 , · · · , xn ) =0,(8)..
.
fn (x1 , x2 , · · · , xn ) =0,(9)f (x) = 0

(10)

まとめてベクトルの形で書くとを満たす方程式を解く場合を考えます。二分法のような考え方は一次元の場合は解を含む領域を設定できてよいのですが多次元の場合は応用が難しくなります。一方で Newton 法であれば、多変数関数となっても同じアルゴリズムを拡張して用いることが可能です。
(k)(k)(k)

まず k 回目の近似解 x(k) = {x1 , x2 , · · · , xn } が求まったとします。これを用いて k + 1 回目の近似解

x(k+1) を Newton 法によって求めます。x(k) 周りで f (x) を Taylor 展開すると
fi (x(k) + ∆x(k) ) = fi (x(k) ) +n∑∂fi (x(k) )j=14∂xj(k)∆xj+ O((∆x(k) )2 )(11)

となります。ここで Jacobian を∂fi (x(k) )∂xjJij (x(k) ) ≡(12)と書くとベクトルでf (x(k) + ∆x(k) ) = f (x(k) ) + J(x(k) ) · ∆x(k) + O((∆x(k) )2 )(13)と書けます。この Tayloer 展開を 1 次まで行った f (x) がゼロとなる点として x(k+1) を選ぶと∆x(k) = − J −1 (x(k) ) · f (x(k) ),x(k+1)(k)=x+ ∆x(k)(k)=x−J(14)−1(x(k)) · f (x(k))(15)となります。Jacobian の計算は式 (12) が代数的に計算できない場合は数値微分によって計算し、その逆行列は次元が大きい場合は LAPACK のルーチン DGESV 等を用いて計算します。

5 Broyden 法
Newton 法は多次元の非線形方程式でも使えますが、Jacobian とその逆行列の計算の必要があり、代数的に偏微分が行えない複雑な方程式の場合や、数値微分を行わないと行けない場合はこれらの計算コストが高くなってしまうことがあります。Broyden 法はセカント法を多次元に拡張したもので、Jacobian を計算することなく反復を行うことができます。

Newton 法において Jacobian を行列 B で近似するとします。そうすると、式 (14) はB (k) · ∆x(k) = −f (x(k) )(16)となり、与えられた B (k) を用いて ∆x(k) ≡ x(k+1) − x(k) が決まります。(初期値 B (1) は Jacobian を使うか、単位行列を使う場合もあります) 続いて B (k+1) を多次元セカント法で求めるとすると、これはB (k+1) · ∆x(k) = ∆f (k)(17)となります。ここで ∆f (k) = f (x(k+1) ) − f (x(k) ) です。この式は一次元の場合はセカント法に対応し、B (k+1) は (x(k) , f (x(k) )) と (x(k+1) , f (x(k+1) )) を通る直線の傾きを表します。多次元の場合も同様の意味を持ちますが、一次元の場合と違い、これを解いて一意に B (k+1) を求めることはできません。Broyden 法ではB (k+1) = B (k) +(∆f (k) − B (k) · ∆x(k) ) ⊗ ∆x(k)∆x(k) · ∆x(k)(18)とすることで、式 (17) を満たす行列 B (k+1) が前の反復での値 B (k) および x(k) , x(k+1) , f (x(k) ), f (x(k+1) )の値から計算することができます。Broyden 法の計算に必要なのはこの逆行列で(B (k+1) )−1 = (B (k) )−1 +[∆x(k) − (B (k) )−1 · ∆f (k) ] ⊗ ∆x(k) · (B (k) )−1∆x(k) · (B (k) )−1 · ∆f (k)で与えられます。

5(19)6 演習問題x4x4

関数 f (x) = sin( ) − cos( ) について

(34) f (x) = 0 の 0 < x < 10 の範囲内にある解を二分法で求めよ。
(35) f (x) = 0 の解を割線法で求めよ。
(36) f (x) = 0 の解を Newton 法で求めよ。
関数 f (x, y) = x2 + y 2 − 10, g(x, y) = x − y 3 について

(37) f (x, y) = 0, g(x, y) = 0 の解を多次元 Newton 法を用いて求めよ。
(38) f (x, y) = 0, g(x, y) = 0 の解を Broyden 法を用いて求めよ。

6