第 3 章 連立 1 次方程式の解法|ベクト ルコン
ピュータに適した反復法
前章に述べた SOR 法と ADI 法は典型的な反復法であるが，その点過緩和または線過緩和の計算は並列処理することができないので，科学技術計算に広く用いられているベクトルコンピュータには向かない．一方

Jacobi 過緩和法は並列処理可能であるが収束性が良くない．本章には並列処理可能でかつ収束性も良いい
くつかの反復法について述べる．本章では読者が前章の内容を理解していることを前提に解説する．

3.1

Chebyshev SOR 法

Chebyshev SOR 法 (red-black point relaxation) は SOR 法をベクトルコンピュータ向きに改良したもの
で，２次元問題の 5 点差分方程式，３次元問題の 7 点差分方程式の境界値問題に適用することができる．ここでは２次元 Poisson 方程式の 5 点差分方程式
1

cij ui;1 + c2 u +1 + c3 u ;1 + c4 u +1 + c0 u = f
j

ij

i

j

ij

ij

ij

ij

ij

ij

(3.1)

ij

の Dirichlet 問題から説明する．計算領域は図 3.1 に示すように長方形で，緩和計算は 1の内点に対して行われる．内点を i+j が偶数になる点 (red point，図に

i

1 1

i 

j

1

j

 印で示す) と奇数になる点 (black

point， 印で示す) に分ける．このとき偶数点の緩和はほかの偶数点のデータと無関係に行えるので，偶数点のみを同時に陽的1に緩和することができる．奇数点についても同様のことが言える．Chebyshev SOR
法はすべての偶数点の点過緩和と，すべての奇数点の点過緩和を交互に解が収束するまで反復するもので
ある．

j
j1
f


4 

3

2 

1
j=0
i=0 1 2 3

i1 i

f



図 3.1: Chebyshev SOR 法の偶数点 と奇数点

1 差分方程式を個々に計算するときに陽的 (explicit) に解くといい，その解法を陽解法と言う．また連立方程式として計算するときに陰的 (implicit) に解くといいその解法を陰解法と言う．

1

2

反

復

法（ Chebyshev SOR 法）

以下には，式 (3.1) の Dirichlet 問題を Chebyshev SOR 法で解くプログラムを示す．
!*******************************************************************************! Solution of Dirichlet problem of elliptic pde by red-black point relaxation!*******************************************************************************SUBROUTINE REDBLK1(u,c,f,r,i1,if,j1,jf,me,nef,mo,nof,al,na,naf,rmean,e)
DIMENSION u(0:if,0:jf),c(i1,j1,0:4),f(i1,j1),r(i1,j1),me(2,nef),mo(2,nof)IF(al<=1.)al=1.5 naf=MAX0(naf,100)
e=AMAX1(e,1.E-5)
ne=0 no=0
DO i=1,i1 DO j=1,j1
IF(MOD(i+j,2)==0)THEN!for red points
ne=ne+1 me(1,ne)=i me(2,ne)=j!input ij in me
ELSE!for black points
no=no+1 mo(1,no)=i mo(2,no)=j!input ij in mo
ENDIF
ENDDO ENDDO
na=0 10 na=na+1
DO n=1,nef i=me(1,n) j=me(2,n)!red pointsr(i,j) = c(i,j,1)*u(i-1,j)+c(i,j,2)*u(i+1,j) &!compute residuals+c(i,j,3)*u(i,j-1)+c(i,j,4)*u(i,j+1) &+c(i,j,0)*u(i,j)-f(i,j)
u(i,j) = u(i,j)-al*r(i,j)/c(i,j,0)!point over-relaxation
ENDDO
DO n=1,nof i=mo(1,n) j=mo(2,n)!black pointsr(i,j) = c(i,j,1)*u(i-1,j)+c(i,j,2)*u(i+1,j) &!compute residuals+c(i,j,3)*u(i,j-1)+c(i,j,4)*u(i,j+1) &+c(i,j,0)*u(i,j)-f(i,j)
u(i,j) = u(i,j)-al*r(i,j)/c(i,j,0)!point over-relaxation
ENDDO
rmean=0.
DO i=1,i1 DO j=1,j1
rmean = rmean+ABS(r(i,j))!compute mean residual
ENDDO ENDDO
rmean=rmean/FLOAT(i1*j1)
IF(rmean>e.AND.na<naf)
GOTO 10!decide convergence
END

このサブルーチンは，配列 u に uij の予測値と境界値，配列 c,f に式 (3.1) の係数 ckij と右辺 fij の値，引数 al, naf, e に緩和係数 ，最大反復数 nf ，収束判定用の小さい正数 " の値，i1,if,j1,jf,nr,nb にも

1 = jf ; 1，また nr と nb は偶数点の数 nred または奇数点の数
nblack で，nred は nblack と同じか１つ多くなる．nred + nblack = i1  j1 ．このプログラムでは，まずすべて値を入れて引用する．ただし i1 = if

;1

 j

の偶数点またはすべての奇数点を簡単に呼出せるように，配列 mr, mb を用意し ，これらの配列にあらかじ
め偶数点または奇数点の i j を入れておく．反復計算は，偶数点の点過緩和，奇数点の点過緩和の順に行わ
れ，平均残差値 rmean

 (P j j) ( 1 1 )
rij

= i j

"

または反復回数 n = nf になるまで繰り返される．計算終了

時に，配列 u に関数値 uij ，配列 r に残差値 rij ，引数 na, rmean に n と rmean の値が入る．このプログラ
ムは，すべての偶数点またはすべての奇数点の緩和計算が並列に行えるように作られているので，ベクトル化率だけでなく並列度も十分高くなるっている．なお念のため書けば，この解法は，パソコンなどのスカラー演算のみを行うコンピュータでは SOR 法よりもかなり劣るもので，ベクトルコンピュータを使用する
大規模計算においてはじめて威力を発揮するものである．
問1

上記問題で，図 3.2 に示すように，右側境界 (i = if ) を Neumann 条件，上下境界 (j = 0 jf ) を周期条件に変更した混合境界値問題の Chebyshev SOR 法のプログラムを示せ．

連立 1 次方程式の解法

Neumann 条件 ()@ u=@ x if j

1

ci

fj

uif

=

gj

を考慮した差分方程式;1 + c3f u f ;1 + c4f u f +1 + c0f u f = f f
j

i

i j

j

i

3

j

i j

i

j

i

j

i

j

と周期条件 ui jf = ui0 +  (@ u=@ y )i jf = (@ u=@ y )i0 を考慮した式

1

ci j uif

ui jf;1 f + c2 f u +1 f + c3 f u f ;1 + c4 f u 1 + c0 f u f = f f 

=

j

ui

ij

i

j

ij

ij

ij

i

ij

ij

ij

0 +

を用いることにする．またここでは，jf を偶数に取り，点過緩和の計算を下方の境界点を除外して行うことにする．このとき奇数点と偶数点の数は等しくなる．次のプログラムは上記のサブルーチン REDBLK1 に倣って書かれている．
!***************************************************************************************! Solution of mixed bound value problem of elliptic pde by red-black point relaxation!***************************************************************************************SUBROUTINE REDBLK2(u,c,f,r,if,jf,phi,me,nef,mo,nof,al,na,naf,rmean,e)
DIMENSION u(0:if,0:jf),c(if,jf,0:4),f(if,jf),r(if,jf),me(2,nef),mo(2,nof)IF(al<=1.)al=1.5 naf=MAX0(naf,1000)
e=AMAX1(e,1.E-5)
ne=0 no=0
DO i=1,if DO j=1,jf
IF(MOD(i+j,2)==0)THEN!for red points
ne=ne+1 me(1,ne)=i me(2,ne)=j!input ij in me
ELSE!for black points
no=no+1 mo(1,no)=i mo(2,no)=j!input ij in mo
ENDIF
ENDDO ENDDO
na=0 10 na=na+1
DO n=1,nef i=me(1,n) j=me(2,n)!red points
ip1=MIN0(i+1,if)!Neumann cond
jp1=j+1 IF(j==jf)jp1=1!periodic condr(i,j) = c(i,j,1)*u(i-1,j)+c(i,j,2)*u(ip1,j) &!compute residuals+c(i,j,3)*u(i,j-1)+c(i,j,4)*u(i,jp1) &+c(i,j,0)*u(i,j)-f(i,j)
u(i,j) = u(i,j)-al*r(i,j)/c(i,j,0)!point over-relaxation
ENDDO
DO n=1,nof i=mo(1,n) j=mo(2,n)!black points
ip1=MIN0(i+1,if)!Neumann cond
jp1=j+1 IF(j==jf)jp1=1!periodic cond
!compute residualsr(i,j) = c(i,j,1)*u(i-1,j)+c(i,j,2)*u(ip1,j) &+c(i,j,3)*u(i,j-1)+c(i,j,4)*u(i,jp1) &+c(i,j,0)*u(i,j)-f(i,j)
u(i,j) = u(i,j)-al*r(i,j)/c(i,j,0)!point over-relaxation
ENDDO
FORALL(i=1:if)u(i,0)=u(i,jf)-phi
rmean=0.
DO i=1,if DO j=1,jf
rmean = rmean+ABS(r(i,j))!compute mean residual
ENDDO ENDDO
rmean = rmean/FLOAT(if*jf)
IF(rmean>e.AND.na<naf)
GOTO 10!decide convergence
END

以下には，ポアソン方程式 uxx + uyy =u;2 の図 3.2 の境界条件の課された２つの境界値問題を，出発値= 0 を与えて解いた結果を示す．計算には格子間隔 0.05 の正方形格子が用いられた．解uij

は，" = 10;5

4

反法（ Chebyshev SOR 法）

復で小数点以下 3 桁までほぼ正確に求めることができる．表 3.1 は緩和係数  と最終反復数 n，その時点における平均残差値 rmean の関係を示したものである．解は  = 1:75 程度の大き目の緩和係数に対して良く収束する．また計算領域の一辺のみに関数値の与えられる混合境界値問題では，収束までにかなりの反復を要することが分かる．なおこの混合境界値問題の厳密解は u(x y ) =Dirichlet
u=1

1

x

1
HH

Dirichlet
u=y

Dirichlet
u=y

y =x0= 0; 2 +2 + である．

Dirichletu=0(a) Dirichlet 問題y

periodicHu(x1)= u(x0)+1 H
Neumann
Dirichlet
u
(
x
1)=
u
(
x
0)
u=y
u =0
y

