MATLABによる非線形方程式の解法【第２回】

前回の記事では、１変数の非線形方程式を解く方法について理論（計算のメカニズム）とMATLABでの実装方法を解説しました。

しかし、モデルベース開発（MBD）などでモデルを作成する際には、多変数の非線形方程式（連立非線形方程式）を解きたい場面も多くあります。

本記事では、連立非線形法方程式をニュートン・ラフソン法を用いてMATLABで解く方法を解説します。

前回の記事

本記事はこちらの記事の続きとなります。

2023年1月28日 MATLABによる非線形方程式の解法【第１回】

多変数へ拡張したニュートン・ラフソン法

前回の記事にてニュートン・ラフソン法で非線形方程式を解く場合、下式を繰り返し計算すればよいことを解説しました。

$$ x_{k+1} = x_{k} \,-\, \frac{f(x_k)}{f'(x_k)} \tag{5} $$

実は多変数になっても解き方は全く同じです。多変数の場合は連立方程式になりますので、線形代数学で学習する行列を用いると解くことができます。

カヲル

では、最初から解説するね！

まず、以下のようにn個の非線形方程式があったとします。
（もし、右辺にも項がある場合は移項してゼロにします）

$$ \begin{equation} \left. \, \begin{aligned} & f_1(x_1,x_2,\cdots,x_n) = 0 \\ & f_2(x_1,x_2,\cdots,x_n) = 0 \\ & \qquad \vdots \\ & f_n(x_1,x_2,\cdots,x_n) = 0 \\ \end{aligned} \right\} \tag{11} \end{equation} $$

カヲル

n個の変数があって、n個の重複していない方程式があれば原理的には解くことができるけど、非線形の場合は基本的には数式変形では解くことができないね。（解析的に解けない）

心太

ニュートン・ラフソン法の出番ですね！

ここで、以下のように$x_1,x_2,\cdots,x_n$を$\boldsymbol{x}$ベクトルと置き、ゼロベクトルを$\boldsymbol{0}$と置きます。

$$ \boldsymbol{x} = \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix} , \quad \boldsymbol{0} = \begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0 \end{bmatrix} \tag{12} $$

すると、式（１１）は以下のようになります。

$$ \boldsymbol{f}(\boldsymbol{x}) = \boldsymbol{0} \tag{13} $$

式（５）の$f'(x_k)$に相当する行列を$J$と置くと、第$\boldsymbol{k}$ステップの近似解$\boldsymbol{x}_{k+1}$を求める多変数の場合の繰り返し計算の式は以下となります。

$$ \boldsymbol{x}_{k+1} = \boldsymbol{x}_{k} \,-\, J^{-1}\boldsymbol{f}(\boldsymbol{x_k}) \tag{14} $$

カヲル

行列に割り算は存在しないので、割り算に相当する逆行列というものを用いるよ。

心太

行列$J$の逆行列は$J^{-1}$（ジェーインバース）と書きますね！

ここで、式（１４）でいきなり登場した行列$J$はヤコビ行列（ヤコビアン、Jacobian matrix）と呼ばれ、各関数を各変数で順番に偏微分した成分をもつ行列になります。

