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

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

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

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

Ax=b(a11a12a1na21a22a2nan1an2ann)(x1x2xn)=(b1b2bn)(a11a12a1n0a22a2n00ann)(x1x2xn)=(b1b2bn)(100010001)(x1x2xn)=(b1b2bn)(x1x2xn)=(b1b2bn)x=b\begin{aligned} \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{aligned}

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

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

58
59
60
61
62
63
64
65
66
67
68
69
70
// 前進消去
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 は左上から右下に動くというイメージを持つといいと思います。 ji の中で下の行に動き、 kj の中で右の列に動くイメージです。
繰返しになりますが、この処理で matrixA は上三角行列になります。

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

76
77
78
79
80
81
82
83
84
85
86
87
// 後退消去
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 は右下から左上に動き、 ji の中で上の行に動くイメージです。

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

Parse()
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
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();
0
/* 前進消去 */
72
73
74
Console.WriteLine("前進消去後");
WriteEquations(matrixA, vectorB);
Console.WriteLine();
0
/* 後退消去 */
89
90
91
92
93
94
95
96
97
98
99
100
101
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}")));

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
初期状態
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 は前進消去後に上三角行列となり、後退消去後に単位行列となっていることがわかるかと思います。