y

A
A
A
A

y =x0= 0

1x(b) 混合境界値問題

y1

図 3.2: Chebyshev SOR 法の例題
表 3.1: Chebyshev SOR 法の解の収束状況

Dirichlet 問題


n

1.00
1.25
1.50
1.75
1.80

339
209
118
40
50

mean
10;5
10;5
10;5
10;5
10;5

r

混合境界値問題

mean

n

r

1000 63  10;51000 11  10;5
850
10;5
385
10;5
300
10;5

Chebyshev SOR 法は，３次元の 7 点差分方程式にそのまま拡張することができる．すなわち + + が偶数の点と奇数の点に分ければ ，各々の点過緩和の計算を並列処理することができる．図 3.3(a) に示す ２
次元の 9 点高次差分方程式の場合には，mod( + 3) が 0 の点，1 の点，2 の点に分ければ，点過緩和の計
算を並列処理することができる．また (b) 図の 9 点高次差分方程式の場合には，表に示すように 4 色の点に
i

i

r
r
r r er r r
r
r
(a) 7

j

分ければ ，各色ごとの点過緩和を並列処理することができる．

点高次差分式r r r
r er r
r r r

HH j
i
H
H

偶数

奇数

偶数

シアン マゼンタ奇数イエロー(b) 9 点高次差分式
図 3.3: 2 次元高次差分演算子

黒

j

k

連立 1 次方程式の解法

3.2

5

Chebyshev SLOR 法

Chebyshev SLOR 法 (zebra line overrelaxation) は線過緩和法 (SLOR iterative scheme) をベクトルコン
ピュータ向きに改良したものである．5 点差分方程式 (3.1) の場合には，偶数列の線緩和は互いに独立に行えるので，すべての偶数列を並べて同時に線緩和することができる．奇数列についても同様のことが言える．Chebyshev SLOR 法は偶数列の線過緩和と奇数列の線過緩和の並列処理を交互に解が収束するまで反復するものである．線過緩和については 2.4 節の ADI 法の項を参照されたい．

次に，Chebyshev SLOR 法のプログラムの簡単な例として，前節に示したポアソン方程式の Dirichlet 問

題を解くサブルーチンを示す．
!******************************************************************************! Solution of Dirichlet problem of elliptic pde by zebra line overrelaxation!******************************************************************************SUBROUTINE ZEBRA1(u,c,f,r,i1,if,j1,jf,a,b,n,ie,io,bt,na,naf,rmean,e)
DIMENSION u(0:if,0:jf),c(0:if,0:jf,0:4),f(0:if,0:jf),r(0:if,0:jf),a(io,n,3),b(io,n)IF(bt<.75)bt=.8 naf=MAX0(naf,1000)
e=AMAX1(e,1.E-5)
na=0 10 na=na+1
DO ii=1,ie i=2*ii DO j=0,jf!even columnsa(ii,j+1,1) = c(i,j,3)!setup coef matrixa(ii,j+1,2) = c(i,j,0)*bt
a(ii,j+1,3) = c(i,j,4)b(ii,j+1) = f(i,j)-c(i,j,1)*u(i-1,j)-c(i,j,0)*u(i,j) &!setup rhs*(1.-bt)-c(i,j,2)*u(i+1,j)
ENDDO ENDDOCALL GAUSSZ(a,b,io,n,ie)!solve sets of linear eqns
DO ii=1,ie i=2*ii DO j=0,jfu(i,j) = b(ii,j+1)
ENDDO ENDDO
DO ii=1,io i=2*ii-1 DO j=0,jf!odd columnsa(ii,j+1,1) = c(i,j,3)!setup coef matrixa(ii,j+1,2) = c(i,j,0)*bt
a(ii,j+1,3) = c(i,j,4)b(ii,j+1) = f(i,j)-c(i,j,1)*u(i-1,j)-c(i,j,0)*u(i,j) &!setup rhs*(1.-bt)-c(i,j,2)*u(i+1,j)
ENDDO ENDDOCALL GAUSSZ(a,b,io,n,io)!solve sets of linear eqns
DO ii=1,io i=2*ii-1 DO j=0,jfu(i,j) = b(ii,j+1)
ENDDO ENDDO
rmean=0.
DO i=1,i1 DO j=1,j1
r(i,j) = c(i,j,1)*u(i-1,j)+c(i,j,2)*u(i+1,j) &!compute residuals+c(i,j,3)*u(i,j-1)+c(i,j,4)*u(i,j+1) &+c(i,j,0)*u(i,j)-f(i,j)
rmean = rmean+ABS(r(i,j))/FLOAT(i1*j1)!compute mean residual
ENDDO ENDDO
IF(rmean>e.AND.na<naf)
GOTO 10!decite convergence
END!***** Subroutine GAUSSZ *****! Solution of sets of linear eqns with tri-diagonal matrix by Gaussian elimination
SUBROUTINE GAUSSZ(a,b,m,n,if)
DIMENSION a(m,n,3),b(m,n)
DO k=1,n-1 DO i=1,if

6

反法（ Chebyshev SLOR 法）

復b(i,k)=b(i,k)/a(i,k,2) a(i,k,3)=a(i,k,3)/a(i,k,2)b(i,k+1)=b(i,k+1)-a(i,k+1,1)*b(i,k)
a(i,k+1,2)=a(i,k+1,2)-a(i,k+1,1)*a(i,k,3)
ENDDO ENDDOFORALL(i=1:if)b(i,n)=b(i,n)/a(i,n,2)
DO k=n-1,1,-1 DO i=1,ifb(i,k)=b(i,k)-a(i,k,3)*b(i,k+1)
ENDDO ENDDO
END

このプログラムは偶数列の線過緩和，奇数列の線過緩和，平均残差値による収束判定を繰返し実行するものである．このプログラムでは，i1 = if

;1

1 = jf ; 1，ie,io はそれぞれ偶数列の数 ieven または奇数列

 jの数 iodd で，ieven + iodd = i1 である．この場合の iodd は ieven と同数か１つ多くなる．a と b は iodd 個の３重対角行列の連立 1 次方程式の係数行列と右辺を記憶する配列で，n は各連立 1 次方程式の未知変数の数 n = jf +1 である．偶数列の計算では，ieven 個の n 元連立１次方程式 (2.23a) の係数と右辺を作業用配列
a,b に並べて入力し ，サブルーチン GAUSSZ を用いて並べて解く．奇数列の計算も同様である．GAUSSZ は，

1.4 節に述べた１つの３重対角行列の連立 1 次方程式を解くサブルーチン GAUSS3 を，if 個の方程式を並べて解くように拡張したものである．並列度を上げるためには，i に関する DO 文を元のサブルーチンの DO
文の内側に入れなければならない．
問2

2.4 節の混合境界値問題の解を Chebyshev SLOR 法で求めよ．

この混合境界値問題のプログラムを上記サブルーチン ZEBRA1 に倣って作ったものを以下に示す．
!*********************************************************************************! Sol of mixed bound value problem of elliptic pde by zebra line overrelaxation!*********************************************************************************SUBROUTINE ZEBRA2(u,c,f,r,i0,i1,if,j1,jf,a,b,m,n,bt,na,naf,rmean,e)
REAL*4 u(0:if,0:jf),c(0:if,0:jf,0:4),f(0:if,0:jf),r(0:if,0:jf),a(m,n,3),b(m,n)
IF(bt<.75)bt=.8 naf=MAX0(naf,1000)
e=AMAX1(e,1.E-5)
ie0=(i0+1)/2 io0=i0-ie0 ie=if/2+1 io=if+1-ie
je=jf/2+1 jo=jf+1-je
na=0 10 na=na+1
DO ii=2,ie i=2*ii-2 ip1=i+1 DO j=0,jf!even columns
IF(i==if)ip1=if-2!ibn=2a(ii,j+1,1) = c(i,j,3)!setup coef matricesa(ii,j+1,2) = c(i,j,0)*bt
a(ii,j+1,3) = c(i,j,4)b(ii,j+1) = f(i,j)-c(i,j,1)*u(i-1,j)-c(i,j,0)*u(i,j) &!setup rhs*(1.-bt)-c(i,j,2)*u(ip1,j)
ENDDO ENDDOCALL GAUSSMZ(a,b,m,n,2,ie0,jf+1,3,3)!solve sets of linear eqns
CALL GAUSSMZ(a,b,m,n,ie0+1,ie,jf+1,1,1)
DO ii=2,ie i=2*ii-2 DO j=0,jfu(i,j) = b(ii,j+1)
ENDDO ENDDO
DO ii=1,io i=2*ii-1 ip1=i+1 DO j=0,jf!odd columns
IF(i==if)ip1=if-2!ibn=2a(ii,j+1,1) = c(i,j,3)!setup coef matricesa(ii,j+1,2) = c(i,j,0)*bt
a(ii,j+1,3) = c(i,j,4)b(ii,j+1) = f(i,j)-c(i,j,1)*u(i-1,j)-c(i,j,0)*u(i,j) &!setup rhs*(1.-bt)-c(i,j,2)*u(ip1,j)
ENDDO ENDDO

