放物型偏微分方程式の解法４ まとめ [ネコ騙し数学]
放物型偏微分方程式の解法４ まとめ
熱伝導方程式
を、陽解法を用いて離散化すると、
（純）陰解法を用いると、
となる。
そして、クランク・ニコルソン法の場合は、（２）、（３）式の辺々を足しあわせ２で割った
である。
ここで、
また、（２）、（３）、（４）式は
と一つの式表すことができ、
になる。
（６）式を
と変形し、
とおくと、（７）式は
となり、陽解法、陰解法、クランク・ニコルソン法と解法の区別をすることなく、三重対角行列を係数に持つ連立方程式として（１）の近似解を求めることができる。
この方針にしたがって作ったプログラムは以下の通り。
parameter (nx_max=20, nt_max=100)
real t(0:nx_max,0:nt_max)
real kappa,l
nt = 10; nx = 8
l=4.0
dt = 0.25
dx = l/nx
dx2 = dx*dx
kappa=0.5
write(6,*) '解法を入力：陽解法 1, 陰解法 2, クランク・ニコルソン 3'
read(5,*) iflag
if(iflag.eq.1) then
alpha=0.
else if (iflag.eq.2) then
alpha = 1.
else
alpha=0.5
end if
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
call solver(t,nx,nt,r,alpha)
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
!close(1)
100 format(12(f8.5, 1x))
end
! 厳密解
function exact(x,t)
pi = acos(-1.0)
s=0.
do i=1,30
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
! 【注意】 ここから下はいじってはいけない
! 差分法を用いて解く本体
subroutine solver(t,nx,nt,r,alpha)
real t(0:20,0:100)
real a(nx), b(nx), c(nx), d(nx)
do j=1, nt
do i=1,nx-1
a(i)=-alpha*r
b(i)=1+2.*alpha*r
c(i)=-alpha*r
d(i)=t(i,j-1)+(1-alpha)*r*(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
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
陽解法の時にも連立方程式を解いているので、陽解法のとき、少し無駄な計算をすることになるけれど、陽解法、陰解法、クランク・ニコルソン法といった解法の区別をすることなく計算できる。
また、αを一種の重みと考え、αを0、1/2、1だけではなく、0≦α≦1の任意のαを選び計算することもできる。
ただ、
にとると、
の値によって、伝播誤差のために、振動解が得られる。
例によって、（１）のκ=1/2とし、条件
とし、Δt=0.1、Δx=0.4に選んで計算した結果を以下に示す。
陽解法
陰解法
クランク・ニコルソン法
Δt、Δxをこの程度にとれば、陽解法、陰解法、クランク・ニコルソン法ともに、（１）の厳密解とよく一致することが分かる。
また、Δt=0.2、Δx=0.4とし、陽解法で解いた結果は次の通り。
このとき、
となるので、この条件で陽解法を用いて解こうとすると、伝播誤差のために、計算を進めるたびに、数値解の振動が激しくなり、厳密解からの逸脱が大きくなることが分かる。
