8. 非線形連立方程式

8.1 概要

一般の \(n\) 変数非線形関数 \(f_i()\) の \(n\) 連立方程式の実数のゼロ点を求める問題を考えます.
\[
\begin{gather}
& f_1(x_1, x_2, \dots, x_n) = 0 \\
& f_2(x_1, x_2, \dots, x_n) = 0 \\
& \vdots \\
& f_n(x_1, x_2, \dots, x_n) = 0 \\
\end{gather}
\] \(x_1, x_2, \dots, x_n\) を要素とする列ベクトルを \(\boldsymbol{x}\) と表すことにし, 同様に関数の列ベクトルを \(\boldsymbol{f}\) と表すことにすれば, これらはまとめて次のようにベクトル表現で書くことができます.
\[
\boldsymbol{f(x)} = 0
\]

8.2 非線形連立方程式の解法

8.2.1 ニュートン法

1 変数関数に対するニュートン法を \(n\) 変数関数の (\(n\) 元連立方程式に拡張することができます.

\(\boldsymbol{J(x)}\) はヤコビ行列を表すものとします.
\[
\boldsymbol{J(x)} =
\begin{pmatrix}
{\partial}f_1/{\partial}x_1 & {\partial}f_1/{\partial}x_2 & \cdots & {\partial}f_1/{\partial}x_n \\
{\partial}f_2/{\partial}x_1 & {\partial}f_2/{\partial}x_2 & \cdots & {\partial}f_2/{\partial}x_n \\
\vdots & \vdots & \ddots & \vdots \\
{\partial}f_n/{\partial}x_1 & {\partial}f_n/{\partial}x_2 & \cdots & {\partial}f_n/{\partial}x_n \\
\end{pmatrix}
\] \(\boldsymbol{f(x)}\) をテイラー展開すると次のようになります.
\[
\boldsymbol{f(x_{k+1})} = \boldsymbol{f(x_k + d_k)} = \boldsymbol{f(x_k)} + \boldsymbol{J(x)d_k} + \dots
\] 2 次以上の項を無視するとニュートン法の反復式が得られます.
\[
\begin{align}
& \boldsymbol{d_k} = -\boldsymbol{J(x_k)^{-1}}\boldsymbol{f(x_k)} \\
& \boldsymbol{x_{k+1}} = \boldsymbol{x_k} + \boldsymbol{d_k} \\
\end{align}
\] \(||\boldsymbol{d_k}||\) が十分に小さくなったら終了します.

実際に計算するときにはヤコビ行列の逆行列を求めることはせずに次の連立一次方程式を解きます.
\[
\boldsymbol{J(x_k)}\boldsymbol{d_k} = -\boldsymbol{f(x_k)}
\] 1 変数のときと同様に, 求める解に十分に近い初期値を与えないと収束しないことがあります.

8.2.2 準ニュートン法

ヤコビ行列の計算が難しかったり時間がかかったりすることがあるため, これを他のもので近似する方法です.

代表的なものとしては次のブロイデン法があります.
\[
\boldsymbol{B_{k+1}} = \boldsymbol{B_k} + (\boldsymbol{y} – \boldsymbol{B_ks})\boldsymbol{s^T}/(\boldsymbol{s^Ts})
\] ここで, \(\boldsymbol{y} = \boldsymbol{f(x_{k+1})} – \boldsymbol{f(x_k)}, \boldsymbol{s} = \boldsymbol{x_{k+1}} – \boldsymbol{x_k}\) です.

毎回ヤコビ行列を再計算する代わりに, 漸化式を使って更新した \(\boldsymbol{B_k}\) を使用することにより計算効率の向上を図るものです.

8.2.2 パウエルのハイブリッド法

準ニュートン法と最急降下法を組み合わせた方法です.
\[
(\boldsymbol{J(x_k)^TJ(x_k)} + \gamma_k\boldsymbol{I})\boldsymbol{d_k} = -\boldsymbol{J(x_k)^Tf(x_k)}
\] パラメータ \(\gamma_k = 0\) であれば準ニュートン法になります. また, \(\gamma_k\) の値が大きくなると最急降下法に近づいていきます.

8.3 XLPack を使った非線形連立方程式の解き方

非線形連立方程式を解くためには, パウエルのハイブリッド法による VBA サブルーチン Hybrd1 を使うことができます. これは微分 (ヤコビ行列) を差分近似で計算し, ユーザーによる微分計算を不要にしたサブルーチンです. Hybrd1 は XLPack ソルバーから使うこともできます.

非線形連立方程式では実数解が存在しないことがあり, その場合失敗します. また, 複数の実数解が存在する場合には, どの解が求められるかは初期値に依存します.

例題

次の 2 変数非線形連立方程式を解く.
\[
\begin{align}
& f_1(x_1, x_2) = 4x_1^2 + x_2^2 – 16 = 0 \\
& f_2(x_1, x_2) = x_1^2 + x_2^2 – 9 = 0 \\
\end{align}
\] プロットして関数の形を見ると次のようになります.

これより, (0, 0) を中心に対称に 4 つの実数解があることがわかります.

8.3.1 VBA プログラムを使用した解き方 (1)

Hybrd1 を使ったプログラム例を示します.

Sub F(N As Long, X() As Double, Fvec() As Double, IFlag As Long)
    Fvec(0) = 4 * X(0) ^ 2 + X(1) ^ 2 - 16
    Fvec(1) = X(0) ^ 2 + X(1) ^ 2 - 9
End Sub

Sub Start()
    Const NMax = 10
    Dim N As Long, X(NMax) As Double, Fvec(NMax) As Double, XTol As Double
    Dim Info As Long, I As Long
    '--- Initialization
    N = 2
    XTol = 0.000000000001 '1e-12
    For I = 0 To 4
        '--- Input data
        X(0) = Cells(6 + I, 1): X(1) = Cells(6 + I, 2)
        '--- Compute zeros of system of equations
        Call Hybrd1(AddressOf F, N, X(), Fvec(), XTol, Info)
        '--- Output zeros
        Cells(6 + I, 3) = X(0): Cells(6 + I, 4) = X(1): Cells(6 + I, 5) = Info
    Next
End Sub

目的関数を外部サブルーチンとして用意し, 初期値と要求精度を指定して Hybrd1 を呼び出します. Hybrd1 は, 初期値として与えられた近似解を出発点に反復計算を行い, その近傍のゼロ点を求めます. ここでは, 5 つの異なる初期値を与えたときにどのような解が得られるかを見るプログラムとしています.

まず, 4 つあるゼロ点それぞれを目指して 4 つの初期値 (1, 1), (-1, 1), (1, -1), (-1, -1) を試してみます. さらに, 中心にある (0, 0) を初期値にするとどうなるか見てみます.

このプログラムを実行すると, 次の結果が得られました.

それぞれ, 近くの解に収束していますが, (0, 0) を出発点としたときには左上の解に収束しました.

8.3.2 VBA プログラムを使用した解き方 (2)

リバースコミュニケーション版 (RCI) の Hybrd1_r を使ったプログラム例を示します.

Sub Start()
    Const NMax = 10
    Dim N As Long, X(NMax) As Double, Fvec(NMax) As Double, XTol As Double
    Dim Info As Long, I As Long
    Dim XX(NMax - 1) As Double, YY(NMax - 1) As Double, IRev As Long
    '--- Initialization
    N = 2
    XTol = 0.000000000001 '1e-12
    For I = 0 To 4
        '--- Input data
        X(0) = Cells(6 + I, 1): X(1) = Cells(6 + I, 2)
        '--- Compute zeros of system of equations
        IRev = 0
        Do
            Call Hybrd1_r(N, X(), Fvec(), XTol, Info, XX(), YY(), IRev)
            If IRev = 1 Or IRev = 2 Then
                YY(0) = 4 * XX(0) ^ 2 + XX(1) ^ 2 - 16
                YY(1) = XX(0) ^ 2 + XX(1) ^ 2 - 9
            End If
        Loop While IRev <> 0
        '--- Output zeros
        Cells(6 + I, 3) = X(0): Cells(6 + I, 4) = X(1): Cells(6 + I, 5) = Info
    Next
End Sub

目的関数を外部関数として与えるのではなく, IRev = 1 または 2 のときに XX() の値を使って関数値を計算し YY() に入れて再度 Hybrd1_r を呼び出します. RCI の詳細についてはこちらを参照してください.

このプログラムを実行すると上と同じ結果が得られます.

8.3.3 ソルバーを使用した解き方

XLPackソルバーアドインの「非線形連立方程式」を使って解くこともできます. B9 および C9セルに数式 (=4*B8^2+C8^2-16 および =B8^2+C8^2-9) が入力されています.

ソルバーについてはこちらも参照ください.

8.3.4 VBA プログラムを使用した解き方 (3)

例題 2

例題 1 の方程式を次のように少し変形してみます.
\[
\begin{align}
& f_1(x_1, x_2) = 4x_1^2 + (x_2 – 2)^2 – 16 = 0 \\
& f_2(x_1, x_2) = x_1^2 + x_2^2 – 9 = 0 \\
\end{align}
\] 関数の形は次のようになります.

図のように解が 2 つしかなくなりました.

上の例題プログラムにおいて目的関数 Sub F() を少し変更します.

Sub F(N As Long, X() As Double, Fvec() As Double, IFlag As Long)
    Fvec(0) = 4 * X(0) ^ 2 + (X(1) - 2) ^ 2 - 16
    Fvec(1) = X(0) ^ 2 + X(1) ^ 2 - 9
End Sub

上と同じ初期値を使って計算してみると, 次の結果が得られました.

(1, 1) と (-1, 1) から出発した場合には近くの解に収束しています. ところが, (1, -1) と (-1, -1) の場合にはむしろ遠くの解に収束しました. また, (0, 0) の場合には収束せず解を求めることができませんでした. 初期値が解から遠すぎる場合にはうまく目的の解を求めることができないことがあります.

計算が失敗する場合には, 初期値が解から遠すぎないか, また, 方程式が実数解を持っているかをチェックする必要があります.

公開日: 2025/03/31 (更新日: 2025/05/26)