$$ J = \frac{\partial \boldsymbol{f}(\boldsymbol{x}_k)}{\partial \boldsymbol{x}} = \begin{bmatrix} \frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \frac{\partial f_1}{\partial x_3} & \cdots & \frac{\partial f_1}{\partial x_n} \\ \frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \frac{\partial f_2}{\partial x_3} & \cdots & \frac{\partial f_2}{\partial x_n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \frac{\partial f_n}{\partial x_1} & \frac{\partial f_n}{\partial x_2} & \frac{\partial f_n}{\partial x_3} & \cdots & \frac{\partial f_n}{\partial x_n} \end{bmatrix}_{x=x_k} \tag{15} $$

【関連記事】多変数関数の微分

偏微分（$\frac{\partial f}{\partial x}$）の計算に関しては、こちらの記事を参考にしてみてください。

2022年6月21日多変数関数の微分のイメージを掴もう！

前回の記事と同様に式（１４）を使って$k=1,2,\cdots$と繰り返し計算を行い、求める計算精度まで到達したら計算を止めます。

多変数の場合の収束条件の式は以下となります。

$$ \sum_{i=1}^n \left( f_i(\boldsymbol{x}_{k+1}) \right)^2 \lt \delta^2 \tag{16} $$

ここで$\delta$は相対誤差と呼ばれるもので、求める計算精度に合わせて設定する「小さな正数」となります。

MATLABによる実装

では、実際にMATLABを用いてニュートン・ラフソン法で方程式を解いてみましょう。

前回同様に答え合わせが簡単にできるほうがよいので、以下のものを例題にします。

例題

以下の回路に流れる電流$I_1$、$I_2$、$I_3$と電圧$E_2$を求めよ。

この例題は非線形の問題ではないので、以下のように簡単に解けます。

$R_1$と$R_2$の合成抵抗は$\frac{R_2 \cdot R_3}{ R_2 + R_3}= 9 \, \Omega$なので、抵抗分圧により$E_2=90\,V$となります。

また、$R_2$、$R_3$と$R_1$の合成抵抗は$9 + 1 = 10\,\Omega$となり、$E_1$と$E_3$との間の電位差は$100\,V$ですのでオームの法則より$I_1=10\,A$となります。

キルヒホッフの法則より$I_1$と$(I_2+I_3)$は等しく、抵抗比で分流されるので$I_2=9\,A$、$I_3=1\,A$となります。

これを遠回りとなりますが、今回はニュートン・ラフソン法で解いていきます。先ずは連立方程式を立てます。

$$ \begin{equation} \left. \, \begin{aligned} & R_1 I_1 = E_1 – E_2 \\ & R_2 I_2 = E_2 – E_3 \\ & R_3 I_3 = E_2 – E_3 \\ & I_1 = I_2 + I_3 \\ \end{aligned} \right\} \tag{17} \end{equation} $$

各定数を代入し、$f(x) = 0$の形に変形します。

$$ \begin{equation} \left. \, \begin{aligned} & f_1(I_1,I_2,I_3,E_2) = \frac{R_1 I_1}{E_1 \,-\, E_2} \,-\, 1 \quad (=0) \\ & f_2(I_1,I_2,I_3,E_2) = \frac{R_2 I_2}{E_2} \,-\, 1 \quad (=0) \\ & f_3(I_1,I_2,I_3,E_2) = \frac{R_3 I_3}{E_2} \,-\, 1 \quad (=0) \\ & f_4(I_1,I_2,I_3,E_2) = I_1 \,-\, I_2 \,-\, I_3 \quad (=0) \\ \end{aligned} \right\} \tag{18} \end{equation} $$

ここで、式（１２）に倣い、$\boldsymbol{x}$を以下のように定義します。

$$ \boldsymbol{x} = \begin{bmatrix} E_2 \\ I_1 \\ I_2 \\ I_3 \end{bmatrix} \tag{19} $$

また、式（１３）に倣うと以下のようになります。

$$ \boldsymbol{f}(\boldsymbol{x}) = \begin{bmatrix} f_1(\boldsymbol{x}) \\ f_2(\boldsymbol{x}) \\ f_3(\boldsymbol{x}) \\ f_4(\boldsymbol{x}) \end{bmatrix} = \begin{bmatrix} \frac{R_1 I_1}{E_1 \,-\, E_2} \,-\, 1 \\ \frac{R_2 I_2}{E_2} \,-\, 1 \\ \frac{R_3 I_3}{E_2} \,-\, 1 \\ I_1 \,-\, I_2 \,-\, I_3 \\ \end{bmatrix} =\boldsymbol{0} \tag{20} $$

さらに、式（１５）のヤコビ行列を計算します。