連立 1 次方程式の解法
CALL GAUSSMZ(a,b,m,n,1,io0,jf+1,3,3)
CALL GAUSSMZ(a,b,m,n,io0+1,io,jf+1,1,1)
DO ii=1,io i=2*ii-1 DO j=0,jfu(i,j) = b(ii,j+1)
ENDDO ENDDO
GOTO 11na=na+1
DO jj=1,je j=2*(jj-1) jm1=j-1 jp1=j+1 DO i=0,if
IF(i<i0)THEN IF(j==0)jm1=jf-1 IF(j==jf)jp1=1
ELSE
IF(j==0)jm1=0 IF(j==jf)jp1=jf
ENDIF
a(jj,i+1,1) = c(i,j,1)a(jj,i+1,2) = c(i,j,0)*bt
a(jj,i+1,3) = c(i,j,2)b(jj,i+1) = f(i,j)-c(i,j,3)*u(i,jm1)-c(i,j,0)*u(i,j) &*(1.-bt)-c(i,j,4)*u(i,jp1)
ENDDO ENDDOCALL GAUSSMZ(a,b,m,n,1,je,if+1,1,2)
DO jj=1,je j=2*(jj-1)
FORALL(i=0:if)u(i,j)=b(jj,i+1)
ENDDO
DO jj=1,jo j=2*jj-1 DO i=0,ifa(jj,i+1,1) = c(i,j,1)a(jj,i+1,2) = c(i,j,0)*bt
a(jj,i+1,3) = c(i,j,2)b(jj,i+1) = f(i,j)-c(i,j,3)*u(i,j-1)-c(i,j,0)*u(i,j) &*(1.-bt)-c(i,j,4)*u(i,j+1)
ENDDO ENDDOCALL GAUSSMZ(a,b,m,n,1,jo,if+1,1,2)
DO jj=1,jo j=2*jj-1
FORALL(i=0:if)u(i,j)=b(jj,i+1)
ENDDO
11 rmean=0.
DO i=1,if-1 im1=i-1 ip1=i+1 DO j=0,jf jm1=j-1 jp1=j+1
IF(i<i0)THEN IF(j==0)jm1=jf-1 IF(j==jf)jp1=1
ELSE
IF(j==0)jm1=0 IF(j==jf)jp1=jf
ENDIF
IF(i==if.AND.j>0.AND.j<jf)ip1=if-2
r(i,j) = c(i,j,1)*u(im1,j)+c(i,j,2)*u(ip1,j) &+c(i,j,3)*u(i,jm1)+c(i,j,4)*u(i,jp1) &+c(i,j,0)*u(i,j)-f(i,j)
rmean=rmean+ABS(r(i,j))/FLOAT(if*jf)
ENDDO ENDDO
IF(rmean>e.AND.na<naf)
GOTO 10END

7!solve sets of linear eqns

!even rows
!ib1=ibn=3!ib1,ibn=1!setup coef matrices!setup rhs!solve sets of linear eqns

!odd rows!setup coef matrices!setup rhs!solve sets of linear eqns

!ib1=ibn=3!ib1,ibn=1!ibn=2!compute residuals!compute mean residual!decide convergence!***** Subroutine GAUSSMZ *****! Sol of sets of linear eqns with modified tri-diagonal matrix based on Gaussian elimination
SUBROUTINE GAUSSMZ(a,b,m1,n1,is,ie,n,ib1,ibn)
REAL*4 a(m1,n1,3),b(m1,n1)
IF(ib1==2)THEN!ib1=2
FORALL(i=is:ie)a(i,1,1)=a(i,1,1)/a(i,1,2)
FORALL(i=is:ie)a(i,2,3)=a(i,2,3)-a(i,2,1)*a(i,1,1)
ENDIF
DO k=1,n-1IF(ib1==3.AND.k==n-2) &!ib1=ibn=3
FORALL(i=is:ie)a(i,k,3)=a(i,k,3)+a(i,k,1)

8

反

復

法（ Chebyshev SLOR 法）

FORALL(i=is:ie)b(i,k)=b(i,k)/a(i,k,2) a(i,k,3)=a(i,k,3)/a(i,k,2)b(i,k+1)=b(i,k+1)-a(i,k+1,1)*b(i,k)
a(i,k+1,2)=a(i,k+1,2)-a(i,k+1,1)*a(i,k,3)
ENDFORALL
IF(ibn==2.AND.k==n-2)THENFORALL(i=is:ie)b(i,n)=b(i,n)-a(i,n,3)*b(i,k)
FORALL(i=is:ie)a(i,n,1)=a(i,n,1)-a(i,n,3)*a(i,k,3)
ENDIF
IF(ib1==3)THEN
IF(k<n-2)THENFORALL(i=is:ie)a(i,k,1)=a(i,k,1)/a(i,k,2)
FORALL(i=is:ie)a(i,k+1,1)=-a(i,k+1,1)*a(i,k,1)
ENDIF
IF(k>1.AND.k<n-2)THENFORALL(i=is:ie)b(i,n)=b(i,n)-a(i,n,3)*b(i,k)
FORALL(i=is:ie)a(i,n,1)=a(i,n,1)-a(i,n,3)*a(i,k,1)
FORALL(i=is:ie)a(i,n,3)=-a(i,n,3)*a(i,k,3)
ELSEIF(k==n-2)THENFORALL(i=is:ie)b(i,n)=b(i,n)-a(i,n,3)*b(i,k)
FORALL(i=is:ie)a(i,n,1)=a(i,n,1)-a(i,n,3)*a(i,k,3)
ENDIF
ENDIF
ENDDO
FORALL(i=is:ie)b(i,n)=b(i,n)/a(i,n,2)
DO k=n-1,1,-1
FORALL(i=is:ie)b(i,k)=b(i,k)-a(i,k,3)*b(i,k+1)
IF(ib1==3.AND.k<n-2) &FORALL(i=is:ie)b(i,k)=b(i,k)-a(i,k,1)*b(i,n-1)
ENDDO
IF(ib1==2)FORALL(i=is:ie)b(i,1)=b(i,1)-a(i,1,1)*b(i,3)
END!ibn=2!ib1=ibn=3!ib1=ibn=3!ib1=2

このプログラムでは，偶数 i 列の線過緩和，奇数 i 列の線過緩和，偶数 j 行の線過緩和，奇数 j 行の線過緩和，平均残差値による収束判定が繰り返される．i 列の線過緩和では，連立 1 次方程式は周期条件の課されている部分と Dirichlet 条件の課されている部分に分けて解かれる．また偶数 j 行の線過緩和でもこれらの境界条件を考慮するなど 必要である．偶数列または奇数列の連立 1 次方程式はサブルーチン GAUSSMZ を用い並列処理される．GAUSSMZ は 2.4 節のサブルーチン GAUSSM を並列処理できるように書換えたもので，ここに出てくるすべての境界条件に対応できる．その引数 ie0,io0,ie,io,je,jo は周期条件の課されている区間の偶数列の数 ieven0 ，奇数列の数 iodd0 ，全緩和区間の偶数列の数 ieven ，奇数列の数 iodd ，偶数行の数 jeven ，奇数行の数 jodd である．また m1,n1 は同時に解く連立 1 次方程式の最大数，その未知数の最大数，is,ie は同時に線緩和する初めと終わりの列番号または行番号，ib1,ibn については 2.4 節参照，メ
インプログラムについても 2.4 節参照．

以下には，ポアソン方程式 uxx + uyy =;2 の Dirichlet 問題とラプラス方程式値問題を，図 3.4 の境界条件のもとで出発値 u= 0 の混合境界= 0 を与えて解いた結果を示す．計算には格子間隔 0.05 の正方形格子を用いた．次表は " = 10;5 における係数

uxx

+

uyy

と最終反復数 n の関係を示したものである．なお収束しないときには nf = 1000 における平均残差値 rmean が示してある．i-column は i 列の線過緩和のみの場合，j-row は j 行の線過緩和のみの場合である．i 列の線過緩和と j 行の線過緩和を交互に行う ADI の場

合は，予想に反し反復数は多くなった．括弧内の数字は i 列と j 行の線過緩和を 4 回ずつ交互に行った場合

で反復数は少なくなる．この表から i 列の線過緩和または j 行の線過緩和のど ちらかで十分なことが分かる．

Chebyshev SLOR 法は，図 3.3(b) に示す２次元の 9 点高次差分方程式に対してはほとんどそのまま適用

連立 1 次方程式の解法

Dirichlet
u=1

1

1

periodicu(x1)= u(x0)+1
u (x1)= u (x0)Dirichlet
u = y Dirichlet
u=y

y =x0= 0

Dirichletu=0(a) Dirichlet 問題y

HH

Neumannu = ;1yA
AAA

y =x0= ;1

1

Dirichlet
u = 1;x

HH

Dirichlet
u=y

9

0

x

Dirichletu = ;x(b) 混合境界値問題

1

図 3.4: Chebyshev SLOR 法の例題
表 3.2: Chebyshev SLOR 法の解の収束状況

0.75
0.80
1.00
1.25
1.50

Dirichlet 問題
混合境界値問題i-column
i-column j-row
ADI391
41  10;5 diverge
260 (104)24
99
100
306 (152)
157
426
430
536 (464)
313
853
856
926 (888)
;
5
;
5
470
4  104  105  10;5(5  10;5)

することができる．(a) 図の 9 点高次差分方程式の場合には， mod (i 3) が 0 の列，1 の列，2 の列に分け線緩和を並列処理するように書き変えることが必要である．３次元の 7 点差分方程式の場合には，xy 面
内の Chebyshev SLOR 法によるの計算を 方向に掃引することも考えられるが，Chebyshev SLOR 法とChebyshev SOR 法を組み合わせ， 方向に並ぶ点列を + が偶数のものと奇数のものに分け線緩和を並列zzに行えばいっそう高並列度になる．

i

j

10

反

3.3

復

法（ 多重格子法）

多重格子法

多重格子 (multigrid，MG) 法は，スカラーコンピュータの時代に，楕円型偏微分方程式の境界値問題の効率的な数値解法として考案されたもので，その後多くの問題に適用できるように拡張されてきた．多重格子法は，格子間隔 h の基本格子のほかに，それに重なる格子間隔 2h 4h : : : の格子を用い緩和計算を行うものである．前述のように，短波長の残差の波は容易に緩和できるが，長波長のものは容易でない．多重格子法は，異なる寸法の格子を用い，すべての波長の残差を効率よく緩和しようとするものである．多重格子法では，格子間隔 h の細かい格子の残差は限定集計 (restriction) によって１つずつ粗い格子のものに移され，粗い格子の修正値は延長補間 (prolongation) によって１つずつ細かい格子のものに移される．また
点緩和は適宜 各格子で行われるが，過緩和の必要はなく専ら Gauss 反復法が用いられる．最も粗い格子では直接法も用いられる．限定集計，延長補間，Gauss 反復法はすべて並列処理可能で多重格子法はベクトル
コンピュータにも適する計算技法である．

3.3.1

楕円型線形微分方程式の多重格子法;2 に対して多重格子法を説明する．格子間隔 の正方形格子の場合

始めに Poisson 方程式 uxx + uyy =hには，差分方程式はui;1 + u +1 + u ;1 + u +1 ; 4u = ;2h2
j

i

j

ij

ij

ij

となり，その残差の式は
rij= ;1 + +1 +
ui

j

ui

j

uij

;1 + u +1 ; 4u +2h2ijijとなる．今出発値 uij = 0 を与えて計算を始めることにすればその初期残差は rij = 2h2 となる．また格子
間隔 H = 2h の粗い格子に対しては初期残差は rIJ = 2H 2 = 8h2 となる．これより初期残差値は，格子セルの面積に比例することが分かる．格子間隔
ui

