ラプラス方程式の数値解法２

ラプラス方程式の数値解法２

 

ラプラス方程式

　　

に差分法を適用すると、

　　

という連立差分方程式が得られる。

前回、この連立方程式をガウス・ザイデル法とSOR法といった反復法を用いて解いた。しかし、ガウス・ザイデル法、SOR法ともに、１回の反復で計算の影響が隣接する格子にしか及ばないので、収束が遅く、収束解を得るまでに多くの反復計算を要する。そこで、この反復計算を減らす方法について、今回、考えることにする。

その前に、Patankharの流儀にならって（３）式を

　　

と書き換えることにする。

ここで、

　　

であり、

　　

である。

WはWest、EはEast、SはSouth、NはNorthを省略したものくらいの意味。

 

ガウス・ザイデル反復法は、計算領域の各点に対して

　　

と計算する方法である。このため、反復計算ごとの、の更新値には、これに隣接するの４点の値しか影響を与えない。このため、ガウス・ザイデル法の収束は遅くなる。

そこで、

　　

と変形する。すると、左辺は３項係数行列の方程式形となり、

　　

とすることによって、TDMA（トマスのアルゴリズム）を用いてj列の（の推測値）を一気に計算することができる。

j列のすべてについてこのように計算することによって、１回の反復計算の収束速度を速め、反復回数を減らすことができる。

 

一方、SOR法は計算領域の各点に対して

　　

と計算する。

ここで、ωは収束を早めるための加速係数で1<ωであり、は反復前のの値である。

（８）式は

　　

と書き換えることができる。

したがって、（７）式の係数を

　　

とすることによって、さらに収束を早めることができる。

 

 

parameter (nx=100,ny=100)
real u(0:nx,0:ny)
real aw(nx,ny),ae(nx,ny),as(nx,ny),an(nx,ny),ap(nx,ny)
real a(nx),b(nx),c(nx),d(nx)

u=0.
a=0.; b=0.; c=0.; d=0;

aw=0. ; ae=0. ; as=0. ; an=0. ; ap=0.

pi=acos(-1.0)
dx=1./nx; dx2=dx*dx
dy=1./ny; dy2=dy*dy

limit=400
alp=1.3

! 境界条件 (y=1)
do i=0,nx
    u(i,ny)=sin(pi*i*dx)
end do

! 係数の計算
do i=1,nx-1
    do j=1,ny-1
        aw(i,j)=1./dx2
        ae(i,j)=1./dx2
        as(i,j)=1./dy2
        an(i,j)=1./dy2
        ap(i,j)=aw(i,j)+ae(i,j)+as(i,j)+an(i,j)
    end do
end do

do iter=1, limit
    e_max=0.0
    do j=1,ny-1
        do i=1,nx-1
            a(i)=-aw(i,j)
            b(i)=ap(i,j)/alp
            c(i)=-ae(i,j)
            d(i)=as(i,j)*u(i,j-1)+an(i,j)*u(i,j+1)+(1.-alp)*ap(i,j)/alp*u(i,j)
        end do
        d(1)=d(1)+aw(1,j)*u(0,j)
        d(nx-1)=d(nx-1)+ae(nx-1,j)*u(nx,j)
       
        call tdma(a,b,c,d,nx-1)
       
        do i=1,nx-1
            cor=u(i,j)-d(i)
            e_max=amax1(e_max,abs(cor))
            u(i,j)=d(i)
        end do
    end do
   
    if (e_max.lt.1.e-6) exit
end do

write(*,*) iter,e_max

do j=0,ny,10
    write(*,100) (u(i,j),i=0,nx,10)
end do

open(1,file='Laplace.dat', status = 'replace')
do i=0,nx
    do j=0,ny
        write(1,110) i*dx, j*dy, u(i,j)
    end do
    write(1,*)
end do
close(1)

100 format (10(f8.5,1x),f8.5)
110 format(2(f8.5,1x),f8.5)

end

!  TDMA
subroutine tdma(a,b,c,d,n)
real a(n), b(n), c(n),d(n)

! 前進消去
do i=1,n-1
    ratio=a(i+1)/b(i)
    b(i+1)=b(i+1)-ratio*c(i)
    d(i+1)=d(i+1)-ratio*d(i)
end do

! 後退代入
d(n)=d(n)/b(n)
do i=n-1,1,-1
    d(i)=(d(i)-c(i)*d(i+1))/b(i)
end do

end

 

 

　　

を、x、y方向ともに１００分割し、Δx=Δy=0.01とし、ω=1.3、さらに収束条件をとした場合、SOR法の反復回数は３３７６回であるのに対し、直接法TDMAとSOR法を組み合わせた計算法の反復回数は３５９回と反復回数を約1/10に抑えることができる。

 

PCの計算速度がべらぼうに速くなった現在でも、これくらいの計算量になると、計算速度の差を十分に体感できる！！

ただし、計算格子の数が５０×５０以下ならば、反復回数はSOR法の方が多いけれど、サブルーティンを使っていないのでむしろSOR法の方が計算時間は短いのではないだろうか(^^ゞ

 

今回紹介した計算法よりもさらに収束を早めるADI法が存在するけれど、基本的な考え方は今回紹介した方法と同様であり、プログラムの複雑さの割にたいして速くなるわけでもないので、割愛するにゃ。