C#で簡単なガウス・ジョルダンの消去法

C#でガウス・ジョルダンの消去法を実行してみます。C言語などで書いた場合とほとんど差はありません。
実際にはガウスの消去法の方が計算の効率が良いので、基本的にはこちらが使われることは少ないと思います。

簡単な…といういうことで、軸を選ぶ処理などは実装していません。

ガウス・ジョルダンの消去法では、解くべき連立一次方程式を%%$%\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}%$%%として、%%$%\boldsymbol{A}%$%%を前進消去と呼ぶプロセスにより上三角行列にし、さらに後退消去と呼ぶプロセスにより単位行列化することで%%$%\boldsymbol{x}%$%%を求めます。式で表すと下のようになります。

##$#\begin{align*}
\boldsymbol{A}\boldsymbol{x} & =\boldsymbol{b}\\
\left(\begin{array}{cccc}
a_{11} & a_{12} & \cdots & a_{1n}\\
a_{21} & a_{22} & \cdots & a_{2n}\\
\vdots & \vdots & \ddots & \vdots\\
a_{n1} & a_{n2} & \cdots & a_{nn}
\end{array}\right)\left(\begin{array}{c}
x_{1}\\
x_{2}\\
\vdots\\
x_{n}
\end{array}\right) & =\left(\begin{array}{c}
b_{1}\\
b_{2}\\
\vdots\\
b_{n}
\end{array}\right)\\
\left(\begin{array}{cccc}
a_{11} & a_{12} & \cdots & a_{1n}\\
0 & a_{22}^{\prime} & \cdots & a_{2n}^{\prime}\\
\vdots & \vdots & \ddots & \vdots\\
0 & 0 & \cdots & a_{nn}^{\prime}
\end{array}\right)\left(\begin{array}{c}
x_{1}\\
x_{2}\\
\vdots\\
x_{n}
\end{array}\right) & =\left(\begin{array}{c}
b_{1}\\
b_{2}^{\prime}\\
\vdots\\
b_{n}^{\prime}
\end{array}\right)\\
\left(\begin{array}{cccc}
1 & 0 & \cdots & 0\\
0 & 1 & \cdots & 0\\
\vdots & \vdots & \ddots & \vdots\\
0 & 0 & \cdots & 1
\end{array}\right)\left(\begin{array}{c}
x_{1}\\
x_{2}\\
\vdots\\
x_{n}
\end{array}\right) & =\left(\begin{array}{c}
b_{1}^{\prime\prime}\\
b_{2}^{\prime\prime}\\
\vdots\\
b_{n}^{\prime\prime}
\end{array}\right)\\
\left(\begin{array}{c}
x_{1}\\
x_{2}\\
\vdots\\
x_{n}
\end{array}\right) & =\left(\begin{array}{c}
b_{1}^{\prime\prime}\\
b_{2}^{\prime\prime}\\
\vdots\\
b_{n}^{\prime\prime}
\end{array}\right)\\
\boldsymbol{x} & =\boldsymbol{b}^{\prime\prime}
\end{align*}#$##

数式で表すと簡単なようにも見えますが、処理自体は3重ループも出てくる大きな処理です。ソースコード等はこちらから

まず、前進消去と呼ばれるプロセスから見ていきましょう。

			// 前進消去
			for (int i = 0; i < Dimension - 1; i++)
			{
				for (int j = i + 1; j < Dimension; j++)
				{
					var s = matrixA[j, i] / matrixA[i, i];
					for (int k = i; k < Dimension; k++)
					{
						matrixA[j, k] -= matrixA[i, k] * s;
					}
					vectorB[j] -= vectorB[i] * s;
				}
			}

ここでは、iは左上から右下に動くというイメージを持つといいと思います。jはiの中で下の行に動き、kはjの中で右の列に動くイメージです。
繰返しになりますが、この処理でmatrixAは上三角行列になります。

次は後退消去と呼ばれるプロセスを見ていきます。

			// 後退消去
			for (int i = Dimension - 1; i >= 0; i--)
			{
				vectorB[i] /= matrixA[i, i];
				matrixA[i, i] = 1; // matrixA[i, i] /= matrixA[i, i]
				for (int j = i - 1; j >= 0; j--)
				{
					var s = matrixA[j, i]; // s = matrixA[j, i] / matrixA[i, i]
					matrixA[j, i] = 0; // matrixA[j, i] -= s
					vectorB[j] -= vectorB[i] * s;
				}
			}

ここでは、iは右下から左上に動き、jはiの中で上の行に動くイメージです。

全体としては下のプログラムになります。(一部省略)

		static void Solve()
		{
			const int Dimension = 4;

			double[,] matrixA = new double[Dimension, Dimension]
			{
				{ 2, 3, 1, 4 },
				{ 4, 1, -3, -2 },
				{ -1, 2, 2, 1 },
				{ 3, -4, 4, 3 },
			};

			double[] vectorB = new double[Dimension]
			{
				10,
				0,
				4,
				6
			};

			double[] solutionX;

			Console.WriteLine("初期状態");
			WriteEquations(matrixA, vectorB);
			Console.WriteLine();

			/* 前進消去 */

			Console.WriteLine("前進消去後");
			WriteEquations(matrixA, vectorB);
			Console.WriteLine();

			/* 後退消去 */

			Console.WriteLine("後退消去後");
			WriteEquations(matrixA, vectorB);
			Console.WriteLine();

			// solutionX = new double[Dimension];
			// for (int i = 0; i < Dimension; i++)
			// {
			// 	solutionX[i] = vectorB[i] / matrixA[i, i];
			// }
			solutionX = vectorB;

			Console.WriteLine("解");
			Console.WriteLine(string.Join("\n", solutionX.Select(x => $"{x,8:F4}")));

計算結果は次の通りです。

初期状態
  2.0000   3.0000   1.0000   4.0000  |  10.0000
  4.0000   1.0000  -3.0000  -2.0000  |   0.0000
 -1.0000   2.0000   2.0000   1.0000  |   4.0000
  3.0000  -4.0000   4.0000   3.0000  |   6.0000

前進消去後
  2.0000   3.0000   1.0000   4.0000  |  10.0000
  0.0000  -5.0000  -5.0000 -10.0000  | -20.0000
  0.0000   0.0000  -1.0000  -4.0000  |  -5.0000
  0.0000   0.0000   0.0000 -30.0000  | -30.0000

後退消去後
  1.0000   0.0000   0.0000   0.0000  |   1.0000
  0.0000   1.0000   0.0000   0.0000  |   1.0000
  0.0000   0.0000   1.0000   0.0000  |   1.0000
  0.0000   0.0000   0.0000   1.0000  |   1.0000

解
  1.0000
  1.0000
  1.0000
  1.0000

matrixAは前進消去後に上三角行列となり、後退消去後に単位行列となっていることがわかるかと思います。

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です