x

=

h

y

=

h=

の長方形格子の場合には，差分方程式は;1 ; 2u + u +1 + 2 (u ;1 ; 2u + u +1) = ;2h2
j

ij

i

j

ij

ij

(3.2)

ij

また残差の式はrij= ;1 ; 2 + +1 + 2 (
ui

j

uij

ui

j

uij

;1 ; 2u + u +1 )+2h2
ij

(3.3)

ij

となる．多重格子法では，まず既知の h 格子の残差値 rij に対し H 格子の残差値 rIJ を見積もることが必要である．
上記の初期残差の例では，残差値が計算領域を通して一様で，また計算領域内の残差値を集計したものは格子の大きさによらず一定であった．この性質を一般化すれば h 格子と H 格子の残差の関係は次のようにな

る．図 3.5(a) は多重格子を示したもので，h 格子の格子点を 印，H 格子の格子点を印で，また h 格子の

格子点番号を ij ，H 格子のものを IJ で示してある．計算領域内の残差の集計値が一定になるように，残差値 rij はすべて rIJ に，ri+1j は 2 分して rIJ と rI +1J に，ri+1j +1 は 4 分して rIJ  rI +1J  rI +1J +1  rI J +1に集計することにする．このとき，残差値 rIJ は図に

 印を付けた格子点の残差値

rij  : : :を太線枠内の面積を勘案?して集計したものになる．このような操作を限定集計 (restriction) と呼び，式で表せば次のよう

になる．
rIJ

=

rij

+ 12 ( ;1 + +1 +
ri

j

ri

j

rij

1

;1 + r +1 )+ 4 (r ;1 ;1 + r +1 ;1 + r ;1 +1 + r +1 +1 )
ij

i

j

i

j

i

j

i

j

(3.4)

連立 1 次方程式の解法

j +2
J +1
j +1

e

j
J

j ;1
j ;2
J ;1

i;2
I ;1

s

s

s

s

se

s

s

s

s

i;1

             
a/4

e

a/4

-a/2

-1-a
-2-2a
-1/2

-a/2

a/4

a/2
a

a/4

1/4
1/2
1
1/4

e

i +1

i
I
(a)

a/2
a
-1/2

i +2
I +1

図 3.5: 長方形格子の限定集計

(b)

11

1/4
1/2
1
1/4

a= 2

式 (3.3) を限定集計すれば ，左辺はかなり長い式になるが，その係数を整理したものは (b) 図のようになる．

すなわち，この図の位置関係は (a) 図と同じで，図中の円またはオーバル内の上段の数字になる．また下段の数字はこの係数を H 格子の格子点に集めたものである．これより h 格子のラプラス演算子を限定集計したものは，当然のことながら，H 格子のラプ ラス演算子で近似できることが分かる．結局 式 (3.3) を限定集計した H 格子の式は次のようになる2 ．

= ;1 ; 2

rIJ

uI

uIJ

J

+

uI

2 (u

+1J +I J;1 ; 2u + u +1)+2H 2
IJ

I J

(3.5)

多重格子法の計算の途中の段階を考え，その u の近似値を uij ，また解を uij + dij とする．この dij は理想的な修正値で，すべての残差を 0 にするものである．したがって

0 = ( ;1 + ;1 ) ; 2( + )+( +1 + +1 )
+ 2 ( ;1 + +1 ) ; 2( + )+( +1 + +1 ) +2 2
ui

j

uij

di

j

dij

uij

dij

ui

uij

dij

uij

j

di

j

dij

h

この式に残差の式 (3.3) を用いれば次の h 格子の修正値 dij の差分方程式が導かれる．

;1 ; 2d + d +1 + 2 (d +1 ; 2d + d +1 )+2h2 = ;r

di

j

ij

i

j

ij

ij

ij

(3.6)

ij

この式を限定集計すれば ，次の H 格子の修正値 dIJ の差分方程式が導かれる．
dI;1 ; 2d + d +1 + 2 (d +1 ; 2d + d +1)+2h2 = ;r
J

IJ

I

J

I J

IJ

I J

IJ

(3.7)

通常この差分方程式はガウス反復法で解かれ H 格子の uIJ の修正値 dIJ が求められる．この修正値は h 格子の修正値 dij にそのまままたは Lagrange 補間によって移される．

== (; ;3 +9
2 = (3 1 +6 3 ; 5 ) 8

dij

dIJ 

di

di

dij

dij

di

di

dij

;1 +9d +1 ; d +3 )=16
ij

ij

=

2 なお差分方程式と残差の式がu ;1 ; 2u + u +1 + u ;1 ; 2u + u +1 = ;2
x2
y2
u
u
u + u +1 +2;
1 ; 2u + u +1
r =+ ;1 ; 2y
2
x2
i

ij

j

ij

i

j

i

ij

j

ij

i

j

ij

ij

ij

ij

ij

のようになっている場合には，残差の限定集計ではなく重み平均を取るべきである．

(3.8)

12

反

復

法（ 多重格子法）

多重格子法では H 格子の修正値からこのようにして h 格子の修正値を求める操作を延長補間 (prolongation)という．多重格子法では，各格子レベルで緩和計算を行い，格子サイズに合った残差の波の成分を緩和することになる．小サイズの格子の修正値は，その格子で求めた修正値にそれよりも大きいサイズの格子で求めた修正値を延長補間したものすべてを合算したものになる．
以上の説明をもとに，次に多重格子法を一般的に説明する．線形微分方程式は

( )=

L u

(3.9)

f

またその近似差分方程式は
h

L

( )=
h

u

f

(3.10)

hのように書くことができる．上添え字 h は格子間隔を意味するものとする．その第 n 近似解を uh(n) とす

れば ，式 (3.10) の残差は
r( ) = Lh (uh(n) ) ; f hh n(3.11)となる．今 式 (3.10) の残差を零にする理想的な uh(n) の修正値を dh とする．すなわち Lh (uh(n)+dh ) = f h ．
hLは線形演算子であるから，この式と式 (3.11) から，次の細かい格子の修正値 dh の差分方程式が得られる．
hL( )=; ( )
h

d

h n

r

(3.12)

ここで格子間隔 H = 2h の粗い格子を考え，式 (3.12) の限定集計を取れば次式が得られる．
H

( )=;

H

=

L

H

d

r

H

(3.13)

ただし
r

H

Ih

H

Ih r

h

(3.14)

は限定集計演算子 (restriction operator) と呼ばれるものである．式 (3.13) は H 格子の修正値 dH の差分方程式で，dH の値は緩和計算によって求められる．修正値 dH から h 格子の修正値 dh は延長補間によって求められる．
d

h

=

h

H

IH d

(3.15)h
ただし IHは延長補間演算子 (prolongation operator) と呼ばれるものである．u の第 n +1 近似値は

( +1) = uh(n) + d~hh nu(3.16)となる．ただし d~h は，h 格子の緩和計算で求めた修正値に，2h 4h : : : 格子の緩和計算で求めた修正値を延長補間したものを加えた修正値である．
多重格子法の計算の手順，すなわち残差の計算，限定集計，修正値の計算，延長補間，修正値の加算をどのような順番で反復するのかは各種のものが提案されている．その最適手順は解くべき問題，微分方程式と境界条件，格子や数値解法によるものと思われる．図 3.6 に２つの 1 サイクルの計算手順を示す．この図の説明は次項にする．

mesh 1mesh 2mesh 3

e

er

;?

;

?
?


;
;


;

;
mesh 4 ?

;?
;

連立 1 次方程式の解法

er r e r e



;
;?

;
;


?
?




;
;




;
?
;



(a) MGM1eee r e determine residual



;

;


(b) MGM2

13

restriction
 determine correction
 prolongation;;
apply correction?

r

図 3.6: 多重格子法の計算手順

3.3.2

Poisson 方程式の Dirichlet 問題の多重格子法のプログラム

ここには多重格子法のプログラムの一例を示す．次の Poisson 方程式の Dirichlet 問題

= ;20
(0
2 0
(0 ) = (2 ) = ( 0) = ( 1 2) = 0

uxx
u

+

 y

uyy

< x <

u

 y

u x

u x



< y <

1 2)

(3.17a)(3.17b)::を考える．計算領域は長方形で，計算格子は４重の正方形格子でそれぞれを mesh 1, mesh 2, mesh 3, mesh

4 と呼ぶことにする．mesh 1 は格子間隔 = 0 05 の基本格子である．プログラム中の記号の意味はだいたh:

い分かると思うが，各格子の残差と修正値は次のように定義する．

mesh 1：(41  25)
mesh 2：(21  13)
mesh 3： (11  7)
mesh 4： (6  4)

R1 = r1 ,D1 = d1

R3 = r3 ,D3 = d3

R2 = r2 ,R4 = r4 ,D2 = d2
D4 = d4

また a と b は，修正値 d4 を Gauss 消去法で計算するために連立 1 次方程式の係数行列と右辺を入れる配列である．この問題は対称であるが，ここでは全領域を計算することにする．サブルーチン MGM1 は図 3.6(a)に示す標準的な計算手順，また MGM2 は (b) 図に示す簡略化した手順のプログラムである．
program MAIN
REAL u(0:40,0:24),f(0:40,0:24)
DATA u/1025*0./ f/1025*-20./
DATA i1,j1,h/40, 24, .05/ naf,e/100, 1.E-5/
OPEN(20,FILE='OUTPUT.dat')
CALL MGM2(u,f,i1,j1,h,na,naf,rmean,e)
WRITE(20,'(I5, F6.3)')(u(i,j),i=0,i1)na,rmean
DO j=j1,0,-1
WRITE(20,'(41F6.3)')(u(i,j),i=0,i1)
ENDDO
CLOSE(20)STOP
END!********** Subroutine MGM1
SUBROUTINE MGM1(u,f,i1,j1,h,na,naf,rmean,e)
REAL u(0:i1,0:j1),f(0:i1,0:j1),a(24,-4:4),b(24), &r1(0:40,0:24),r2(0:20,0:12),r3(0:10,0:6),r4(0:5,0:3), &d1(0:40,0:24),d2(0:20,0:12),d3(0:10,0:6),d4(0:5,0:3)
i2=i1/2 i3=i2/2 i4=i3/2
j2=j1/2 j3=j2/2 j4=j3/2
kf=(i4+1)*(j4+1) lf=j4+1
na=0 100 na=na+1
n =5 101 n=n-1! Determine residuals and restrictions!predictor+bound values, rhs!call MGM1 or MGM2
!output na,rmean!output u

