C# で簡単なガウスの消去法

Tech
Share

C# でガウスの消去法をやってみましたが、ガウス・ジョルダン法と同じく C 言語で書いたものとあまり変わらないものになってしまいました。

何もしていません。あくまでシンプルな処理を紹介しているため、 0 による除算の確認などはしていません。また、単純に高速で解を求めるというようりも、数式に沿った処理をしています。

ガウスの消去法では、解くべき連立一次方程式を Ax=b\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b} として、 A\boldsymbol{A} を前進消去と呼ぶプロセスにより上三角行列にします。ここから、後退代入と呼ぶプロセスにより x\boldsymbol{x} の要素を 1 つずつ解いていきます。式で表すと下のようになります。

Ax=b(a11a12⋯a1na21a22⋯a2n⋮⋮⋱⋮an1an2⋯ann)(x1x2⋮xn)=(b1b2⋮bn)(a11a12⋯a1n0a22′⋯a2n′⋮⋮⋱⋮00⋯ann′)(x1x2⋮xn)=(b1b2′⋮bn′)(x1x2⋮xn)=(b1−∑i=2na1ixia11b2′−∑i=3na2i′xia22′⋮bn′ann′)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}{c} x_{1}\\ x_{2}\\ \vdots\\ x_{n} \end{array}\right) & =\left(\begin{array}{c} \frac{b_{1}-\sum_{i=2}^{n}a_{1i}x_{i}}{a_{11}}\\ \frac{b_{2}^{\prime}-\sum_{i=3}^{n}a_{2i}^{\prime}x_{i}}{a_{22}^{\prime}}\\ \vdots\\ \frac{b_{n}^{\prime}}{a_{nn}^{\prime}} \end{array}\right)\\ \boldsymbol{x} & =\boldsymbol{b}^{\prime\prime} \end{aligned}

ガウス・ジョルダンの消去法よりも式が複雑に見えますが、プログラムにすると、ガウスの消去法の方が簡単で、処理の速さでも実用的です。ソースコード等はこちらから

ガウス・ジョルダンの消去法と全く同じですが、前進消去と呼ばれるプロセスから見ていきます。

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
// 後退代入
for (int i = Dimension - 1; i >= 0; i--)
{
    var s = vectorB[i];
    for (int j = i + 1; j < Dimension; j++)
    {
        s -= matrixA[i, j] * solutionX[j];
    }
    solutionX[i] = s / matrixA[i, i];
}

これは、下から順に solutionX の要素を求めて行く処理です。ある要素は、その下の要素を使うことにより、解が求まります。

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

Solve()
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 = new double[Dimension];

Console.WriteLine("初期状態");
WriteEquations(matrixA, vectorB);
Console.WriteLine();
0
/* 前進消去 */
72
73
74
Console.WriteLine("前進消去後");
WriteEquations(matrixA, vectorB);
Console.WriteLine();
0
/* 後退代入 */
87
88
89
90
91
92
Console.WriteLine("後退代入後");
WriteEquations(matrixA, vectorB);
Console.WriteLine();

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

後退代入後
  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
  1.0000
  1.0000
  1.0000

matrixA は前進消去後に上三角行列となっていますが、ガウス・ジョルダンの消去法と違って後退代入後には変わっていないのがわかります。