連立方程式の解法ー直接法・反復法ー

def solve():
    pivot = 0
    p = 0
    for r in range(0,dof_total):
        pivot = Kc[r][r]
        for c in range(r,dof_total):
            Kc[r][c] = Kc[r][c] / pivot
        F[r] = F[r] / pivot
        for rr in range(r+1, dof_total):
            p = Kc[rr][r]
            for cc in range(r,dof_total):
                Kc[rr][cc] = Kc[rr][cc] - p * Kc[r][cc]
            F[rr] = F[rr] - p * F[r]
    for r in range(dof_total-1,-1,-1):
        U[r] = F[r]
        if r != dof_total-1:
            for c in range(r+1,dof_total):
                U[r] = U[r] - Kc[r][c] * U[c]

ガウス・ジョルダン法

ガウス・ジョルダン法は、ガウスの消去法の変形である。
ガウスの消去法では、\(\mathbf{A}\mathbf{x} = \mathbf{b}\)の連立方程式を\(\mathbf{U}\mathbf{x} = \mathbf{b^*}\)の上三角の係数行列に変形させた。
ガウス・ジョルダン法では、\(\mathbf{A}\mathbf{x} = \mathbf{b}\)の連立方程式を次式のように変形させる。
\begin{equation*}
\left[\begin{array}{cccc}
1 & 0& \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ 0 & 0 & \cdots &0 \\ 0 & 0 & \cdots & 1\\
\end{array}\right]
\left[\begin{array}{c}
x_1 \\ x_2 \\ \cdots\\ x_n
\end{array}\right] =
\left[\begin{array}{c}
b_1^** \\ b_2^** \\ \cdots\\ b_n^**
\end{array}\right]
\tag{4}\end{equation*}
このように変形できれば、解は\(\mathbf{x} = \mathbf{b^{**}}\)となる。

LU分解法

\(\mathbf{A}\mathbf{x} = \mathbf{b}\)の連立1次方程式を解くのに、行列\(\mathbf{A}\)を\(n \times n\)下三角行列\(\mathbf{L}\)と\(n \times n\)上三角行列\(\mathbf{U}\)との積\(\mathbf{LU}\)に分解し、連立1次方程式を次式に変形する。
$$\mathbf{LU}\mathbf{x} = \mathbf{b}\tag{5}$$
\begin{equation*}
\mathbf{L} = \left[\begin{array}{ccccc}
l_{11} & 0 & 0 &  \cdots & 0 \\ l_{21} & l_{22} & 0 & \cdots & 0 \\ l_{31} & l_{32} & l_{33} & \cdots & 0 \\ \cdots & \cdots & \cdots & \cdots & 0 \\ l_{n1} & l_{n2} & l_{n3} & \cdots & l_{nn} \\
\end{array}\right],
\mathbf{U} =
\left[\begin{array}{ccccc}
u_{11} & u_{12} & u_{13} &  \cdots & u_{1n} \\ 0 &u_{22} & u_{23} & \cdots & u_{2n} \\ 0 & 0 & u_{33} & \cdots & u_{3n} \\ 0 & 0 & 0 & \cdots & \cdots \\ 0 & 0 & 0 & \cdots & u_{nn} \\
\end{array}\right]
\tag{6}\end{equation*}
この時、
$$\mathbf{L}\mathbf{y} = \mathbf{b}\tag{7}$$
$$\mathbf{U}\mathbf{x} = \mathbf{y}\tag{8}$$
という関係式が成り立つ。はじめに、式(7)の\(\mathbf{y}\)を求める。次に、\(\mathbf{y}\)を式(8)に代入して、\(\mathbf{x}\)を求める。
つまり、行列\(\mathbf{A}\)を\(n \times n\)下三角行列\(\mathbf{L}\)と\(n \times n\)上三角行列\(\mathbf{U}\)に分解することができれば、あとはガウスの消去法の後退代入と同様の操作で、解を求めることができる。
\(\mathbf{LU}\)分解の方法として、
●\(\mathbf{L}\)の対角成分が\(l_{kk} = 1 ( k = 1, 2, \cdots, n )\)となるように求めるドゥリットル法
●\(\mathbf{U}\)の対角成分が\(u_{kk} = 1 ( k = 1, 2, \cdots, n )\)となるように求めるクラウト法
の2つがある。いずれも係数行列\(\mathbf{A}\)が正則であることが条件である。

反復法