14

反

復

法（ 多重格子法）

FORALL(i=1:i1-1,j=1:j1-1)r1(i,j) = &u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1)-4.*u(i,j)-h*h*f(i,j) !comp r1
IF(n==1)GOTO 11
CALL RESTRC(r1,i1,j1,r2,i2,j2)
IF(n==2)GOTO 12!restriction(r1 → r2)
CALL RESTRC(r2,i2,j2,r3,i3,j3)
IF(n==3)GOTO 13!restriction(r2 → r3)
CALL RESTRC(r3,i3,j3,r4,i4,j4)
IF(n==4)GOTO 14!restriction(r3 → r4)! Determin corrections by point relaxation
11 FORALL(i=1:i1-1,j=1:j1-1)d1(i,j)=r1(i,j)/4. GOTO 21!comp d1
12 FORALL(i=1:i2-1,j=1:j2-1)d2(i,j)=r2(i,j)/4. GOTO 22!comp d2
13 FORALL(i=1:i3-1,j=1:j3-1)d3(i,j)=r3(i,j)/4. GOTO 23!comp d3
! Determine corrections by Gaussian elimination
14 DO i=0,i4 DO j=0,j4 k=i*lf+j+1!comp d4
FORALL(l=-lf:lf)a(k,l)=0. a(k,0)=1.
IF(i>0.AND.i<i1.AND.j>0.AND.j<j1)THENa(k,-lf)=1. a(k,-1)=1. a(k,0)=-4. a(k,1)=1. a(k,lf)=1.
ENDIF
b(k)=-r4(i,j)
ENDDO ENDDOCALL GAUSSB(a,b,kf,-lf,lf,det)
FORALL(i=0:i4,j=0:j4)d4(i,j) = b(i*lf+j+1)! Prolongations and apply corrections
CALL PROLON(d4,i4,j4,d3,i3,j3)!prolongation(d4 → d3)
23 CALL PROLON(d3,i3,j3,d2,i2,j2)!prolongation(d3 → d2)
22 CALL PROLON(d2,i2,j2,d1,i1,j1)!prolongation(d2 → d1)
21 FORALL(i=1:i1-1,j=1:j1-1)u(i,j)=u(i,j)+d1(i,j)!apply d1! Decide convergence by mean residual
rmean=0.
DO i=1,i1-1 DO j=1,j1-1
rmean = rmean+ABS(r1(i,j))!comp mean residual
ENDDO ENDDO
rmean = rmean/FLOAT(i1*j1)
IF(n>1)GOTO 101IF(rmean>e.AND.na<naf)
GOTO 100!decide convergence
END!***** Subroutine RESTRCSUBROUTINE RESTRC(r1,i1,j1,r2,i2,j2)
REAL r1(0:i1,0:j1),r2(0:i2,0:j2)
DO ii=1,i2-1 i=2*ii DO jj=1,j2-1 j=2*jj
r2(ii,jj) = r1(i,j)+(r1(i-1,j)+r1(i+1,j)+r1(i,j-1)+r1(i,j+1))/2. &+(r1(i-1,j+1)+r1(i+1,j-1)+r1(i-1,j-1)+r1(i+1,j+1))/4.
ENDDO ENDDO
END!***** Subroutine PROLONSUBROUTINE PROLON(d2,i2,j2,d1,i1,j1)
REAL d2(0:i2,0:j2),d1(0:i1,0:j1)
FORALL(i=1:i2-1,j=1:j2-1)d1(2*i,2*j)=d2(i,j)
DO ii=1,i2-1 i=2*ii
d1(i,1)=(3.*d1(i,0)+6.*d1(i,2)-d1(i,4))/8.
FORALL(j=3:j1-3:2)d1(i,j)=(-d1(i,j-3)+9.*d1(i,j-1)+9.*d1(i,j+1)-d1(i,j+3))/16.
d1(i,j1-1)=(-d1(i,j1-4)+6.*d1(i,j1-2)+3.*d1(i,j1))/8.
ENDDO
DO j=1,j1-1
d1(1,j)=(3.*d1(0,j)+6.*d1(2,j)-d1(4,j))/8.
FORALL(i=3:i1-3:2)d1(i,j)=(-d1(i-3,j)+9.*d1(i-1,j)+9.*d1(i+1,j)-d1(i+3,j))/16.
d1(i1-1,j)=(-d1(i1-4,j)+6.*d1(i1-2,j)+3.*d1(i1,j))/8.
ENDDO
END

連立 1 次方程式の解法!***** Subroutine GAUSSBSUBROUTINE GAUSSB(a,b,n,m1,m2,det)
DIMENSION a(n,m1:m2),b(n)det=1.
cycle_1: DO k=1,ndet=a(k,0)*detb(k)=b(k)/a(k,0) IF(k==n)CYCLE cycle_1
j2=MIN0(m2,n-k)
FORALL(j=1:j2)a(k,j)=a(k,j)/a(k,0)
i2=MIN0(n,k-m1)
DO i=k+1,i2 j0=k-ib(i)=b(i)-a(i,j0)*b(k)DO jk=1,j2 j=jk+k-ia(i,j)=a(i,j)-a(i,j0)*a(k,jk)
ENDDO
ENDDO
ENDDO cycle_1
DO i=n-1,1,-1 j2=MIN0(m2,n-i)
DO j=1,j2 k=i+jb(i)=b(i)-a(i,j)*b(k)
ENDDO
ENDDO
END
!********** Subroutine MGM2
SUBROUTINE MGM2(u,f,i1,j1,h,na,naf,rmean,e)
REAL u(0:i1,0:J1),f(0:i1,0:j1),a(24,-4:4),b(24), &r1(0:40,0:24),r2(0:20,0:12),r3(0:10,0:6),r4(0:5,0:3), &d1(0:40,0:24),d2(0:20,0:12),d3(0:10,0:6),d4(0:5,0:3)
i2=i1/2 i3=i2/2 i4=i3/2
j2=j1/2 j3=j2/2 j4=j3/2
kf=(i4+1)*(j4+1) lf=j4+1
na=0 100 na=na+1! Determine residuals and restrictions
FORALL(i=1:i1-1,j=1:j1-1)r1(i,j) = &u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1)-4.*u(i,j)-h*h*f(i,j)
CALL RESTRC(r1,i1,j1,r2,i2,j2)
CALL RESTRC(r2,i2,j2,r3,i3,j3)
CALL RESTRC(r3,i3,j3,r4,i4,j4)! Determine corrections by Gaussian elimination
DO i=0,i4 DO j=0,j4 k=i*lf+j+1
FORALL(l=-lf:lf)a(k,l)=0. a(k,0)=1.
IF(i>0.AND.i<i1.AND.j>0.AND.j<j1)THENa(k,-lf)=1. a(k,-1)=1. a(k,0)=-4. a(k,1)=1. a(k,lf)=1.
ENDIF
b(k)=-r4(i,j)
ENDDO ENDDOCALL GAUSSB(a,b,kf,-lf,lf,det)
FORALL(i=0:i4,j=0:j4)d4(i,j) = b(i*lf+j+1)! Prolongation, corect corrections by point relaxation
CALL PRORLX(d4,i4,j4,d3,i3,j3,r3)!prolongation(d4 → d3),d3
CALL PRORLX(d3,i3,j3,d2,i2,j2,r2)!prolongation(d3 → d2),d2
CALL PRORLX(d2,i2,j2,d1,i1,j1,r1)!prolongation(d2 → d1),d1
FORALL(i=1:i1-1,j=1:j1-1)u(i,j)=u(i,j)+d1(i,j)!apply d1! Decide convergencermean=0.
DO i=1,i1-1 DO j=1,j1-1
rmean = rmean+ABS(r1(i,j))

15

16

反

復

法（ 多重格子法）

ENDDO ENDDO
rmean = rmean/FLOAT(i1*j1)
IF(rmean>e.AND.na<naf)
GOTO 100
END!***** Subroutine PRORLXSUBROUTINE PRORLX(d2,i2,j2,d1,i1,j1,r1)
REAL d2(0:i2,0:j2),d1(0:i1,0:j1),r1(0:i1,0:j1)

この間にサブルーチン PROLON と同じ延長補間に関する 11 行の文が入る]
FORALL(i=1:i1-1,j=1:j1-1) &r1(i,j)=r1(i,j)+d1(i-1,j)+d1(i+1,j)+d1(i,j-1)+d1(i,j+1)-4.*d1(i,j)
FORALL(i=1:i1-1,j=1:j1-1)d1(i,j)=d1(i,j)+r1(i,j)/4.
END

メインプログラムではサブルーチン MGM1, MGM2 のいずれかを CALL しなければならない．サブルーチン

MGM1 は，図 3.6(a) に示す計算手順によるもので，下図のフローチャートのように組まれている．一般に緩和計算では残差の計算と修正値の計算は対をなすが，MGM1 では常に基本格子で残差 r1 を計算し，限定集計

によって r2  r3  r4 を求めている．また各格子で修正値 d4  d3  d2 を計算し延長補間によって基本格子の d1を求め解 u に加えている．最後に r1 から直接 d1 求め解 u に加えている．n は緩和計算を行う格子を示す指標である．図 (a) の手順は解が十分に収束するまで反復される．na はサイクル数で，計算は残差の平均

値が十分小さくなったとき rmean"，または所定のサイクル数に達したとき na = naf に終了する．プログラムの内容はその中のコメントと下図から理解できるものと思う．

START
n =0

?

n = n +1 
; Convergent ? @
;
@ Yes
n=5
@
;
@ or n = n ? ;YesHHn = n;1 ?
H
Hn > 1 ?
STOP
Compute r1
Apply d1
6
H
H Yes- Compute d1
H
Hn = 1 ?
Prolongation
Restriction
6
H
H Yes- Compute d2
H
Hn = 2 ?
Prolongation
Restriction
6
H
H Yes- Compute d3
H
Hn = 3 ?
Prolongation
Restriction- Compute d46
図 3.7: MGM1 のフローチャート
a

a

a

a

af

e
e
e

サブルーチン MGM2 はサブルーチン MGM1 の計算手順をコンパクトに纏めたものである．MGM1 では粗い

