京大理学部 物理科学課題演習 C1 （その４）

#連立一次方程式の解法

　ここでは連立一次方程式の解き方について考えます。この問題は基礎的な数値計算の教科書には必ず現れるものです。世の中には偏微分方程式で定式化される問題がたくさんあります。例えば流体に関する問題は全て偏微分方程式で記述されます。そこで、巨大な連立一次方程式を解く必要性が生じるわけです。記憶領域を節約しながら高速に連立一次方程式を解く方法がいろいろと開発されています。また数値計算法の本の中には、連立一次方程式（行列）に関することだけを扱ったものも出版されています。連立一次方程式を解く方法だけで本が一冊書けるわけですから、ここでその全てを解説するのは不可能です。そこで、最も基本的な手法のひとつである**ガウス・ジョルダンの消去法**について説明することにします。

　次のような連立一次方程式を考えます。

```math
\left.
\begin{array}{l}
  a_{11}x_1 + a_{12}x_2 + \cdots + a_{1n}x_n = b_1 \\
  a_{21}x_1 + a_{22}x_2 + \cdots + a_{2n}x_n = b_2 \\
  \cdots \\
  a_{n1}x_1 + a_{n2}x_2 + \cdots + a_{nn}x_n = b_n \\
\end{array}
\right\}
```

これを行列を用いて表現すると次のように書けます。

```math
\left(
\begin{array}{cccc}
  a_{11} & a_{12} & \ldots & a_{1n} \\
  a_{21} & a_{22} & \ldots & a_{2n} \\
  \vdots & \vdots & \ddots & \vdots \\
  a_{m1} & a_{m2} & \ldots & a_{mn}
\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)
```

消去法というのは、方程式の両辺に $0$ でない数をかけたり、割ったり、また方程式を他の方程式に加えたり、引いたりしても解が変わらないという性質を利用して、上の連立次方程式の係数を次々に消去して、

```math
\left(
  \begin{array}{cccc}
     & 1 &  & {\Large{0}} \\
     &  & 1 \\
     &  &  & \ddots \\
     & {\Large{0}} &  &  & 1
  \end{array}
\right)
\left(
\begin{array}{c}
  x_1 \\
  x_2 \\
  \vdots \\
  x_n
\end{array}
\right) =
\left(
\begin{array}{c}
  c_1 \\
  c_2 \\
  \vdots \\
  c_n
\end{array}
\right)
```

の形に変形することです。さて、ガウス・ジョルダンの消去法について説明します。まず係数行列と右辺の列ベクトルを合わせて、次のような拡大行列 $\tilde{A}$ を作ります。

```math
\tilde{A} =
\left(
\begin{array}{ccccc}
  a_{11} & a_{12} & \ldots & a_{1n} & b_1 \\
  a_{21} & a_{22} & \ldots & a_{2n} & b_2 \\
  \vdots & \vdots & \ddots & \vdots & \vdots \\
  a_{n1} & a_{n2} & \ldots & a_{nn} & b_n 
\end{array}
\right)
```

そして第 1 行を $a_{11}$ で割り、第 1 行を除く全ての行に対して、第 1 行を $a_{i1}\ (i=2, 3, \cdots, n)$ 倍したものを第 i 行から引きます。このようにして得られた行列を $\tilde{A}_1$ とします。

```math
\tilde{A}_1 =
\left(
\begin{array}{ccccc}
  1 & a_{12}^{(1)} & \ldots & a_{1n}^{(1)} & b_1^{(1)} \\
  0 & a_{22}^{(1)} & \ldots & a_{2n}^{(1)} & b_2^{(1)} \\
  \vdots & \vdots & \ddots & \vdots & \vdots \\
  0 & a_{n2}^{(1)} & \ldots & a_{nn}^{(1)} & b_n^{(1)} 
\end{array}
\right)
```

ここで

```math
a_{1j}^{(1)} = \frac{a_{1j}}{a_{11}},\ \ \ \ \ b_{1}^{(1)} = \frac{b_{1}}{a_{11}}\\
a_{ij}^{(1)} = a_{ij} - a_{i1}a_{1j}^{(1)}, \ \ \ \ \ b_{i}^{(1)} = b_i - a_{i1}b_{1}^{(1)}
```

ただし、$i = 2,3,\cdots,n$ 、 $j = 1,2,\cdots,n$ です。次にこの行列の第 2 行を $a_{22}^{(1)}$ で割り、第 2 行を除く全ての行に対して、第 2 行を $a_{i2}\ (i=1, 3, 4, \cdots, n)$ 倍したものを第 i 行から引きます。その結果、行列 $\tilde{A}_2$ を得ます。

