放物型偏微分方程式の数値解法２ [ネコ騙し数学]
放物型偏微分方程式の数値解法２
前回に引き続き、差分法を用いた熱伝導方程式の数値解法について述べることにする。
この方程式の時間微分を
右辺の空間微分を
と近似すると、（１）式は
となる。
とおくと、
という３重対角行列を係数を持つ連立方程式が得られる。
（１）の解をT=T(x,t)とし、初期条件を
境界条件を
とする。
そして、
とすると、初期条件より
境界条件より
となる。
したがって、
という連立方程式が得られ、逐次的に、３重対角行列アルゴリズム（TDMA）を用いて解くことによって、微分方程式（１）の近似解を求めることができる。
このような解法を陰解法という。
陰解法は、前回紹介した陽解法とは異なり、無条件安定である。
陽解法との比較参照のために、前回と同じ問題
問題
を、初期条件
境界条件
のもとで、Δx=1、Δt=1として解け。
を解かせたものは次の通り。
この問題の場合、厳密解と比較すると、陰解法による近似解は厳密解よりも大きいことがわかる。
さらに、前回の陽解法による結果を加え比較したものが次のグラフである。
込み入っていて、少しわかりづらいと思うのだけれど、厳密解は陰解法と陽解法の間に位置している。
ということで、陰解法と陽解法によって得られた数値解を足して２で割ると厳密解に近い結果が得られることが予想できる。
この予想は正しく、陽解法と陰解法の中間的な解法、クランク・ニコルソン法という解法が存在する。
Δt=1、Δx=1という非常に粗い計算でありながら、クランク・ニコルソン法は、厳密解とよく一致していることが分かる。
ただし、この事実をもって、陽解法や陰解法は精度が悪いと思わないで欲しい。
Δt=0.1程度にとれば、陽解法、陰解法ともに良好な結果を得ることができる。
陽解法はともかく、陰解法は熱や流れなどの数値計算ではよく使われている、優れた解法である。
問題の初期条件を
境界条件を
と変え、陰解法を用いてΔt=0.5、Δx=0.4として計算した結果を以下に示す。
このようにあらい格子でも、厳密解をよく再現していることが分かる。
そして、t→∞の定常解
に収束することが確認できる。
誤差は最大で１０％くらい。
Δt=0.1、Δx=0.4にすると、誤差は最大で約３％。
陰解法よりより高精度な計算ができる、クランク・ニコルソン法は、次回ということで。
計算に使用したFortranプログラムは次の通り。
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)+x/l
end do
! boundary condition
do j=1, nt
t(0,j)=0.
t(nx,j)=1.
end do
do j=1, nt
do i=1,nx-1
a(i)=-kappa/dx2
b(i)=2.*kappa/dx2 + 1/dt
c(i)=-kappa/dx2
d(i)=t(i,j-1)/dt
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(*,*) ' *** result *** '
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
write(*,*)
write(*,*) 'error(percent)'
do j=1,nt
write(*,100) j*dt, ((t(i,j) - exact(i*dx,j*dt))/exact(i*dx,j*dt)*100, i=1,nx-1)
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+x/4.
end