mesh 4 で計算した修正値 4 は延長補間によって一気に mesh 1 の修正値 1 まで移したが，MGM2 では延長補間によって mesh 3 の修正値 3 に移した時点で，先に限定集計によって求めた残差値 3 がその影響を受
d

d

d

r

けることに注目し ，
rij

=

rij

+ ;1 + +1 +
di

j

di

j

dij

;1 + d +1 ; 4dijij

連立 1 次方程式の解法

17のように修正する．この式の右辺の残差 r3 は修正前の u に対する残差で，また左辺の残差 r3 は修正後の+ 3 に対する残差である．サブルーチン PRORLX では，まず PROLON と同じ延長補間によって 4 から 3を求め，次に 3 を修正し，それからこの修正された 3 に対し点緩和の修正値4 を求めこれを延長補間3で移した に加えている．この一連の操作は mesh 2 と mesh 1 に関しても同様に行われる．MGM2 ではこ
u

d

d

r

r

d

rij =

d

のようにして，１つのループの中で各格子レベルで残差の修正，修正値の計算と加算が効率よく行われることになる．

計算の結果は次表に示すように MGM1 と MGM2 による計算はともにサイクル数 na = 8 程度で収束し ，昔のベクトルコンピュータによる計算では，MGM2 の CPU 時間は MGM1 のものの約 1/2 である．
MGM1

n
1
2
3
4
5
6
7
8

rmean
0.009063
0.002929
0.000608
0.000251
0.000068
0.000032
0.000013
0.000007

a

3.3.3

ms
20
39
57
75
93
111
128
145

MGM2

rmean
0.009046
0.002901
0.000602
0.000248
0.000068
0.000031
0.000014
0.000007

ms
8
17
26
34
41
49
57
65

準線形 2 階楕円型微分方程式の混合境界値問題への適用

式 (3.9) が 2 独立変数 2 階の準線形楕円型微分方程式1の場合について述べる．

( )

L u

auxx

+

buxy

+

cuyy

+ =0

(3.18)d

ただし a b c d は一般に x y u ux  uy の関数である．準線形の場合には，a b c d の中の uh を暫時 定数とみなし ，線形の場合の式をそのまま用いて計算を進め，随時 a b c d の中の uh を新しいものに換えていく方法が取られる．多くの問題では，微分方程式中の最高階数の微分が主役を演じるので，この方法で収束解を得ることができる．

式 (3.18) の近似差分方程式は一般に0+ 1 ;1 + 2 +1 + 3 ;1 + 4+ 5 ;1 ;1 + 6 +1 ;1 + 7 ;1 +1 + 8
cij uij

cij uicij uijjcij uicij uijcij uijcij uijjcij uij+1+1j +1 = fij

cij ui(3.19)のように書かれる．ただし c0ij + c1ij + c2ij + c3ij + c4ij = 0，c5ij + c6ij + c7ij + c8ij = 0 で，左辺の係数は優対角の条

件を満足またはあまり逸脱しないものとする．また右辺 fij は格子間隔の２乗 h2 またはこれに相当の因子を含むものとする．式 (3.19) の場合には残差の式 (3.11) は(n) = c0 u(n) + c1 u(n) + c2 u(n) + c3 u(n) + c4 u(n)ij ijij i;1j
ij i+1jij ij ;1
ij ij +1
(
n)
(
n)
(
n)
n)
5
6
7
+cij ui;1j;1 + cij ui+1j;1 + cij ui;1j+1 + c8ij u(i+1
; fij
j +1

rij

(3.20)

また修正値を求める式 (3.12) は0+ 1 ;1 + 2 +1 + 3 ;1 + 4+ 5 ;1 ;1 + 6 +1 ;1 + 7 ;1 +1 + 8
cij dijcij dicij dijjcij dicij dijjcij dijcij dicij dij

j

+1

(n)
+1j +1 = ;rijcij di(3.21)

1 準線形とは最高階数に関して線形であることをいう．例えば 2 階の微分方程式の場合には 2 階微分の線形結合，1 階の連立微分方程式の場合には各式が 1 階微分の線形結合であることをいう．

18

反

復

法（ 多重格子法）となる．次にこの式に上記の限定集計の式 (3.4) を適用し，粗い格子の修正値を求める式を導出する．その c2  c3 
4
1
3
0
5
8
6
7
2
c  c + c + c =2  0 c  c  ;c  ;c のように近似可能で ，また通常なめらかに変化することを考式は初めかなり冗長なものになる．しかしながら，係数 c は式 (3.18) の 2 階微分に由来し c1

慮すれば ，限定集計した式の dij  di;1j  di;2j  di;1j ;1 の係数は次のように簡単になる．

0 +(c1 + c2 + c3 + c4 )=2+(c5+ c6i;1j+1 + c7i+1j;1 + c8i;1j;1 )=4  c0ij =2
i+1j
i;1j
ij +1ij ;1i+1j +1
0
1
3
4
5
7
0
1
ci;1j =2+ cij +(ci;1j +1 + ci;1j +1 )=4+(cij +1 + cij ;1 )=2  ci;1j =4+ cij =21
571
ci;1j =2+(ci;1j +1 + ci;1j ;1 )=4  ci;1j =2
0
1
3
5
5
ci;1j ;1 =4+(cij ;1 + ci;1j )=2+ cij  cij
cij;2 ;1 と d ;2 ;2 の係数はそれぞれ c1;1 ;1 =4+ c5;1 =2，c5;1 ;1 =4 となる．ほかの係数についても同

di

j

i

j

i

j

i

j

i

j

様である．また延長補間で得られた修正値 dh がかなりなめらかに変化することに注意すれば ，近似的に;1 ;1  c5 d ;1 ;1 =2+(c5 d ;2 ;2 + c5 d )=4555c ;1 d ;2 ;1  (c ;1 d ;2 ;2 + c ;1 d ;2 )=25

cij di
i

j

i

j

i

ij

j

j

i

i

ij

j

i

j

j

i

ij

ij

i

j

j

のように分配し ，また線形の場合と同様に;1 + c0;1 d ;1 + c0;1 d ;1 + c0;1 d ;1 )=4  c0 d 
1
1
1
1
c ;1 d ;2 =2+(c ;1 ;1 d ;2 ;1 + c ;1 +1 d ;2 +1 )=4  c ;1 d ;1 
(c5 + c5;1 + c5 ;1 + c5;1 ;1 )d ;2 ;2  c5;1 ;1 d ;1 ;1