$$ \begin{eqnarray} J = \frac{\partial \boldsymbol{f}(\boldsymbol{x}_k)}{\partial \boldsymbol{x}} &=& \begin{bmatrix} \frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \frac{\partial f_1}{\partial x_3} & \cdots & \frac{\partial f_1}{\partial x_n} \\ \frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \frac{\partial f_2}{\partial x_3} & \cdots & \frac{\partial f_2}{\partial x_n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \frac{\partial f_n}{\partial x_1} & \frac{\partial f_n}{\partial x_2} & \frac{\partial f_n}{\partial x_3} & \cdots & \frac{\partial f_n}{\partial x_n} \end{bmatrix}_{x=x_k} \\ &=& \begin{bmatrix} \frac{R_1 I_1}{(E_1 \,-\, E_2)^2} & \frac{R_1}{E_1 \,-\, E_2} & 0 & 0 \\ -\,\frac{R_2 I_2}{E_2^2} & 0 & \frac{R_2}{E_2} & 0 \\ -\,\frac{R_3 I_3}{E_2^2} & 0 & 0 & \frac{R_3}{E_2} \\ 0 & 1 & -1 & -1 \end{bmatrix}_{x=x_k} \tag{21} \end{eqnarray} $$

カヲル

もし、少し複雑な微分が出てきたときは、数式処理ソフトを使うと簡単だよ。
無料のものとしてはMaxima（マキシマ）があるよ。

以上をもとにMATLABでスプリクトを書くと以下のようになります。

$f_1(\boldsymbol{x})$～$f_4(\boldsymbol{x})$と$J$は２７行目以降で関数として定義し、メインのルーチンから呼び出すようにしています。

また、式（１６）の収束条件の式はMATLABのnorm関数を用いて、記述を簡略化しています。

clear;clc;

% 定数の設定
R1  = 1;    R2  = 10;   R3  = 90;
E1  = 100;
CON = [R1 R2 R3 E1]';

% 初期の近似値を設定
E2  = 5;
I1  = 5;  I2  = 5; I3  = 5;
x   = [E2, I1, I2, I3]';

for idex = 1:1:100          % 最大反復数を制限
    % 元の関数
    F   = [f1(CON,x) f2(CON,x) f3(CON,x) f4(CON,x)]';
    
    if norm(F) <= 10^-7     % 収束判定          
        break;
    end
    
    % 近似値計算
    x   = x - inv( JCB_M(CON,x) ) * F;
end

disp(x);                    % 出力

% ========== 関数定義 ==========
function y = f1(c, x)
    R1  = c(1,1);   E1  = c(4,1);
    E2  = x(1,1);   I1  = x(2,1);
    
    y   = R1 * I1 / (E1 - E2) - 1;
end

function y = f2(c, x)
    R2  = c(2,1);
    E2  = x(1,1);   I2  = x(3,1);
    
    y   = R2 * I2 / E2 - 1;
end

function y = f3(c, x)
    R3  = c(3,1);
    E2  = x(1,1);   I3  = x(4,1);
    
    y   = R3 * I3 / E2 - 1;
end

function y = f4(c, x)
    I1  = x(2,1);   I2  = x(3,1);   I3  = x(4,1);
    
    y   = I1 - I2 - I3;
end

function j = JCB_M(c, x)
    R1  = c(1,1);   R2  = c(2,1);   R3  = c(3,1);   E1  = c(4,1);
    E2  = x(1,1);   I1  = x(2,1);   I2  = x(3,1);   I3  = x(4,1);
    
    %行列の成分計算    
    df1_dE2     = R1 * I1 / (E1 - E2)^2;
    df1_dI1     = R1 / (E1 - E2);
    df1_dI2     = 0;
    df1_dI3     = 0;
    
    df2_dE2     = - R2 * I2 / E2^2;
    df2_dI1     = 0;
    df2_dI2     = R2 / E2;
    df2_dI3     = 0;
    
    df3_dE2     = - R3 * I3 / E2^2;
    df3_dI1     = 0;
    df3_dI2     = 0;
    df3_dI3     = R3 / E2;
    
    df4_dE2     = 0;
    df4_dI1     = 1;
    df4_dI2     = -1;
    df4_dI3     = -1;
    
    % ヤコビ行列定義    
    j   = [df1_dE2  df1_dI1  df1_dI2  df1_dI3;
           df2_dE2  df2_dI1  df2_dI2  df2_dI3;
           df3_dE2  df3_dI1  df3_dI2  df3_dI3;
           df4_dE2  df4_dI1  df4_dI2  df4_dI3];        
end

上記をMATLABで実行すると以下がコマンドウィンドウに表示されます。

90.0000
10.0000
9.0000
1.0000

$E_2 = 90\,V$、$I_1=10\,A$、$I_2=9A$、$I_3=1A$となり、正しく計算できていることが分かります。

プログラミング関連記事