Gauss-Seidel法と二次元Laplace方程式の数値解法について

これはなに

　おはこんばんにちはてれじょんです。9月に入ると陸の孤島ではクソデカ積乱雲くんや激高気温くんが多く観測され、ようやく夏という風情が出てきました(まじでなに？)。8月もまあまあ暑かったですが、陸の孤島は9月もクソ暑くなりそうなので泣いています。太平洋高気圧、嘘だよな……？
　さて今回はGauss-Seidel法による連立方程式の反復解法について書いて、さらにそれを用いてLaplace方程式を数値的に解いていきたいと思います。私が参考にした本ではLaplace方程式の数値解を求めるときにSOR法を用いていましたが、収束条件の証明をよく理解していない(したがって最適なパラメーターのとり方がよく分かってない)のでGauss-Seidel法を使うことにします。

これはなに
Gauss-Seidel法 is 何
二次元Laplace方程式を解く(理論の説明)
解くぜ！
解析結果
最後に
参考文献

Gauss-Seidel法 is 何

　前回は、名前がよく似たGauss-Jordan法を使って連立方程式を解きました。詳しい説明は前回の記事に譲りますが、これは手計算と同じように行基本変形を繰り返していって解を求める方法でした。一方ここでは行列の変形は行わず、未知変数と同じ数の漸化式をつくり、適当な初期値から始めて逐次値を代入し、誤差項がある値 $\epsilon >0$ より小さくなったら計算を打ち切るという方法で近似解を求めていきます。あくまで求めているのは近似解なので、これを反復法と言ったりもするらしいです。前回扱ったのは直接法です。

　具体的な操作を見ていきましょう。まず一般的に、 $n$ 次実正方行列 $A \in Mat_n(\mathbb{R})$ と $n$ 行縦ベクトル $\boldsymbol{b} \in \mathbb{R}^n$ を用いて $A\boldsymbol{x}=\boldsymbol{b}$ と書ける $n$ 元一次連立方程式の解き方を考えます。適当な初期値 $\boldsymbol{x}^{(0)}=(x_1^{(0)},x_2^{(0)},\cdots,x_n^{(0)})$ を作り、漸化式

