ver. 2020/7/26五次以上の高次の多項式や、複雑な関数では代数解が存在しない場合が多くあります。今回は、一般的な非線形方程式f(z) = 0の解xを数値的に求める方法について説明します。
1 ニ分法
まずはぁが一次元の場合、ƒ(ｪ) ＝ 0を考えます。一次元の問題の場合はまずは関数のグラフを書いてみて解のおおよその位置を検討することが大事です。 f(x)は図1のような関数であったとします。求めたいのはx軸との交点でのxの値です。二分法(bisection method)ではまず解の存在する区間を見つけ、その区間のサイズを狭めて近似解に収束させます。関数f(x)の連続な領域で、x = x1,x2の二点でf(x1)とf(x2)の符号が異符号であるものが見つかったとします。図１の例では● x 1, x 2でf (x 1) < 0, f (x 2) > 0であり、ƒ(х)が連続であることから、解はx 1 ＜ x ＜ x 2にあります。続いてこの中点x 3 ＝ (x 1 + x 2 )/ 2を考えます。 f ( x 3 )の符号は正負のどちらかですが、 ● f(x 1 )とf(x 3 )の符号が同じであれば、x 3のほうがx 1よりも解に近いため解はx 3 ＜ x ＜ x 2の範囲に存在● f(x 2)とf(x 3)の符号が同じであれば、x 3のほうがx 2よりも解に近いため解はx 1 < x < x 3の範囲に存在どちらの場合でももともとの解の存在する区間x  < x < x 2を半分に制限することができました。プログラムの中でx 1またはx 2をx 3で置き換え、再び中点を考えてその符号がどちらと同じかを調べます。図1ではf(x 3) > 0であるため、プログラムの中でxをx 3で置き換えます。統いてf(x 4) < 0であるため、x 1をx 4で置き換え。ƒ(x 5) > 0よりx 3をx 5で置き換えます。このように中点を新しい境界の値とすることで解の存在する領域を半分ずつ狭めてゆき、反復法によって解を求めます。二分法は他の方法に比べて収束に時間がかかりますが(収束の精度はk回行って(x 2 − x 1)/2 k )、確実に収束するため信頼性が高いです。
2　はさみうち法・割線法(セカント法)
二分法ではx 1, x 2の中点を用いて解の存在する領域を狭めていきましたが、中点ではなく、[x 1, x 2]内の任意の点でも構いません。はさみうち法(false position method/regular falsi method)ではは( x 1 , f ( x 1 )), ( x 2 , f ( x 2 ))を通る直線がx軸と交わる点を中点の代わりに用います。このような直線は
と書けますので、x軸との交点x 3は
となります。ƒ(x 3)の符号をみて、x 1,x 2のどちらかを二分法と同様に置き換えることで解を求めます。はさみうち法は二分法より速く収束することもありますが、遅いときもあります。割線法(secant method)は、最新の2点を使って直線を引く方法です。これまでに反復によってx 1, x 2, · · x nが求まっているとき、 (x n − 1, f (x n − 1)), (x n , f (x n ))を通る直線がx軸と交わる点をx +1とします。
はさみうち法ではƒ(x,+1)の符号を見てx,−1, x,のどちらをx,+1で置き換えるかを判定しますが、割線法では必ずx n −1で置き換えます。割線法は二分法やはさみうち法よりも早く収束しますが、割線法では新しい点ːπ+1はːπ –1とːπの間に必ずあるわけではありません。そのため、特に引いた直線がː軸に平行に近い場合はx n +1がもともと考えていた空間[ x 1, x 2]から大きくはずれてしまう可能性があります。そのような場合は反復を続けても収束解が得られないこともあります。例えば[x 1, x 2]から外れた点が反復で得られた場合は二分法に戻すなどの対策が必要です。
3
Newton法法
Newton法(Newton-Raphson法)は微分を用いて次の反復点を決定します。適当な値x,から計算を開始し、反復法によって真の値に近づけていきます。いまx,での値を使ってx,+1での値を求めます。f(x)をx = x nでTaylor展開すると
となります。△ xの一次までで近似し、 f ( x ) = 0とすると
となります。つまり、xn+1として
としてxの値を更新します。これは図3のようにx = x,でのf(x)の接線がx軸と交わる点をx,+1とすることに対応します。 Newton法ではx,での微保数の値を計算する必要があります。 f(x)が多項式などの場合で解析的にƒ′(z)が計算できる場合はそれを用いればよいですが、ƒ(z)が複雑で直接の微分が不可能である場合(例えば行列要素がrの関数になっている行列の固有値や固有ベクトルを用いてf (x)が求められる場合、など)は数値微分を用いてƒ′(ｪ,)を評価します。このとき数値微分でやったような複雑な公式を用いる必要はなく、 [ƒ(x, + δ) − ƒ(x,)]/δ、あるいは[ƒ(x, + δ) − ƒ(x, − δ)]/(2δ)などで近似するので十分です。後者のほうがよい近似ですが、一点あたりf(x)を3回計算する必要があります。Newton法は割線法よりも早く収束するため非常によく使われますが、初期値の選び方を間違えると失敗します。特に微係数が0に近くなる場合(ƒ(ｪ)がｪ軸に対して平行に近い場合)は(5)式からΔτが非常に大きくなるため、ｪ軸との交点が離れた場所に飛んでしまいます。階段関数に近いような振る舞いを示す関数の場合はNewton法では収束させることができず、二分法を用いる必要があります。
4　多次元のNewton法
多変数の場合(x = {x 1,x 2,. · x n})を考えます。 f(x) = {f 1(x),f 2(x),. · f n(x)}をn個の関数として
まとめてベクトルの形で書くと
を満たす方程式を解く場合を考えます。二分法のような考え方は一次元の場合は解を含む領域を設定できてよいのですが多次元の場合は応用が難しくなります。一方でNewton法であれば、多変数関数となっても同じアルゴリズムを拡張して用いることが可能です。まずた回目の近似解x(k) = {z(k),x(k),...,x(k)}が求まったとします。これを用いてk + 1回目の近似解x (k+1)をNewton法によって求めます。 x (k)周りでf (x)をTaylor展開すると
となります。ここでJacobianを
と書くとベクトルで
と書けます。このTayloer展開を１次まで行ったƒ(α)がゼロとなる点としてæ(k+l)を選ぶと
となります。Jacobianの計算は式(12)が代数的に計算できない場合は数値微分によって計算し、その逆行列を求めて式(14)に用います。逆行列の数値計算法についてはこの授業では扱いませんが、興味のある方は前年度の資料を参照してください。
5    Broyden法
Newton法は多次元の非線形方程式でも使えますが、Jacobianとその逆行列の計算の必要があり、代数的に偏微分が行えない複雑な方程式の場合や、数值微分を行わないといけない場合はこれらの計算コストが高くなってしまうことがあります。Broyden法はセカント法を多次元に拡張したもので、Jacobianを計算することなく反復を行うことができます。 Newton法でのJacobian J i ( x (k ))を行列B (k )で近似するとします。そうすると、式(14)は
となり、与えられたB(k)を用いてΔ式(k) ≡式(k+1) −式(k)が決まります。(初期値B(1)はJacobianを使うか、単位行列を使う場合もあります)続いてB(k+l)を多次元セカント法で求めるとすると、これは
となります。ここでΔƒ(k) = ƒ(x(k+1)) − ƒ(x(k))です。この式は一次元の場合はセカント法に対応し、 B(k+l)はは(x(k), f(x(k)))と(x(k+l), f(x(k+l)))を通る直線の傾きを表します。多次元の場合も同様の意味を持ちますが、一次元の場合と違い、これを解いて一意にB(k+1)を求めることはできません。Broyden法では
(ただし(α ⑧ b)ij  = α,bj)とすることで、式(17)を満たす行列B(k+l)が前の反復での値B(k)およびx (k), x (k+1), f ( x (k )), f ( x (k +1))の値から計算することができます。Broyden法の計算に必要なのはこの逆行列で
で与えられます。
演習問題
6
第9回分の演習問題のサンプルプログラムはFortran, C言語それぞれ$ wget https://wwwnucl.ph.tsukuba.ac.jp/~hinohara/compphys2-20/src/fortran9.tgz $ wget https://wwwnucl.ph.tsukuba.ac.jp/~hinohara/compphys2-20/src/c9.tgzとすることでダウンロードできます。解凍は$ tar zxvf fortran9.tgz $ tar zxvf c9.tgz − cos ( − )について関数f ( x ) = sin ( x ) (22)　f(x) = 0の0 < x < 10の範囲内にある解を二分法で求めよ。 (23)   f(x) = 0の解を割線法で求めよ。 (24)   f ( x ) = 0の解をNewton法で求めよ。関数f(x,y) = x² + y² − 10, g(x,y) = x − y³について(25)   f(x,y) = 0, g(x,y) = 0の解を多次元Newton法を用いて求めよ。 (26)   f(x,y) = 0, g(x,y) = 0の解をBroyden法を用いて求めよ。