直接法は、マトリクスの成分値を変更することで解を求める方法であったが、反復法では、マトリクスの成分値に変更を加えずに\(\mathbf{x}\)に適当な値を代入していき正解に漸近させていく。初期値ベクトル\(\mathbf{x}^{(0)}\)から出発し、\(\mathbf{x}^{(1)}\)，\(\mathbf{x}^{(2)}\)，\cdots，\(\mathbf{x}^{(k)}\)と順に代入していく。計算を打ち切るための収束判定は、次式である。
$$\left|x_i^{(k+1)} – x_i^{(k)}\right| \lt \epsilon (i = 1, 2, \cdots, n)\tag{9}$$
式(9)を満たす時、\(\mathbf{x}^{(k)}\)は十分に解に近づいたことを意味する。また、\(\mathbf{x}^{(i)}\)を決めるために、次式のような反復公式がある。
$$\left|x_i^{(k+1)} = x_i^k\right| + e_i^{(k)} (i = 1, 2, \cdots, n)\tag{10}$$
この\(e_i^{(k)}\)を算出する方法として、いくつかの方法が存在する。今回は、基本的な3つの方法を紹介する。

ヤコビ法

ヤコビ法の\(e_i^{(k)}\)を算出は次式となる。
\begin{align*}e_i^{(k)}
&=(b_i – a_{i1}x_1^{(k)} – a_{i2}x_2^{(k)} – \cdots – a_{in}x_n^{(k)}/a_{ii}\\
&= \displaystyle\frac{1}{a_{ii}} \left[b_i – \sum_{j=1}^{n}a_{ij}x_j^{(k)}\right]\tag{11}\end{align*}
ヤコビ法は、対角要素が非対角要素と比較して十分に大きいときに成り立つ。

式だけではイメージがしにくいと思うので、例を挙げてみる。
簡単な連立方程式として、次式を考える。
\begin{equation*}\left[\begin{array}{cc}
5& -1 \\1 & 5\\
\end{array}\right]
\left[\begin{array}{c}
x \\ y
\end{array}\right] =
\left[\begin{array}{c}
4 \\ 6
\end{array}\right]\tag{12}\end{equation*}
式(12)を式(11)に代入すると、
$$e_x^{(k)} = (4 – 5x^{(k)} + y^{(k)} )/ 5 ,   x^{(k+1)} = x^{(k)} + e_x^{(k)}$$
$$e_y^{(k)} = (6 – x^{(k)} – 5y^{(k)} )/ 5 ,   y^{(k+1)} = y^{(k)} + e_y^{(k)}$$
となる。これを繰り返し計算することで、解に近づいていく。

ガウス・ザイデル法

ヤコビ法の反復式において、\(x_i^{(k+1)}\)が求められたら、それを以後の\(x_j^{(k+1) (j = i+1, i+2, \cdots, n)}\)の計算に使うことで、ヤコビ法よりも速く収束すると考えられる。
この方法をガウス・ザイデル法をいう。\(e_i^{(k)}\)の算出は次式となる。
$$e_i^{(k)} = \displaystyle\frac{1}{a_{ii}} \left[b_i – \sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)} – \sum_{j=1}^{n}a_{ij}x_j^{(k)}\right]\tag{13}$$
ガウス・ザイデル法はヤコビ法と同様で、対角要素が非対角要素と比較して十分に大きいときに成り立つ。

式(12)を例に用いて、解いてみると次式になる。
$$e_x^{(k)} = (4 – 5x^{(k)} + y^{(k)} )/ 5 ,   x^{(k+1)} = x^{(k)} + e_x^{(k)}$$
$$e_y^{(k)} = (6 – x^{(k+1)} – 5y^{(k)} )/ 5 ,   y^{(k+1)} = y^{(k)} + e_y^{(k)}$$
この式をみると、ヤコビ法と違い\(e_y^{(k)}\)の式に\(x_i^{(k+1)}\)が用いられていることが分かる。

SOR法

ガウス・ザイデル法で得られた\(x_i^{(k+1)}\)を\(\underline{x_i}^{(k+1)}\)と置き換え、これと\(x_i^{(k)}\)との差に係数を乗じたものを修正量とする方法をSOR法または逐次緩和加速法という。
$$x_j^{(k+1)} = x_j^{(k)} + w(\underline{x_i}^{(k+1)} – x_j^{(k)}\tag{14}$$
ここで、\(w\)は緩和係数とよび、\(w \gt 1\)の時にはガウス・ザイデル法で得られる修正量よりも多く修正を加えることになり、一般に収束が速くなる。\(w = 1\)の時はガウス・ザイデル法と一致する。\(w\)を大きくすれば良いというものではなく、最適な値が存在している。