nonliniar

非線型連立方程式

Octave では非線型連立方程式の解を得る事ができます。Octave で解ける非線型連立方程式は次の型をしている必要があります。この場合未知数 x はベクトルになります。

F (X) = 0

fsolve の書式は次のようになっています

[X, INFO] = fsolve (FCN, X0)

F(X) 型の関数と初期値ベクトル X0 が与えられたとき fsolve は F(X) == 0 型の方程式を解いて解 X と解の収束情况の情報 INFO を返します。INFO は数値で返されますが、perror ("fsolve", 1) で情報の内容を知ることができます。

fsolve_options (OPT, VAL)関数は、fsolve 関数の設定値を見たり変更したりするのに使います。オプション名だけを引数にすると、その値を表示します。引数なしで実行した場合、全てのオプションの名称と値を表示します。

非線型連立方程式を解く fsolve という関数は MINPACK の hybrd サブルーチンに基づいています。それでは、次の例題を解いてみましょう

     -2x^2 + 3xy   + 4 sin(y) = 6
      3x^2 - 2xy^2 + 3 cos(x) = -4

この二元連立方程式を F(X) = 0 の形にするために、x、y を二次元ベクトル X の要素 x1, x2 と考え F(X) の像 Y の要素が上の式の右辺の定数を左辺に移項した式の値を各々 y1、y2 とするような関数と考え、次のように定義します。

     function y = f (x)
       y(1) = -2*x(1)^2 + 3*x(1)*x(2)   + 4*sin(x(2)) - 6;
       y(2) =  3*x(1)^2 - 2*x(1)*x(2)^2 + 3*cos(x(1)) + 4;
     endfunction

そうすると上の関数 f(x) が[x1; x2] に対し [y1; y2] = [0; 0] を返すときの x1、x2が求める解の値となります。従って、この方程式は次のように解くことができます。

     [x, info] = fsolve ("f", [1; 2])

結果は次のようになります。

上の結果で info = 1 は解が収束したことを表しています。info の数字の表す意味は perror("fsolve", 1) で知ることができます。

上の例で定義した関数 F(X) にはどんな意味があるのでしょうか。それを調べるためにまず X ベクトルを表す x1 - x2 平面上の格子点を次のように作成します。

octave:1> function y = f(x)
> y(1) = -2*x(1)^2 + 3*x(1)*x(2) + 4*sin(x(2)) - 6;
> y(2) = 3*x(1)^2 - 2*x(1)*x(2)^2 + 3*cos(x(1)) + 4;
> endfunction
octave:2> a = linspace (1, 3, 5)';
octave:3> e = ones (5,1);
octave:4> b = [a, e * a(1)];
octave:5> for i = (2:5)
> b = [b; [a, e * a(i)]];
> endfor
octave:6> gplot b with points

次に Y ベクトルを表す y1 - y2 平面上に上の格子点の像をプロットします。

octave:7> c = f(b(1,:))';
octave:8> for i = (2:25)
> c = [c; f(b(i,:))'];
> endfor
octave:9> gset grid
octave:10> gplot c with points

これを見ると写像 F(X) が薄いゴム膜の上に書いた格子点がゴムを変型させたときに移動して格子構造が変型するようなものだと言うことが分かります。F(X) = [0,0] の解は Y 平面上の原点の X 平面上の逆像です。

上の連立方程式の解空間を別の角度からながめてみましょう。いま g(x, y) = -2*x^2 + 3*x*y + 4*sin(y) - 6 という関数と h(x, y) = 3*x^2 - 3*x*y^2 + 3*cos(x) + 4というふたつの関数を考えます。すると上の連立方程式は少なくとも g(x, y) = h(x, y) を満たさないと行けません。この式を満足する点は、z = g(x, y)という曲面と、z = h(x, y)という曲面が交わる曲線上にあります。gnuplot で z = g(x, y) と z = h(x, y) を表示すると次のようになります。

gnuplot> set hidden3d
gnuplot> splot -2*x**2 + 3*x*y + 4*sin(y) - 6, 3*x**2 - 2*x*y**2 + 3*cos(x) + 4

このグラフでふたつの曲面の交わる曲線上に F(X) = [0,0] の解を表す点があります。しかし良く見るとふたつの曲面が交わってできる曲線はふたつあることが分かります。手前の曲線は z = 0 平面の下にあるため F(X) = [0,0] を満たすことはありません。しかし向こう側の曲線は z = 0 平面と交わることができるようです。

したがって、X0 = [10, -10] という初期値を与えてやると手前側の曲線に収束しようとして F(X) = [0, 0] の解が得られないことが予想されます。そこで Octave で計算してみるとやはり解が収束しないことが分かります。

octave:1> function y = f(x)
> y(1) = -2*x(1)^2 + 3*x(1)*x(2) + 4*sin(x(2)) - 6;
> y(2) = 3*x(1)^2 - 2*x(1)*x(2)^2 + 3*cos(x(1)) + 4;
> endfunction
octave:2> [x, info] = fsolve ("f", [10, -10])
x =

   0.12407
  -4.81885

info = 3

octave:3> perror ("fsolve", 3)
iteration is not making good progress

しかし、初期値を [10; 10] にすると向こう側の曲線に収束するので F(X) = [0; 0] の解を導くことができます。

octave:4> [x, info] = fsolve("f", [10,10])
x =

  2.5922
  2.0412

info = 1

octave:6> perror ("fsolve", 1)
solution converged to requested tolerance

この解は最初の解とは違いますが、これは向こう側の曲線が z = 0 平面を２回横切るからです。このように Octave と gnuplot を使うと解析学を視覚的に理解することができます。( gnuplot で次のようにグラフを描かせるとこのあたりの事情がよくわかります。)

gnuplot> set hidden3d
gnuplot> splot [-1:6][-5:5] -2*x**2+3*x*y+4*sin(y)-6, 3*x**2-2*x*y**2+3*cos(x)+ 4, 0