連立一次方程式と行列
簡単な連立一次方程式は手計算で解くことも可能ですが、巨大な連立方程式を手計算で解くのは難題です。例えば授業の後半で習う偏微分方程式を陰解法で数値的に解く問題では、巨大な係数行列を持つ連立一次方程式を解くのと同じ作業を行うことになります。このような問題を手計算で行うのは不可能です。流体に関する偏微分方程式を解く場合にも、巨大な連立一次方程式を解く必要が生じます。巨大な連立一次方程式を高速に解く方法は色々開発されており、数値計算の本で、連立一次方程式(行列)に関することだけを扱った本なども出版されています。
ここではまず、その基礎となる、最も基本的な手法の1つであるガウス・ジョルダンの消去法に関する演習を行います。
ガウス・ジョルダンの消去法
次のような連立一次方程式を考えます。
$$ 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\cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots $$ $$ a_{n1}x_1 + a_{n2}x_2 + \cdots + a_{nn}x_n = b_n $$
これを行列を用いて表現すると、次のように書けます。
$$\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)\ \ \ \ \ \ \ \ \ \ (1) $$
消去法というのは、方程式の両辺に0でない数を掛けたり、割ったり、また方程式を他の方程式に加えたり、引いたりしても解が変わらないという性質を利用して、式(1)の連立方程式の係数を次々に消去して、
$$\left(\begin{array}{cccc} 1 & & & {\Huge 0}\\ & 1 & & \\ & & \ddots & \\ {\Huge 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)$$
の形に変形することです。
さて、ガウス・ジョルダンの消去法を説明します。まず係数行列と右辺の列ベクトルを合わせて、次のような拡大行列$\hat{A}$を作ります。
$$\hat{A}=\left(\begin{array}{cccccc} a_{11} & a_{12} & a_{13} & \cdots & a_{1n} & b_1\\ a_{21} & a_{22} & a_{23} & \cdots & a_{2n} & b_2\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\ a_{n1} & a_{n2} & a_{n3} & \cdots & a_{nn} & b_n \end{array}\right)$$
そして第1行を$a_{11}$で割ります。
$$\left(\begin{array}{ccccc} 1 & \displaystyle{\frac{a_{12}}{a_{11}}} & \displaystyle{\frac{a_{13}}{a_{11}}} & \cdots & \displaystyle{\frac{a_{1n}}{a_{11}}} & \displaystyle{\frac{b_1}{a_{11}}} \\ a_{21} & a_{22} & a_{23} & \cdots & a_{2n} & b_2\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\ a_{n1} & a_{n2} & a_{n3} & \cdots & a_{nn} & b_n \end{array}\right)$$
次に、第1行を除くすべての行に対して、第1行を$a_{i1}(i=2,3,\cdots,n)$倍したものを第$i$行から引きます。このようにして得られた行列を$\hat{A}_1$とします。
$$\hat{A}_1=\left(\begin{array}{ccccc} a_{11}^{(1)} & a_{12}^{(1)} & \cdots & a_{1n}^{(1)} & b_1^{(1)}\\ a_{21}^{(1)} & a_{22}^{(1)} & \cdots & a_{2n}^{(1)} & b_2^{(1)} \\ \vdots & \vdots & \ddots & \vdots & \vdots\\ a_{n1}^{(1)} & a_{n2}^{(1)} & \cdots & a_{nn}^{(1)} & b_n^{(1)} \end{array}\right)$$
ここで
$$a_{1j}^{(1)}= \displaystyle{\frac{a_{1j}}{a_{11}}},\ \ \ b_{1}^{(1)}= \displaystyle{\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}^{(1)}(i=1,3,4,\cdots,n)$倍したものを第$i$行から引きます。その結果、行列$\hat{A}_2$を得ます。
$$\hat{A}_2=\left(\begin{array}{ccccc} 1 & 0 & \cdots & a_{1n}^{(2)} & b_1^{(2)}\\ 0 & 1 & \cdots & a_{2n}^{(2)} & b_2^{(2)} \\ \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & \cdots & a_{nn}^{(2)} & b_n^{(2)} \end{array}\right)$$
ここで
$$a_{2j}^{(2)}= \displaystyle{\frac{a_{2j}^{(1)}}{a_{22}^{(1)}}},\ \ \ b_{2}^{(2)}= \displaystyle{\frac{b_{2}^{(1)}}{a_{22}^{(1)}}},\ \ \ a_{ij}^{(2)}=a_{ij}^{(1)}-a_{i2}^{(1)}a_{2j}^{(2)}, \ \ \ b_{i}^{(2)}=b_{i}^{(1)}-a_{i2}^{(1)}b_{2}^{(2)}$$
ただし$i=1,3,\cdots,n$、$j=2,3,\cdots,n$です。
この操作を順次$n$回繰り返して
$$\hat{A}_n=\left(\begin{array}{ccccc} 1 & 0 & \cdots & 0 & b_1^{(n)}\\ 0 & 1 & \cdots & 0 & b_2^{(n)} \\ \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & \cdots & 1 & b_n^{(n)} \end{array}\right)$$
になったとき、
$$x_i=b_i^{(n)}$$
が求める解となります。
ここで、次の連立一次方程式を解くプログラム例を示します。
$$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)$です。
1 program jordan
2 implicit none
3
4 integer,parameter:: n = 3
5 double precision A(n,n), B(n)
6 integer i,j,k
7
8 data A /2.0d0, 2.0d0, 1.0d0, 1.0d0, 3.0d0, 1.0d0, 1.0d0, 1.0d0, 3.0d0/
9 data B /2.0d0, 4.0d0, -1.0d0/
10
11 call simeq( A, B, n )
12
13 write(*,'(i1,f5.1)') (k,B(k),k=1,n)
14 write(*,'(3f5.1)') ((A(i,j),j=1,n),i=1,n)
15
16 end program jordan
17
18 subroutine simeq( A, B, n )
19 implicit none
20
21 double precision A(n,n), B(n)
22 integer n,i,j,k
23
24 do k = 1, n
25 do j = k + 1, n
26 a(k,j) = a(k,j) / a(k,k)
27 end do
28 b(k) = b(k) / a(k,k)
29
30 do i = 1, n
31 if ( i .ne. k ) then
32 do j = k + 1, n
33 a(i,j) = a(i,j) - a(i,k) * a(k,j)
34 end do
35 b(i) = b(i) - a(i,k) * b(k)
36 end if
37 end do
38
39 end do
40
41 end subroutine
上のプログラムを用いると正しい解が得られますが、最終的な係数行列は単位行列にはなっていません。これは、以下の理由によります。上述の説明のように行列の要素を1や0にしていく際、1や0になった要素はその後の作業に影響を及ぼしません。例えば、0を引くという作業は何も引かないのと同じです。何も引かないのと同じであれば、要素を0にする必要はありません。さらに、要素を0にする必要がないのであれば、対角要素を1にする必要もありません。上記のプログラム例では、この余分な作業をしないようにしてあります。無駄な計算をしないのが効率的なプログラムですが、多少プログラムは読みにくくなってしまいます。
○課題○ 上のプログラムで次の連立一次方程式が解けるかどうか確かめなさい。
$$0x_1 + x_2 + x_3 = 2$$ $$2x_1 + 3x_2 + x_3 = 6$$ $$x_1 + x_2 + 3x_3 = 5$$
○課題○ それでは次の連立一次方程式はどうでしょうか。
$$2x_1 + 3x_2 + x_3 = 6$$ $$0x_1 + x_2 + x_3 = 2$$ $$x_1 + x_2 + 3x_3 = 5$$
最初の課題では、$a_{11}=0$であるため、第1行を$a_{11}$で割る操作のところでエラー(0による除算)が生じるため、正確な解が得られません。一般に計算の途中で対角要素が0になると計算を続けられなくなります。ところが、後者の課題のように、方程式の順番を入れ替えるとうまく計算ができるようになります。このように$a_{kk}= 0$となったときに、方程式の順番を入れ替えることにより、$a_{kk}\neq 0$とする操作のことをピボットの選択といいます。
一般には$k$行列目の方程式を軸に消去を行う場合、$a_{ik}(i=k\sim n)$の中から絶対値が最大のものを選び、それを含む行と$k$行を入れ替えてから消去を行います。これを部分ピボット選択と言います。絶対値の小さい$a_{kk}$は桁落ち計算の結果である可能性もあるので、ピボット選択で絶対値が最大のものを選ぶことにより、$a_{kk}= 0$を防ぐとともに、桁落ちの誤差伝搬を抑える効果もあります。
一方、$a_{ik}(i=k\sim n)$の中からだけでなく、$a_{ij}(i=k\sim n,\ j=k\sim n)$の中から絶対値が最大のものを選び(それを$a_{imax,jmax}$としましょう)、$k$行と$imax$行、$k$列と$jmax$列を入れ替える方法を完全ピボット選択と言います。
部分ピボット選択の方は、単なる方程式の入れ替えですから簡単です。しかし完全ピボット選択の方は、列の入れ替えを含むため、プログラムがかなり複雑になります。というのは、列の入れ替えは変数の順番の入れ替えに相当するため、最後に得られた解が最初の順番と変わってしまうからです。したがって、何列目と何列目を入れ替えたのかを記憶しておいて、最後に解の順番を並べ替えなければなりません。
こちらに簡単な部分ピボット選択のプログラム例を置いておきます。このプログラムを用いれば、上述の最初の課題の連立一次方程式を解くことができます。
さて、このプログラムは係数行列が正則な場合にしか対応していません。係数行列が非正則な場合には、さらに別の解法が必要となります。詳しくは、例えば『Numerical Recipes』の教科書などを参考にしてください。