楕円形方程式(ラプラス方程式)の数値解法 [数値解析]
楕円形方程式(ラプラス方程式)の数値解法
ラプラス方程式
の差分法を用いた数値解法を考える。
だから、(1)は
と近似できる。
(3)式は優対角条件
を満たしており、連立方程式(3)を反復法(ガウス・ザイデル法)を用いて解くことができる。
特にΔx=Δyのとき(2)式は
問
を、Δx=Δy=0.1として、(4)を用いて数値的に解け。
計算結果
計算に使用したプログラム
parameter (nx=10,ny=10)
real u(0:nx,0:ny)
u=0.0
limit = 1000
pi=acos(-1.0)
dx=1./nx
dy=1./ny
! 境界条件 (y=1)
do i=0,nx
u(i,ny)=sin(pi*i*dx)
end do
do iter=1, limit
res_max=0.
do i = 1, nx-1
do j=1, ny-1
res=u(i,j)-0.25*(u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1)) ! 残差
res_max=amax1(abs(res),res_max)
u(i,j)= 0.25*(u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1))
end do
end do
if (res_max.lt.1.e-6) exit ! 収束判定
end do
do j=0,ny
write(*,100) (u(i,j),i=0,nx)
end do
100 format (10(f8.5,1x),f8.5)
end
Δx≠Δyの場合でも計算できる、より汎用的なプログラムは次の通り。なお、このプログラムでは、反復計算の収束を早めるためにSOR法を用いている。
parameter (nx=10,ny=10)
real u(0:nx,0:ny)
real aw(nx,ny),ae(nx,ny),as(nx,ny),an(nx,ny),ap(nx,ny)
u=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=1000
omega=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
res_max=0.0
do i=1,nx-1
do j=1,ny-1
sum = aw(i,j)*u(i-1,j)+ae(i,j)*u(i+1,j)+as(i,j)*u(i,j-1)+an(i,j)*u(i,j+1)
res = sum-ap(i,j)*u(i,j)
res_max=amax1(res_max,abs(res)/ap(i,j))
u(i,j)=u(i,j)+omega*res/ap(i,j)
end do
end do
if (res_max.lt.1.e-6) exit
end do
do j=0,ny,1
write(*,100) (u(i,j),i=0,nx,1)
end do
100 format (10(f8.5,1x),f8.5)
end
コメント 0