8. 非線形連立方程式
注 – 例題ワークシートをダウンロードすることができます. 1. 例題ワークシートの入手と使用方法 をご覧ください.
多変数の非線形連立方程式の XLPack による実数解の求め方の例を説明します.
例題 1
次の非線形連立方程式を解く.
f1(x1, x2) = 4x12 + x22 - 16 = 0 f2(x1, x2) = x12 + x22 - 9 = 0
プロットして関数の形を見ると次のようになります.
これより, (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(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) を出発点としたときには左上の解に収束しました.
VBA プログラムを使用した解き方 (2)
リバースコミュニケーション版の VBA サブルーチン 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 を呼び出します. リバースコミュニケーション版の詳細については こちら を参照してください.
このプログラムを実行すると上と同じ結果が得られます.
ソルバーを使用した解き方
XLPackソルバーアドインの「非線形方程式」を使って解くこともできます. B9 および C9セルに数式 (=4*B8^2+C8^2-16 および =B8^2+C8^2-9) が入力されています.
ソルバーについては こちら も参照ください.
例題 2
例題の方程式を次のように少し変形してみます.
f1(x1, x2) = 4x12 + (x2 - 2)2 - 16 = 0 f2(x1, x2) = x12 + x22 - 9 = 0
関数の形は次のようになります.
図のように解が2つしかなくなりました. この関数について上と同じ初期値を使って計算してみると, 次の結果が得られました.
(1, 1) と (-1, 1) から出発した場合には近くの解に収束しています. ところが, (1, -1) と (-1, -1) の場合にはむしろ遠くの解に収束しました. また, (0, 0) の場合には収束せず解を求めることができませんでした. 初期値が解から遠すぎる場合にはうまく目的の解を求めることができないことがあります. できるだけよい初期値を与えることが重要です.
