SSブログ

楕円形方程式(ラプラス方程式)の数値解法 [数値解析]

楕円形方程式(ラプラス方程式)の数値解法

 

ラプラス方程式

  

の差分法を用いた数値解法を考える。

  

だから、(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

 


nice!(0)  コメント(0) 

nice! 0

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

この広告は前回の更新から一定期間経過したブログに表示されています。更新すると自動で解除されます。