```math
\tilde{A}_2 =
\left(
\begin{array}{ccccc}
  1 & 0 & \ldots & a_{1n}^{(2)} & b_1^{(2)} \\
  0 & 1 & \ldots & a_{2n}^{(2)} & b_2^{(2)} \\
  \vdots & \vdots & \ddots & \vdots & \vdots \\
  0 & 0 & \ldots & a_{nn}^{(2)} & b_n^{(2)} 
\end{array}
\right)
```

この操作を順次 $n$ 回繰り返して

```math
\tilde{A}_n =
\left(
\begin{array}{ccccc}
  1 & 0 & \ldots & 0 & b_1^{(n)} \\
  0 & 1 & \ldots & 0 & b_2^{(n)} \\
  \vdots & \vdots & \ddots & \vdots & \vdots \\
  0 & 0 & \ldots & 1 & b_n^{(n)} 
\end{array}
\right)
```

になったとき、

$$x_i = b_i^{(n)}$$

が求める解となります。

　以下に、次の連立一次方程式を解くためのプログラム例を示します。

```math
2x_1 \ +\ \  x_2 \ + \ \ x_3\  =\  2 \\
2x_1\  + 3x_2\  + \ \ x_3 \ = \  4 \\
\ \ \ \ \  x_1\ +\ \ x_2\ + 3x_3 \  = -1
```

答は $(x_1, x_2, x_3) = (1, 1, -1)$ です。

```c
#include <stdio.h>
#define N 3

void simeq(double *apt, double *bpt, int n);

int main(void)
{
    double a[N][N] = {{2.0, 1.0, 1.0},
                      {2.0, 3.0, 1.0},
                      {1.0, 1.0, 3.0}};
    double b[N] = {2.0, 4.0, -1.0};

    int i, j, k;

    simeq(a, b, N);

    for (k = 0; k < N; k++) {
        printf("x%d = %4.1lf\n", k, b[k]);
    }
    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++) {
            printf("%4.1lf", a[i][j]);
        }
        printf("\n");
    }

    return 0;
}

/* 連立一次方程式を解くための関数（ガウス・ジョルダンの消去法） */
void simeq(double *apt, double *bpt, int n)
{
    int i, j, k;

    for (k = 0; k < n; k++) {
        for (j = k + 1; j < n; j++) {
            *(apt + k * n + j) = *(apt + k * n +j) / *(apt + k * n + k);
        }
        *(bpt + k) = *(bpt + k) / *(apt + k * n + k);        
        for (i = 0; i < n; i++) {
            if (i != k) {
                for (j = k + 1; j < n; j++) {
                    *(apt + i * n + j) = *(apt + i * n + j) - *(apt + i * n + k) * *(apt + k * n + j);
                }
                *(bpt + i) = *(bpt + i) - *(apt + i * n + k) * *(bpt + k);
            }
        }
    }

}
```

　上のプログラムを使えば正しく答は求められるのですが、最終的な係数行列が単位行列になっていません。それはなぜでしょうか。解法のところで説明したように、行列の要素を $1$ や $0$ にしていくわけですが、$1$ や $0$ になった要素が、その後他の方程式にどのような影響を与えるかについて考えてみましょう。例えば $0$ を引くということは、何も引かないということと同じです。何も引かないのなら、別に要素を $0$ にする必要はありません。$0$ にする必要がないのなら、対角要素を $1$ にする必要もありません。無駄な計算をしないのが効率的なプログラムです。（ただし多少プログラムは読みにくくなってしまいましたが）

##課題

　上のプログラムで次の連立一次方程式が解けるかどうか確かめなさい。

```math
0x_1 \ +\ \  x_2 \ + \ \ x_3\  =\  2 \\
2x_1\  + 3x_2\  + \ \ x_3 \ = \  6 \\
\ \ \ x_1\ +\ \ x_2\ + 3x_3 \  =\  5
```

　それでは、次の連立一次方程式ではどうでしょうか。

```math
2x_1\  + 3x_2\  + \ \ x_3 \ = \  6 \\
0x_1 \ +\ \  x_2 \ + \ \ x_3\  =\  2 \\
\ \ \ x_1\ +\ \ x_2\ + 3x_3 \  =\  5
```

　前者では $a_{11} = 0$ であるため、エラー（$0$ による除算）が生じてしまい正確な解が得られません。一般に計算の途中で対角要素が $0$ になると計算を続けられなくなります。ところが、後者のように方程式の順番を入れ替えるとうまく計算ができるようになります。このように $a_{kk} = 0$ となったときに、方程式の順番を入れ替えることによって $a_{kk} \neq 0$ となるようにする操作のことを**ピボットの選択**とよびます。これについては[次回](https://qiita.com/sasakitakanori/items/b1a0c54b6a566166d80e)解説します。