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

Tech
Share

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

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

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

Ax=b(a11a12⋯a1na21a22⋯a2n⋮⋮⋱⋮an1an2⋯ann)(x1x2⋮xn)=(b1b2⋮bn)(a11a12⋯a1n0a22′⋯a2n′⋮⋮⋱⋮00⋯ann′)(x1x2⋮xn)=(b1b2′⋮bn′)(10⋯001⋯0⋮⋮⋱⋮00⋯1)(x1x2⋮xn)=(b1′′b2′′⋮bn′′)(x1x2⋮xn)=(b1′′b2′′⋮bn′′)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 は左上から右下に動くというイメージを持つといいと思います。 j は i の中で下の行に動き、 k は j の中で右の列に動くイメージです。 繰返しになりますが、この処理で 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 は右下から左上に動き、 j は i の中で上の行に動くイメージです。

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

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 は前進消去後に上三角行列となり、後退消去後に単位行列となっていることがわかるかと思います。