8. 非線形連立方程式 (多変数)

多変数の非線形連立方程式のXLPackによる実数解の求め方の例を示します。

例題 1

次の非線形連立方程式を解く。

    f1(x1, x2) = 4x12 + x22 - 16 = 0
    f2(x1, x2) = x12 + x22 - 9 = 0

プロットして関数の形を見ると次のようになります。

08-01

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

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

VBAサブルーチン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(5 + I, 1): X(1) = Cells(5 + I, 2)
        '--- Compute zeros of system of equations
        Call Hybrd1(AddressOf F, N, X(), Fvec(), XTol, Info)
        '--- Output zeros
        Cells(5 + I, 3) = X(0): Cells(5 + I, 4) = X(1): Cells(5 + I, 5) = Info
    Next
End Sub

目的関数を外部サブルーチンとして用意し、初期値と要求精度を指定してHybrd1を呼び出します。Hybrd1は、初期値として与えられた近似解を出発点に反復計算を行い、その近傍のゼロ点を求めます。

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

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

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

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

リバースコミュニケーション版のVBAサブルーチンHybrd1_rを使ったプログラム例を示します。

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
    Dim XX(NMax - 1) As Double, YY(NMax - 1) As Double, IRev As Long, IFlag As Long
    '--- Initialization
    N = 2
    XTol = 0.000000000001 '1e-12
    For I = 0 To 4
        '--- Input data
        X(0) = Cells(5 + I, 1): X(1) = Cells(5 + 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
                IFlag = IRev
                Call F(N, XX(), YY(), IFlag)
            End If
        Loop While IRev <> 0
        '--- Output zeros
        Cells(5 + I, 3) = X(0): Cells(5 + I, 4) = X(1): Cells(5 + I, 5) = Info
    Next
End Sub

目的関数を外部関数として与えるのではなく、IRev = 1または2のときにXX()の値を使って関数値を計算しYY()に入れて再度Hybrd1_rを呼び出します。リバースコミュニケーション版の詳細についてはこちらを参照してください。

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

例題 2

問題の方程式を次のように少し変形してみます。

    f1(x1, x2) = 4x12 + (x2 - 2)2 - 16 = 0
    f2(x1, x2) = x12 + x22 - 9 = 0

関数の形は次のようになります。

08-02

図のように解が2つしかなくなりました。この関数について上と同じ初期値を使って計算してみると、次の結果が得られました。

(1, 1)と(-1, 1)から出発した場合には近くの解に収束しています。ところが、(1, -1)と(-1, -1)の場合にはむしろ遠くの解に収束しました。また、(0, 0)の場合には収束せず解を求めることができませんでした。初期値が解から遠すぎる場合にはうまく目的の解を求めることができないことがわかります。やはり、できるだけよい初期値を与えることが重要です。

公開日: 2015/10/11 (更新日: 2019/10/13) | XLPack 入門編
Top