2 線形方程式と反復法

さて、というわけで、2次元の長方形領域上でのポアソン方程式の差分近似は、 連立線形1次方程式(14)になるので、これを何らかの方法 で解けばいい。

空間1次元の場合には、出てくる行列が3重対角という特別な性質を持ったもの であったので、自由度（要素数）に比例する計算量やメモリー使用量で方程式 を解くことができた。

2次元の場合はどうなるかというと、そんなにうまくはいかないということが わかる。もちろん、2次元の場合でも、 という ように添字をふりなおして(以下、添字が1つしかないときにはこのように付け変えたものを考える ということにする)、以下のように行列の形に書くことはできる。

ここで、 は のベクトル表示である。 しかし、 A の要素を考えてみると、対角要素、その両側の他に、 だけ離れたところに 0 でない要素が現れる。

このために、計算量としては ( として)になり、1次 元のときよりずっと計算量が多い。例えば 、つまり正方形領域と すれば計算量は結局になってしまう。

もちろん、これでも一般の場合のガウスの消去法の計算量である より は少ないが、しかし非常に多いことには変わりはない。また、3次元の場合に はこの状況はいっそう悪化する。

このために、ガウスの消去法のようにまともに行列を解くのではなく、適当な 近似を繰り返していくことで解を求めることはできないかということが昔から 考えられてきた。以下、その方法のいくつかを紹介する。これらを反復法とい う。

2.1 ガウス反復

なんらかの方法で適当な近似解 （ここで上つき添字 i は i 番 目の近似解ということで、別に の i 乗という意味ではない。下につけ ると他の添字と混乱するので上につけてみた）があるとする。これをもうちょっ とましな解にできないかということを考える。

一つの方法は、行列 A を形式的に以下のように分解することである。

ここで、 D は対角要素 だけをとりだした行列で、 N はそれ以 外の残りである。これを使って、以下のような式を考える

D は対角行列なので、これは解けていて、実際のプログラムとしては、

for (i=1; i<nx-1;i++){
    for (j=1; j<ny-1;j++){
        unew[i][j] = 0.25*(u[i][j-1]+u[i][j+1]+u[i-1][j]+u[i+1][j])
                    - rho[i][j]*0.25*deltax*deltax;
    }
}

この反復で、収束すれば、つまり、 となれば、これら が元の方程式を満たしていることは明らかであろう。収束するかどうかは面倒 なのでここでは省くが、上のような簡単なポアソン方程式の離散化で妙な非線 形項とかがなければ収束するということが証明できる。

実用上は、収束したかどうか判定する必要がある。これには、普通は を計算してそれが適当な設定した値よりも小さく なったらいいことにする。

理論的には、数列は連続する2項間の差が小さくなっても収束に近付いている とは限らないが、ガウス反復の場合には、丸め誤差がなければこれで判定でき ることが証明できる。

が、この方法には2つ問題がある。一つは収束が非常に遅いことであり、もう 一つは u と unew で配列が2ついるのでメモリがたくさん必要に なることである。

2.2 ガウス・ザイデル法

さて、実は、上の2つの問題は一度に解決することができる。上のプログラム を、以下のように変えてしまうのである。

for (i=1; i<nx-1;i++){
    for (j=1; j<ny-1;j++){
        u[i][j] = 0.25*(u[i][j-1]+u[i][j+1]+u[i-1][j]+u[i+1][j])
                    - rho[i][j]*0.25*deltax*deltax;
    }
}

形式的には、これは以下のように書くことができる。 行列 A を、対角成分 D、非対角成分の上半三角分 U、残り L の3個に分ける (A=D+L+U)と、

ということになる。

2.3 SOR

というわけで、ガウス・ザイデルはガウス反復よりはましだが、まだもうちょっ となんとかならないものかという気もする。実際に近似解の挙動を見るとわか るが、どちらの方法でも、誤差の減り方にはパターンがあって、1方向から真 の解に近付いている。

というわけで、以下のようなことが考えられる：ガウス・ザイデル法で、修正 量 に適当な係数 を掛けて水 増ししてやればもうちょっとうまくいくのではないか？

こんな安直な方法で大丈夫なのかと思うであろうが、うまくいってしまうのが すばらしいところである。プログラムとしては、ガウス・ザイデルの時の式を ちょっと変えて、

        u[i][j] += omega*(0.25*(u[i][j-1]+u[i][j+1]+u[i-1][j]+u[i+1][j])
                    - rho[i][j]*0.25*deltax*deltax -u[i][j])

というだけで済む。実際的な問題は、 の値をどうとるかということ である。理論的には でないといけないことがわかっていて、 さらに本当はもちろん線形のポアソン方程式とかなら最適な の値 が計算できるが、ここではそこまでは立ち入らないことにする。

この方法を SOR (Successive Over Relaxation, 本によっては Simultaneous ... となっているものもあるかもしれない)という。