${\displaystyle \begin{eqnarray} \left\{ \begin{array}{l} a_{11}x_1^{(\nu +1)}+a_{12}x_2^{(\nu)}+\cdots+a_{1n}x_n^{(\nu)}&=b_1 \\ a_{21}x_1^{(\nu +1)}+a_{22}x_2^{(\nu+1)}+\cdots+a_{2n}x_n^{(\nu)}&=b_2 \\ \vdots \\ a_{n1}x_1^{(\nu +1)}+a_{n2}x_2^{(\nu+1)}+\cdots+a_{nn}x_n^{(\nu+1)}&=b_n \end{array} \right. \end{eqnarray} }$

に上からに代入していくことで順次 $\boldsymbol{x}^{(1)},\boldsymbol{x}^{(2)},\cdots$ を得ます。何度か反復を行うと $\underset{1\leq i \leq n}{\max} |b_i-\sum_{j=1}^n a_{ij}x_{j}^{(\nu)}|<\epsilon$ となるので、十分 $\epsilon$ が小さくなったときに収束したとみなして解を得ます。これがGauss-Seidel法です。

は？？？？？？？？？

おそらく皆さん「これでなんで収束するって言えるの？」と思ったでしょう。私も思いました。色々調べてみたりもしましたが、大抵の場所では収束性についてさして言及していなかったのでちょっと丁寧に書いてみたいと思います。
　当然この方法は一般の連立方程式に対して収束するものではないので、解を得るためにはいくつかの条件が必要です。そこで出てくるのが以下の2つの定理です。

$Th.1$
$H \in Mat_n(\mathbb{R}),~~\boldsymbol{x}^{(\nu)},\boldsymbol{c} \in \mathbb{R}^n$ とする。このとき、スペクトル半径 $\rho(H)(:=\max{|\lambda|},~\lambdaはHの固有値)$ に対して $\rho(H)<1$ ならば反復法 $\boldsymbol{x}^{(\nu+1)}=H\boldsymbol{x}^{(\nu)}+\boldsymbol{c}$ は $\boldsymbol{x}=H\boldsymbol{x}+\boldsymbol{c}$ の一意解に収束する。

【証明】
$\rho(H)<1$ ならば $E-H$ は逆行列をもち、 $(E-H)^{-1}=E+H+H^2+\cdots$ である。
したがって $\boldsymbol{x}=(E-H)^{-1}\boldsymbol{c}$ より、 $\boldsymbol{x}=H\boldsymbol{x}+\boldsymbol{c}$ は一意解を持つ。それを $\boldsymbol{\alpha}$ とする。このとき上記の反復法は
$\boldsymbol{x}^{(\nu)}-\boldsymbol{\alpha}=H(\boldsymbol{x}^{(\nu-1)}-\boldsymbol{\alpha})=H^\nu(\boldsymbol{x}^{(0)}-\boldsymbol{\alpha})\underset{\nu\rightarrow\infty}{\rightarrow}\boldsymbol{0}$
となるので、題意は示された。

$Th.2$
$A \in Mat_n(\mathbb{R})$ が狭義優対角行列または既約優対角行列ならば、Gauss-Seidel法は任意の初期値に対して $A\boldsymbol{x}=\boldsymbol{b}$ の解に収束する。

【証明】
$D,L,U\in Mat_n(\mathbb{R})$ を
$D=\begin{bmatrix} a_{11} &0 & \cdots & 0\\ 0 & a_{22}&\cdots &0 \\ \vdots&\vdots&\ddots &\vdots \\ 0&0&\cdots&a_{nn} \end{bmatrix}~~,L=\begin{bmatrix} 0 &0 & 0&\cdots & 0\\ a_{21} & 0&0&\cdots &0 \\ a_{31}&a_{32}&0&\cdots&0 \\ \vdots&\vdots&\vdots&\ddots &\vdots \\ a_{n1}&a_{n2}&a_{n3}&\cdots&0 \end{bmatrix},\\U=\begin{bmatrix} 0 &a_{12} & a_{13}&\cdots & a_{1n}\\ 0 & 0&a_{23}&\cdots &a_{2n} \\ 0&0&0&\cdots&a_{3n} \\ \vdots&\vdots&\vdots&\ddots &\vdots \\ 0&0&0&\cdots&0 \end{bmatrix}$

で定義する。このとき、Gauss-Seidel法の漸化式より

$\begin{bmatrix} a_{11} &0 & 0&\cdots & 0\\ a_{21} & a_{22}&0&\cdots &0 \\ a_{31}&a_{32}&a_{33}&\cdots&0 \\ \vdots&\vdots&\vdots&\ddots &\vdots \\ a_{n1}&a_{n2}&a_{n3}&\cdots&a_{nn} \end{bmatrix}\begin{bmatrix} x_1^{(\nu+1)} \\ x_2^{(\nu+1)} \\ x_3^{(\nu+1)} \\ \vdots \\ x_n^{(\nu+1)}\end{bmatrix}+\begin{bmatrix} 0 &a_{12} & a_{13}&\cdots & a_{1n}\\ 0 & 0&a_{23}&\cdots &a_{2n} \\ 0&0&0&\cdots&a_{3n} \\ \vdots&\vdots&\vdots&\ddots &\vdots \\ 0&0&0&\cdots&0 \end{bmatrix}\begin{bmatrix} x_1^{(\nu)} \\ x_2^{(\nu)} \\ x_3^{(\nu)} \\ \vdots \\ x_n^{(\nu+1)}\end{bmatrix}=\boldsymbol{b}$
$\therefore (D+L)\boldsymbol{x}^{(\nu+1)}+U\boldsymbol{x}^{(\nu)}=\boldsymbol{b} \\ \boldsymbol{x}^{(\nu+1)}=-(D+L)^{-1}U\boldsymbol{x}^{(\nu)}+(D+L)^{-1}\boldsymbol{b}$

を得る。故に、 $Th.1$ より $\rho(-(D+L)^{-1}U)<1$ を示せばよく、そのために固有値問題

$-(D+L)^{-1}U\boldsymbol{x}=\lambda \boldsymbol{x} ~~\Leftrightarrow~~-U\boldsymbol{x}=\lambda (D+L)\boldsymbol{x}$

を考え、 $\max{|\lambda|}<1$ を示すことにする。

(i)狭義優対角のとき
$|x_m|=\underset{1\leq j \leq n}{\max}|x_j|$ となるような $m$ を取る。
このとき、 $-U\boldsymbol{x}=\lambda (D+L)\boldsymbol{x}$ の第 $m$ 成分を考えて

$-\sum_{j\geq m+1} a_{mj}x_j=\lambda a_{mm}x_m+\lambda \sum_{j\leq m-1}a_{mj}x_j$
$\therefore ~~ -\lambda a_{mm}x_k=\underset{j\geq m+1}{\sum}a_{mj}x_j+\lambda \sum_{j\leq m-1}a_{mj}x_j$

ここで、 $^\exists|\lambda|\geq 1$ を仮定すると、三角不等式より

$\begin{array}{l} 　|\lambda||a_{mm}||x_m| & \leq |\lambda||\sum_{j\neq m}a_{mj}x_j| \\ &\leq |\lambda|\sum_{j \neq m}|a_{mj}||x_j| \\ &\leq |\lambda||x_m|\sum_{j \neq m} |a_{mj}| \end{array}$

$\therefore ~~ |a_{mm}|\leq \sum_{j\neq m}|a_{mj}|$

となって、 $A$ が狭義優対角行列であることに反する。

(ii)既約優対角のとき
狭義優対角のときと同様に $x_m$ をとり、 $J=\left\{j~|~|x_j|=|x_m|\right\}$ とする。
$|\lambda|\geq 1$ として矛盾を導く。

・ $J=\left\{1,2,\cdots,n\right\}$ のとき
このとき、(i)と全く同様にして
$^\forall k,~ |a_{kk}|\leq \sum_{j\neq k}|a_{kj}|$
が言えるので既約優対角性に反する。

・そうでないとき
$^\exists x_j ~s.t.~|x_m|>|x_j|$

このとき

$|\lambda||a_{mm}||x_m|\leq |\lambda|\sum_{j\neq m}|a_{mj}||x_j|\leq|\lambda|(\sum_{j\in J}|a_{mj}||x_j|+\sum_{j\notin J}|a_{mj}||x_j|)~~\cdots(1)$

が成立する。ところで、 $j \notin J$ ならば $a_{mj}=0$ である。実際、もし $a_{mj}\neq 0$ となるような $j \in J$ が存在するとすれば

$\sum_{j\notin J}|a_{mj}||x_j|<\sum_{j\notin J}|a_{mj}||x_m|$

なので、(1)より

$|\lambda||a_{mm}||x_m|<|\lambda||x_m|\sum_{j\neq m}|a_{mj}| \\ \therefore ~~ |a_{mm}|<\sum_{j\neq m}|a_{mj}$

となって、既約優対角性に反する。したがって、

$m\in J ,~ j\notin J \Rightarrow a_{mj}=0$

が分かり、 $A$ が既約であることに反する。以上より $|\lambda|<1$ であり、題意は示された。

証明長すぎん？？？~~tex打ちがめんどい~~

まあ無事 $Th.2$ が示されたので良しとしましょう。これより、狭義優対角または既約優対角行列に対してはGauss-Seidel法が使えることが分かりました。~~でもまだなんか胡散臭いような気がする。~~

二次元Laplace方程式を解く(理論の説明)

　さて、ようやく本題です。ここではDirichlet境界条件下におけるLaplace方程式 $\frac{\partial^2\phi}{\partial x^2}+\frac{\partial^2\phi}{\partial y^2}=0$ を数値的に解いていこうと思います。簡単のため、 $0\leq x,y,\leq 1$ の正方形領域で考えることとし、境界条件は

${\displaystyle \begin{eqnarray} \phi(x,y)=\left\{ \begin{array}{l} 　0&(x=0,1~or~y=1)\\ \sin(2\pi x)&(y=0) \end{array} \right. \end{eqnarray}}$

とします。当然計算機に解かせるには離散化しなきゃいけません。正方形領域を $x$ 方向に $n$ 等分、 $y$ 方向に $m$ 等分して、その分点をそれぞれ $\left\{x_0=0,x_1=1/n,x_2,x_3,\cdots,x_n=1\right\}$ 、 $\left\{y_0=0,y_1=1/m,y_2,y_3,\cdots,y_m=1\right\}$ とおきます。便宜上、 $\phi(x_i,y_j)=\phi_{i,j}$ と書いて、分点の間隔を $\Delta x=1/n,~\Delta y=1/m$ と置きます。このとき、 $\phi$ をTaylor展開することで

$\begin{array}{l} \phi_{i+1,j}&=\phi(x_i+\Delta x,y_j) \\ &=\phi_{i,j}+\frac{\partial\phi}{\partial x}\Delta x+\frac{1}{2}\frac{\partial^2\phi}{\partial x^2}\Delta x^2+\frac{1}{6}\frac{\partial^3\phi}{\partial x^3}\Delta x^3+O(\Delta x^4) \\ \phi_{i-1,j}&=\phi(x_i-\Delta x,y_j) \\ &=\phi_{i,j}-\frac{\partial\phi}{\partial x}\Delta x+\frac{1}{2}\frac{\partial^2\phi}{\partial x^2}\Delta x^2-\frac{1}{6}\frac{\partial^3\phi}{\partial x^3}\Delta x^3+O(\Delta x^4) \\ \end{array}$

を得ます。片々足してよしなにしてやれば

$\left.\frac{\partial^2\phi}{\partial x^2}+\frac{\partial^2\phi}{\partial y^2}\right\|_{x=x_i,y=y_j}=\frac{\phi_{i+1,j}-2\phi_{i,j}+\phi_{i-1,j}}{\Delta x^2}+\frac{\phi_{i,j+1}-2\phi_{i,j}+\phi_{i,j-1}}{\Delta y^2}+O(\Delta x^2)+O(\Delta y^2)$

を得ます。二次以上の項は無視して良さそうなので無視して、

$\frac{\phi_{i+1,j}-2\phi_{i,j}+\phi_{i-1,j}}{\Delta x^2}+\frac{\phi_{i,j+1}-2\phi_{i,j}+\phi_{i,j-1}}{\Delta y^2}=0 \\ \therefore~~\frac{\Delta y^2}{2(\Delta x^2+\Delta y^2)}(\phi_{i+1,j}+\phi_{i-1,j})+\frac{\Delta x^2}{2(\Delta x^2+\Delta y^2)}(\phi_{i,j+1}+\phi_{i,j-1})-\phi_{i,j}=0$

ここで

$\frac{\Delta y^2}{2(\Delta x^2+\Delta y^2)}=\frac{n^2}{2(n^2+m^2)}:=p,~\frac{\Delta x^2}{2(\Delta x^2+\Delta y^2)}=\frac{m^2}{2(n^2+m^2)}:=q$

とおけば、上の式は

$p(\phi_{i+1,j}+\phi_{i-1,j})+q(\phi_{i,j+1}+\phi_{i,j-1})-\phi_{i,j}=0$

となります。だいぶ見通しが良くなりましたね。これと境界条件をあわせれば、全体で $(n-1)(m-1)+2(n+1)+2(m-1)=(n+1)(m+1)$ 本の連立方程式を成します。未知変数は $(n+1)(m+1)$ 個なのであとはこの連立方程式を解けば良いことが分かります。

解くぜ！

　実装しましょう。解くために $\phi_{i,j}$ を一本のベクトルに格納してみます。うまいこと格納してやれば実装も楽になると思うのですが、いい方法が思い浮かばなかったので雑に格納してみます。上で得た $(n+1)(m+1)$ 本の連立方程式を $A\boldsymbol{x}=\boldsymbol{b}$ と書いて

$\boldsymbol{x}=(\phi_{0,0},\phi_{1,0},\cdots,\phi_{n,0},\phi_{0,1},\phi_{1,1},\cdots,\phi_{n,1},\cdots\cdots,\phi_{0,m},\phi_{1,m},\cdots,\phi_{n,m})^T$

と突っ込んでみます。 $\phi_{i,j}$ は $\boldsymbol{x}$ の第 $i+1+(n+1)j$ 成分になっています。ただし、境界条件上の $\phi$ の値は分かっているのでそこにははじめから値を入れておいて、残りのところには初期値として0を入れます。 $i+1+(n+1)j:=l$ と置いてみると、 $A\boldsymbol{x}=\boldsymbol{b}$ の第 $l$ 成分(ただし $\boldsymbol{x}$ の第 $l$ 成分は未知とします)は

$x_l-p(x_{l+1}+x_{l-1})-q(x_{l+n+1}-x_{l-n-1})=0$

となります。こいつをGauss-Seidel法の漸化式とみなせば実装できますね。 $A$ が $Th.2$ の条件を満たすことの証明は省略します。(多分)サクッとできるので。
んで、実装したコードは以下のようになりました。

module globals
   real(8),parameter :: pi=3.141592653589793d0
end module globals

module subprog
   use globals
   implicit none
contains
   function l(i,j,n) result(k)
      integer,intent(IN) :: i,j,n
      integer k
      k=i+1+(n+1)*j
   end function l
end module subprog

program main
   use globals
   use subprog
   implicit none
   integer :: n,m,i,j,ln,fo=11
   real(8),allocatable :: x(:)
   real(8) p,q,dx,dy,ep,tmp,t1,t2

   ep=0.000001d0
   write(*,*)"input step n(x-axis),m(y-axis)"
   read(*,*)n,m
   allocate(x((n+1)*(m+1)))
   x(1:(n+1)*(m+1))=0.0d0
   p=dble(n)**2/2.0d0/(dble(n)**2+dble(m)**2)
   q=dble(m)**2/2.0d0/(dble(n)**2+dble(m)**2)
   dx=1.0d0/dble(n)
   dy=1.0d0/dble(m)
   
   do i=0,n
      x(l(i,0,n))=sin(2.0d0*pi*dx*dble(i))
   enddo
   
   call cpu_time(t1)

1  do i=1,n-1
      do j=1,m-1
         ln=l(i,j,n)
         x(ln)=p*(x(ln+1)+x(ln-1))+q*(x(ln+n+1)+x(ln-n-1))
      enddo
   enddo

   do i=1,n-1
      do j=1,m-1
         ln=l(i,j,n)
         tmp=abs(x(ln)-p*(x(ln+1)+x(ln-1))-q*(x(ln+n+1)+x(ln-n-1)))
         if(tmp>=ep)then
            goto 1
         else
            if(i/=n-1 .or. j/=m-1)then
               CYCLE
            endif
         endif
      enddo
   enddo

   call cpu_time(t2)

   write(*,*)"cal_time=",t2-t1

   open(fo,file="output.d")
   do i=0,n
      do j=0,m
         ln=l(i,j,n)
         write(fo,*)dx*dble(i),dy*dble(j),x(ln)
      enddo
      write(fo,*)" "
   enddo
   close(fo)

end program main