放物型偏微分方程式の解法３ クランク・ニコルソン法

放物型偏微分方程式の解法３　クランク・ニコルソン法

 

微分方程式

　　

を、陽解法を用いて離散化すると、

　　

となる。

陰解法を用いて離散化すると、

　　

（２）と（３）の平均を取ると、

　　

になる。

　　

とおくと、

　　

になる。

ここで、

とおくと、

　　 

と３重対角行列を係数に持つ連立方程式が得られ、三重対角行列アルゴリズムTDMAを用いてこの連立方程式を解くことが可能になる。

 

大胆に、（６）式のと以外の項を落とすと、

　　

となり、r>0のとき、

　　

となり、クランク・ニコルソン法は無条件安定であることが予想できる。

 

例によって、次の問題をクランク・ニコルソン法で解いてみる。

 

問題

　　

を、初期条件

　　

境界条件

　　

のもとで、Δt=1、Δx=1として解け。

 

 

 

Δt=1、Δx=1という粗い計算でも、良好な計算結果が得られていることが分かる。

 

Δt=0.5、Δx=0.4にすると、厳密解とクランク・ニコルソン法による近似解はほとんど等しくなる。

 

 

この計算に使用したプログラムは次の通り。

 

! crank-nickolson method
parameter (nt = 20, nx = 10)
real t(0:nx,0:nt)
real a(nx),b(nx),c(nx),d(nx)
real kappa,l

l=4.0
dt = 0.5
dx = l/nx
dx2 = dx*dx
kappa=0.5

do i=0, nx
    do j=0, nt
        t(i,j)=0.
    end do
end do

! initial conditon
do i=0, nx
    x=i*dx
    t(i,0)=x*(l-x)
end do

! boundary condition
do j=1, nt
    t(0,j)=0.
    t(nx,j)=0.
end do

r= kappa*dt/dx2
do j=1, nt
    do i=1,nx-1
        a(i)=-r/2
        b(i)=1+r
        c(i)=-r/2
        d(i)=t(i,j-1)+r/2*(t(i-1,j-1)-2*t(i,j-1)+t(i+1,j-1))
    end do
    d(1)=d(1)-a(1)*t(0,j)
    d(nx-1)=d(nx-1)-c(nx-1)*t(nx,j)
    call tdma(a,b,c,d,nx-1)

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

write(*,*) '*** Crank-Nicolson Method ***'
do j=0,nt
    write(*,100) j*dt, (t(i,j), i=0,nx)
end do

write(*,*)
write(*,*) '*** exact solution ***' 
do j=0,nt
    write(*,100) j*dt, (exact(i*dx,j*dt), i=0,nx)
end do

100 format(12(f8.5, 1x))

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

! 厳密解
function exact(x,t)
pi = acos(-1.0)
s=0.
do i=1,10
    s=s+128/((2*i-1)**3*pi**3)*exp(-(2*i-1)**2*pi**2/32*t)*sin((2*i-1)*pi/4*x)
end do
exact = s
end

 

タグ：数値解析 微分方程式 微分積分