( 0;1
ci

j di

i

j

ij

j

i

j

i

j

i

j

i

i

j

i

j

i

i

j

ij

i

j

j

j

i

i

i

i

j

j

i

j

i

j

j

IJ

j

I

i

I

j

IJ

J

J

のように結合することができる．結局限定集計した式は次のようになる．

0

cij dIJ

+ 1;1
ci

j

dI

;1 + c2+1 d +1 + c3 ;1 d ;1 + c4 +1 d +1
J

i

j

I

J

I J

ij

I J

ij

+ 5;1 ;1 ;1 ;1 + 6+1 ;1 +1 ;1 + 7;1 +1 ;1 +1 + 8+1 +1+( 0 + 1 ;1 + 2 +1 + 3 ;1 + 4 +1
+ 5 ;1 ;1 + 6 +1 ;1 + 7 ;1 +1 + 8 +1 +1 ) 2 = ;
ci

dI

j

cij dijcij dici

J

cij di

j

dI

ci

J

cij dijjcij dicij dijjcij dij

dI

j

J

ci

j

dI

+1J +1

cij dijcij di

j

j

=

rIJ

この式の括弧の中は，h 格子の修正値の差分方程対式 (3.21) から;rijになり，この項の大きさは一様残差値に対して rIJ の 1/8 になる．今この部分を無視すれば，粗い格子の修正値を求める式は保存形で次のよう表わされる．

2
3
4;1 (d ;1 ; d )+ c +1 (d +1 ; d )+ c ;1 (d ;1 ; d )+ c +1 (d +1 ; d )+c5;1 ;1 (d ;1 ;1 ; d )+ c6+1 ;1 (d +1 ;1 ; d )+c7;1 +1 (d ;1 +1 ; d )+ c8+1 +1 (d +1 +1 ; d ) = ;r

1

ci

j

I

J

IJ

i

j

I

i

j

I

I

J

J

IJ

i

j

I

J

IJ

J

IJ

i

j

I

J

IJ

i

j

IJ

I J

ij

IJ

ij

I J

(3.22)IJ

IJ

この式に解釈を与えれば，式 (3.18) の 2 階微分の項に関し 1 階差分を，式 (3.19)(3.21) では

f 2 ( +1 ; )gcij ui

j

uij

=h

のように差分近似したのに対し，この式では fc2i+1j (uI +1J ; uIJ )g=H のように差分近似しているということである．限定集計はこの式を用い次々に行われる．また延長補間は式 (3.8) を用いて行われる．

ここで準線形微分方程式 (3.19) の１つの混合境界値問題を考える．その計算領域と計算格子は前節と同じであるが，左側境界に Dirichlet 条件，右側境界に Neumann 条件，上下境界に周期条件が課されるもの
2 なめらかな構造格子に対し写像空間の正方形格子上で計算する場合にはこのような近似は概ね妥当である．

連立 1 次方程式の解法

19とする．ただしここでは微分方程式 (3.18) の項 buxy は上の説明とは異なり d の中に含めるものとする．このとき差分方程式 (3.19) の係数は c0  c1  c2  c3  c4 のみとなり，その連立 1 次方程式は優対角の条件を常に満足する．右側境界では，Neumann 条件 (ux )if j = gj を式 (3.19) に考慮すれば ，差分方程式は( 1f + 2f ) f ;1 + 3f
ci

ci

j

j

ui

j

ci

j

uif j;1 + c4f u f +1 + c f u f = f f
i

j

i j

i

のようになる．あらかじめメインプログラムで fif j に

j

i

j

;2 2f
hci

i

j

gj

j

; 2 2fhci

jgj

を加えておくことが必要である．またこの境界上の限定集計の式は，長波長の残差を想定し rif +1j ;1 = rif j ;1  rif +1j = rif j  rif +1j +1 = rif j +1と置けば ，
rIf J= 1 5 f +( f ;1 +1 5 f ;1 +1 5 f +1 ) 2+( f ;1 ;1 + f ;1 +1 ) 4
: ri

j

ri

j

: ri j: ri j

=

ri

j

ri

j

=

のようになる．上下境界では，周期条件の課される上下の境界点を共に計算している．境界条件に関する文にはコメントを付けてある．
次のプログラムは，上記のサブルーチン MGM2 にならって作られたものである．c,f は差分方程式の係数

0 1
c  c  : : : と右辺 f の配列，a は最も粗い格子の修正値を Gauss 消去法で解くときに用いられる配列である．なお Gauss 消去法による連立 1 次方程式の計算は 1.3 節に述べたサブルーチン GAUSS を引用して行わ

れ，wa はこのときに必要な作業用配列である．MGM2 で用いたサブルーチン PRORLX は，ここでは延長補間のサブルーチン PROLON と Gauss 反復法による緩和計算の部分に分けられる．

!*******************************************************************************! Sol of mixed bound value problem of elliptic pde in rectangul domain by MGM!*******************************************************************************SUBROUTINE MGM3(u,c,f,i1,j1,phi,h,na,naf,rmean,e)
REAL u(0:i1,0:j1),c(0:i1,0:j1,0:4),f(0:i1,0:j1),a(24,25),wa(25), &r1(0:40,0:24),r2(0:20,0:12),r3(0:10,0:6),r4(0:5,0:3), &d1(0:40,0:24),d2(0:20,0:12),d3(0:10,0:6),d4(0:5,0:3)
i2=i1/2 i3=i2/2 i4=i3/2
j2=j1/2 j3=j2/2 j4=j3/2
n=(i4+1)*(j4+1) n1=n+1
na=0 100 na=na+1! Determine residuals r1 and restrictions
DO i=1,i1
ip1=i+1 IF(i==i1)ip1=i-1!Neumann cond
DO j=0,j1 jm1=j-1 IF(j==0)jm1=j1-1!periodic cond
jp1=j+1 IF(j==j1)jp1=1!periodic condr1(i,j) = c(i,j,1)*u(i-1,j)+c(i,j,2)*u(ip1,j) &!comp residuals r1+c(i,j,3)*u(i,jm1)+c(i,j,4)*u(i,jp1) &+c(i,j,0)*u(i,j)-f(i,j)
ENDDO ENDDOCALL RESTRC(r1,i1,j1,r2,i2,j2)!restriction(r1 → r2)
CALL RESTRC(r2,i2,j2,r3,i3,j3)!restriction(r2 → r3)
CALL RESTRC(r3,i3,j3,r4,i4,j4)!restriction(r3 → r4)! Determine corrections d4 by Gaussian elimination
l=j4+1
DO i=0,i4
i0=8*i
cycle_1: DO j=0,J4 j0=8*j k=i*l+j+1
FORALL(m=1:n)a(k,m)=0. a(k,k)=1. a(k,n1)=-r4(i,j)!setup rhs
IF(i==0)CYCLE cycle_1!Dirichlet cond
kml=k-l kpl=k+l im4=i0-4 ip4=i0+4
km1=k-1 kp1=k+1 jm4=j0-4 jp4=j0+4
IF(i==i4)THEN kpl=k-l
ip4=i1ENDIF!Neumann cond
IF(j==0 )THEN km1=k+l-2 jm4=j1-4 ENDIF!periodic cond
IF(j==j4)THEN kp1=k-l+2 jp4=4
ENDIF!periodic cond

20

反

復

法（ 多重格子法）

a(k,kml) = c(im4,j0,1)!setup coef matrixa(k,km1) = c(i0,jm4,3)a(k,k ) =-c(im4,j0,1)-c(ip4,j0,2)-c(i0,jm4,3)-c(i0,jp4,4)a(k,kp1) = c(i0,jp4,4)a(k,kpl) = c(ip4,j0,2)+a(k,kpl)
ENDDO cycle_1 ENDDOCALL GAUSS(a,wa,n,n1,det,1.E-5)!solve linear equations
FORALL(i=0:i4,j=0:j4)d4(i,j) = a(i*l+j+1,n1)!determine d4! Prolongation, determine residuals and corrections by point relaxation
CALL PROLON(d4,i4,j4,d3,i3,j3)!prolongation(d4 → d3)
DO i=1,i3 i0=4*i DO j=0,j3 j0=4*j
im1=i-1 ip1=i+1 jm1=j-1 jp1=j+1
im2=i0-2 ip2=i0+2 jm2=j0-2 jp2=j0+2
IF(i==i3)THEN ip1=i3-1 ip2=i1
ENDIF!Neumann cond
IF(j== 0)THEN jm1=j3-1 jm2=j1-2 ENDIF!periodic cond
IF(j==j3)THEN jp1=1
jp2=2
ENDIF!periodic condr3(i,j) = r3(i,j)+c(im2,j0,1)*(d3(im1,j)-d3(i,j)) &!determine r3+c(ip2,j0,2)*(d3(ip1,j)-d3(i,j)) &+c(i0,jm2,3)*(d3(i,jm1)-d3(i,j)) &+c(i0,jp2,4)*(d3(i,jp1)-d3(i,j))
c0 = -c(im2,j0,1)-c(ip2,j0,2)-c(i0,jm2,3)-c(i0,jp2,4)r3(i,j) = -r3(i,j)/c0!d3 by point relaxation
ENDDO ENDDOFORALL(i=1:i3,j=0:j3)d3(i,j)=d3(i,j)+r3(i,j)!determine d3
CALL PROLON(d3,i3,j3,d2,i2,j2)!prolongation(d3 → d2)
DO i=1,i2 i0=2*i DO j=0,j2 j0=2*j
im1=i-1 ip1=i+1 im0=i0-1 ip0=i0+1
jm1=j-1 jp1=j+1 jm0=j0-1 jp0=j0+1
IF(i==i2)THEN ip1=i2-1 ip0=i1
ENDIF!Neumann cond
IF(j==0 )THEN jm1=j2-1 jm0=j1-1 ENDIF!periodic cond
IF(j==j2)THEN jp1=1
jp0=1
ENDIF!periodic condr2(i,j) = r2(i,j)+c(im0,j0,1)*(d2(im1,j)-d2(i,j)) &!determine r2+c(ip0,j0,2)*(d2(ip1,j)-d2(i,j)) &+c(i0,jm0,3)*(d2(i,jm1)-d2(i,j)) &+c(i0,jp0,4)*(d2(i,jp1)-d2(i,j))
c0 = -c(im0,j0,1)-c(ip0,j0,2)-c(j0,jm0,3)-c(i0,jp0,4)r2(i,j) = -r2(i,j)/c0!d2 by point relaxation
ENDDO ENDDOFORALL(i=1:i2,j=0:j2)d2(i,j)=d2(i,j)+r2(i,j)!determine d2
CALL PROLON(d2,i2,j2,d1,i1,j1)!prolongation(d2 → d1)
DO i=1,i1 DO j=0,j1
im1=i-1 ip1=i+1jm1=j-1 jp1=j+1
IF(i==i1)ip1=i-1!Neumann cond
IF(j==j1)jp1=1!periodic cond
IF(j==0 )jm1=j1-1!periodic condr1(i,j) = r1(i,j)+c(i,j,0)*d1(i,j)&!determine r1+c(i,j,1)*d1(im1,j)+c(i,j,2)*d1(ip1,j) &+c(i,j,3)*d1(i,jm1)+c(i,j,4)*d1(i,jp1)
ENDDO ENDDOFORALL(i=1:i1,j=0:j1)d1(i,j) = d1(i,j)-r1(i,j)/c(i,j,0)!determine d1! Apply corrections
DO i=1,i1
FORALL(j=0:j1)u(i,j)=u(i,j)+d1(i,j)!apply d1
um=u(i,0)+u(i,j1) u(i,0)=(um-phi)/2. u(i,j1)=(um+phi)/2. !periodic cond
ENDDO
! Decide convergencermean=0.

連立 1 次方程式の解法

21

DO i=1,i1 DO j=0,j1
rmean = rmean+ABS(r1(i,j))/(i1*j1)!comp mean residual
ENDDO ENDDO
WRITE(20,'(I3,F12.6)')na,rmean!output converging processIF(rmean>e.AND.na<naf)
GOTO 100!decide convergence
END!***** Subroutine RESTRCSUBROUTINE RESTRC(r1,i1,j1,r2,i2,j2)
REAL r1(0:i1,0:j1),r2(0:i2,0:j2)
DO jj=0,j2 j=2*jj jm1=j-1 jp1=j+1
IF(j==0)jm1=j1-1 IF(j==j1)jp1=1!periodic cond
DO ii=1,i2-1 i=2*ii
r2(ii,jj)=r1(i,j)+(r1(i-1,j)+r1(i+1,j)+r1(i,jm1)+r1(i,jp1))/2. &+(r1(i-1,jm1)+r1(i+1,jp1)+r1(i+1,jm1)+r1(i-1,jp1))/4.
ENDDO
r2(i2,jj)=1.5*r1(i1,j)+(r1(i1-1,j)+1.5*r1(i1,jm1) &!Neumann cond
+1.5*r1(i1,jp1))/2.+(r1(i1-1,jm1)+r1(i1-1,jp1))/4.
ENDDO
END

このサブルーチンについて，!Determine corrections d4 by Gaussian elimination のブロックで akkpl =cip4j +akkpl について，右辺に akkpl を加えることはその上で akm = 0 (m = 1 : : : ) と置いているので一2

般には意味のないことであるが，Neumann 条件の課される右側境界では必要なことに注意．またこの文の前に akkml = c1im4j を置かなければならない．!Apply corretions のブロックで周期条件を確実に満足させる操作は必要なもので，省略すると解が発散またはほとんど 収束しても最後に発散することになる．副サ
ブルーチン PROLON は前項のものと同じである．

最後に，サブルーチン MGM3 を引用する計算例として，極座標系 r  における Poisson 方程式の混合境界値問題

urr + 1r ur + r12 u = 1:0(1 < r < 3 0 <  < =4)u(1 ) = 10+  =(=4) ur (3 ) = ;0:5u(r =4) = u(r 0)+  u (r =4) = u (r 0)を解くことにする．ただし  =(3.23a)(3.23b);1:0．基本格子は r 方向に 40 等分， 方向に 24 等分したものとする．

式 (3.23) の差分方程式は，前の説明とは異なり左辺に 1 階微分の項も含めることにすれば次のようになる．

(1 ; gij )ui;1 j +(1+ gij )ui+1 j + hij ui j;1 + hij ui j+1 ; 2(1+ hij )uij = r2
ただし gij(3.24)

 ( r)=2rij  hij  ( r=rij )2  r = 1=20  = =4=24 である．この差分方程式の係数

1  gij は本来 格子間隔によるが，ここでは粗い格子で修正値を求める式にもそのまま使うことにする．この混合境界値問題のメインプログラムは次のようになる．
program MAIN
REAL u(0:40,0:24),c(0:40,0:24,0:4),f(0:40,0:24)
DATA i1,j1,grad,phi,naf,e/40,24,-.5,-1.,100,1.E-5/
dr=2./FLOAT(i1) dt=3.14159/4./FLOAT(j1)
DO I=1,I1 DO J=0,J1
r=1.+dr*i g=dr/2./r h=dr*dr/(r*dt*r*dt)c(i,j,1)=1.-g c(i,j,2)=1.+g c(i,j,3)=h c(i,j,4)=h c(i,j,0)=-2.*(1.+h)f(i,j)=1.*dr*drENDDO ENDDO

22

反

復

法（ 多重格子法）

FORALL(i=1:i1)f(i, 0)=f(i, 0)+c(i, 0,3)*phi
FORALL(i=1:i1)f(i,j1)=f(i,j1)-c(i,j1,4)*phi
DO j=0,j1
u(0,j)=10.+j/FLOAT(j1)*phi
f(i1,j)=f(i1,j)-grad*2*dr*c(i1,j,2)
FORALL(i=1:i1)u(i,j)=u(0,j)
ENDDO
OPEN(20,FILE='OUTPUT.dat')
CALL MGM3(u,c,f,i1,j1,phi,h,na,naf,rmean,e)
DO j=j1,0,-1
WRITE(20,'(41F8.3)')(u(i,j),i=0,i1)
ENDDO
CLOSE(20)
STOP END

計算の結果は，出発値 u!periodic cond

!Dirichlet cond
!Neumann cond!starting values

 0 を与えた場合にはサイクル数 na = 64 で収束し ，また Dirichlet 条件と同じ周期条件を満足する出発値 u(r ) = 10+  =(=4) を全領域に与えた場合には na = 8 で収束した．

3.4

多重格子法の各種問題への拡張

多重格子法は，差分法と同様に，有限要素法 (FEM, ntie element method) の計算効率をあげるために用いられている．通常 四辺形双 1 次要素または四辺形ラグランジュ要素が用いられ，有限要素方程式 (節点値 ui の連立 1 次方程式) が多重格子法で解かれる．その限定集計と延長補間が差分法と異なるところで，差分法のものを大らかに利用する方法と FEM らし く行う方法がある．前者の方法でも，最も細かい基本
メッシュに対しては，もちろん有限要素方程式から残差を正確に計算しすべての残差が 0 になるまで反復計算を履行するのであるが，その上に重なる粗いメッシュに対しては大らかに計算するということである．
その限定集計と延長補間は下図に示すように大胆に行われる．メッシュのゆがみや間隔の変化が大きいとこ
ろでは，限定集計と延長補間に伴う誤差が大きくなり，サイクル数は増えると思われるが，結果の精度には無関係である．






((s s5 (1
( ((s6 es 2 s1  r h = rh + 21 (rh +rh +rh +rh)
+ 4 (rh + rh + rh + rh )(((s3 s 0s8 (((s7   4  
4rh





( ((
限定集計
e4




((e310sQ (( ((1se 6 sQ QQe 2( ((QQ QQQQdhd=h (=d (hd+hd+hd)=h2+d h +d h)=4
QQQ h h(((





d =d
( ((
2
0

0

1

2

3

4

5

6

7

8

0

10

6

1

2
1

延長補間

2
1

2
1

2
2

2
2

2
3

2
4

連立 1 次方程式の解法

23

後者の方法は，メッシュのゆがみや間隔の変化が大きいところ，また四辺形や六面体以外の要素にもそのまま適用できる．形状関数を Ni とすれば ，限定集計の式は

ri2h = Ih2h rh ]i =Z

Ni2h (x )Njh (x )rjh dxx

e2h(3.25)

また延長補間の式は

dhj = I2hh d2h ]j = Ni2h (x j )d2i h(3.26)となる．ここで蛇足を加えれば ，FEM では要素の節点 x i に節点値 ui (ri  di  : : : ) が定義され，要素内の任意の点 x の u は補間関数 u(x) = Ni (x )ui から求められる．u の微分値，積分値は補間関数を微分または積

分することによって求められる．ある節点の微分値はその節点を共有する要素のものが一般に異なる値を示すことになる．このように異なる値を示すものは形状関数 Ni の重みを付けて積分し平均したものが取られる．延長補間の式 (3.26) は，細かいメッシュh の節点 x j における補間値 dhj は，その点が粗いメッシュ

2h の n 番めの要素内にあれば，その要素の補間関数 Ni2h (x )d2i h から求められるということである．また式(3.25) の Njh (x )rjh は細かいメッシュh の残差の補間値で，この限定集計の式は，粗いメッシュ2h の残差 ri2hがこの補間値に形状関数 Ni2h (x ) の重みを付けて積分したものになることを意味している．粗いメッシュのすべて節点にわたってこの残差の積分値を合計したものは，任意の点の形状関数の合計が 1 であることから，細かい格子の残差値を計算領域にわたって積分したものになる．このことからも，式 (3.25) は限定集計の式として妥当なものであることが分かる．
多重格子法の双曲型方程式への適用例として，Jameson(1979) の遷音速ポテンシャル流れ方程式の解法の概要を紹介する3 ．定常２次元遷音速流れの完全ポテンシャル方程式は，大部分の亜音速流れ域では楕円型であるが一部の超音速流れ域では双曲型になる．超音速域ではある点の状態はその点から上流側に拡がる
円錐状の依存域 (domain of dependence) 内の状態にのみ依存し ，逆にその点の状態はその点から下流側に広がる円錐状の影響域 (domain of inuence) 内の状態にのみ影響を及ぼす．なお亜音速域では依存域，影
響域は流れ場全域となる．離散方程式はこの物理現象を反映したものであるべきで，超音速域で中心差分を上流差分に切替える方法が用いられている．また超音速域に人工粘性を付加する方法も用いられているが，上流差分は中心差分に拡散項を付加したものであることからこれら２つの方法は本質的に同じものである．
この流れの微分方程式は線形化の後L  axx + byy = 0(3.27)となる．この完全ポテンシャル方程式 (full potential equation) は， 形陰解法の式( ; A x2 )( ; B y2 )  = ! L

(3.28)に書換え，ADI 法を適用し解かれるている．ただし A = a=x2  B = b= t2 ， x y は中心差分の演算子，
 4 min(A B )，!  1:5，不等間隔の曲線座標格子を使用．MGM では波長 2 x の残差の波を減衰させれば十分で，それにはをこのように選ぶのが適当である．この方法は亜音速流れには有効であるが，流れが超音速になり，すなわち A < 0，式 (3.27) が双曲型になるところでは不安定になり，このままでは遷音速流れには適用できない．

Jameson は，遷音速ポテンシャル流れの方程式として，式 (3.27) の代わりに双曲形偏微分方程式0

t + 1 xt + 2 yt = axx + byy

(3.29)

3 Jameson, A., Acceleration of trnsonic potential ow calculations on arbitrary meshes by the Multiple grid method,
AIAA 4th CFD Conf., (1979-7), 122{146.

24

反

復

法（ 多重格子法）

を用いる解法を提案している．この式は解が収束すれば @=@t = 0 になるので式 (3.27) と等価になる．式 (3.29)
を ADI 法で解くための  形陰解法の式は次のようになる．

(S ; A x2 )(S ; B y2 )  = !S L
ただし S(3.30)

 0 + 1 x; + 2 y; ， 0  1  2 は 0  1  2 による定数， x;  y; は上流差分の演算子である．

A < 0 のところでは x2 も上流差分のものに切替えることが必要である．MGM の計算手順は，最細格子→

最粗格子→→二番目に細かい格子を 1 サイクルとし ，各格子レベルで式 (3.30) を ADI 法で 1 回計算する．

1 サイクルの計算量は最細格子の計算量の約 1 32 になる．遷音速翼まわりの流れの計算例では，192  32 最細格子，5 レベルの多重格子の場合に 10 サイクルで収束解が得られる．
多重格子法は差分法と FEM では離散方程式の効率的な解法としてかなり用いられているが，境界積分法

ではあまりない．最後に多重格子法のパネル法への適用について述べた Sloo(1982) の解説4を参考までに引用する．Atkinson(1973) は，次の第 2 種 Fredholm 型積分方程式の解法に多重格子法を用いている．
Z 1

f (x) ; K (x  )f ( ) d = g(x)(3.31)0

この方程式の計算式は

f (xi ) ;Xjとなるが，行列 I

Wj K (xi  xj )fj hj = g(xi )(3.32); Wj Kij hj ] = A，ベクトル ff (xj )g = x fg(xi )g = b と置けば連立 1 次方程式

Ax = b

(3.33)になる．1 サイクルの計算は次のようなものである．細かい格子 G2 上で，A2 x2 = b2 を Jacobi 緩和法で解いて x2 を求める．r2 = A2 x2を直接法で解いて

x1 を求める． x2 = I12 x1 を Nystrom 補間すなわち

f (x) = 1 g(x) +; b2 を計算する．粗い格子 G1 に対し ，r1 = I21 r2 を計算する．A1 x1 = r1

Xj

Wj K (x xj )fj hj(3.34)

から計算する．x2 を修正する． この方法では G1 に十分粗い格子を用いているので直接法が適用できる．また式 (3.34) の計算には平滑化操作が入るので，Jacobi の緩和計算を省略することもできる．第 1 種Fredholm 積分方程式は，特異核関数の特異性を分離し ，第 2 種の形に帰着させて解かれる．揚力の作用す
る平板翼まわりの流れに関する積分方程式を，24 パネルを用いて解く場合に，多重格子法では 6 サイクルで収束する．なおこの問題は Jacobi 緩和法のみでは発散する．
er 24 panels;e? er er r; 12 panels; e;;e?;;e; ?
?
6 panels
1 MG cycle -e determine residual? restriction; prolongation(Nystrom) determine correction
r apply correction

4 Sloo, J.W., Requirements and developments shaping a next generation of integral methods, Numerial Methods in
Aeronautical Fluid Dynamics, Roe, P.L., ed., p.517, 1982, Academic Press.