バーガス方程式の解法２ （純）陰解法

バーガス方程式の解法２　（純）陰解法

 

前回に引き続き、バーガース方程式

　　

の数値解法について考える。

 

（１）式の時間微分に前進差分

　　

空間微分に中心差分を用いると、

　　

（１）式は

　　

したがって、

 　　

と、３重対角係数を持つ連立方程式が得られ、TDMAで解くことができる。

 

ノイマンの安定性判定によれば、前回の陽解法とは異なり、無条件安定である。

 

（１）に

を課し、ν=0.01、c=1.0、Δx=0.01、Δt=0.005の場合について、陰解法で解いた結果は次の通り。

 

 

計算結果を比較するために、t=0.5のときの陽解法と陰解法、さらにクランク・ニコルソン法の計算結果を以下に示す。

陽解法（Explicit）と比較すると、陰解法（Implicit）のほうが曲線の裾野が広がっていることがわかる。

また、クランク・ニコルソン法の計算結果を厳密解と考えると、陽解法は（１）式の右辺の拡散項を過小評価し、対して陰解法は過大評価しており、陽解法と陰解法の誤差はほぼ同程度と考えられる。

 

陽解法、陰解法の打切誤差は

　　

程度なので、まぁ、理論通りの結果と考えることもできる。

なお、クランク・ニコルソン法の打切誤差は

　　

なので、陽解法、陰解法より高精度に計算することができる（はずである）。

 

この計算に用いた陰解法のプログラムを以下に示す。

 

 

! 陰解法でBurgers方程式∂u/∂t+c∂u/∂x=ν∂²u/∂x²を解く
parameter(nx=100, nt=200)
real u(0:nx,0:nt)
real a(nx),b(nx),c(nx),d(nx)

a=0; b=0; c=0; d=0

cnyu=0.01 ! cnyu=ν
dx=1./nx
dx2=dx*dx
cu=1.
dt = 0.005

cc = cu*dt/dx !クーラン数
r=cu*dx/cnyu !格子レイノルズ数

u=0

! initial condition :u(x,0)=1-10x in 0<=x<=0.1
do i=0,10
    u(i,0)=1.0-10.0*(i*dx)
end do

! boundary condtion : u(0,t)=1.0
do j=0, nt
    u(0,j)=1.0
end do

do j=1, nt
    do i=1,nx-1
        a(i)= -cc*(1./r+0.5)
        b(i)= (1+2.*cc/r)
        c(i)=-cc*(1./r-0.5)
        d(i)=u(i,j-1)
    end do
    d(1)=d(1)-a(1)*u(0,j) ! 境界条件（x=0）
    d(nx-1)=d(nx-1)-c(nx-1)*u(nx,j) ! 境界条件（x=1）

    call tdma(a,b,c,d,nx-1)

    do i=1,nx-1
        u(i,j)=d(i)
    end do
end do

! 結果の出力
do j=0,nt,20
    write(6,100) j*dt,(u(i,j),i=0,nx,10)
end do

! お絵かきするためにファイルに書き出す
open (1,file='shine2.dat',status='replace')
do i=0,nx
    write(1,*) i*dx, (u(i,j),j=0,nt,20)
end do
close(1)

100 format(11(f8.5,1x),f8.5)
110 format(11(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