So-net無料ブログ作成

放物型偏微分方程式の解法4 まとめ [ネコ騙し数学]

放物型偏微分方程式の解法4 まとめ

 

熱伝導方程式

  

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

  

(純)陰解法を用いると、

  

となる。

そして、クランク・ニコルソン法の場合は、(2)、(3)式の辺々を足しあわせ2で割った

  

 

である。

ここで、

  

 

また、(2)、(3)、(4)式は

と一つの式表すことができ、

  

になる。

(6)式を

  

と変形し、

  

とおくと、(7)式は

  

となり、陽解法、陰解法、クランク・ニコルソン法と解法の区別をすることなく、三重対角行列を係数に持つ連立方程式として(1)の近似解を求めることができる。

 

この方針にしたがって作ったプログラムは以下の通り。

 

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

 

陽解法の時にも連立方程式を解いているので、陽解法のとき、少し無駄な計算をすることになるけれど、陽解法、陰解法、クランク・ニコルソン法といった解法の区別をすることなく計算できる。

また、αを一種の重みと考え、α01/21だけではなく、0≦α≦1の任意のαを選び計算することもできる。

ただ、

にとると、

  

の値によって、伝播誤差のために、振動解が得られる。

 

例によって、(1)のκ=1/2とし、条件

  

とし、Δt=0.1Δx=0.4に選んで計算した結果を以下に示す。

陽解法
Explicit.png

陰解法

Implicit.png

クランク・ニコルソン法

Crank-Nicolson method.png

 

ΔtΔxをこの程度にとれば、陽解法、陰解法、クランク・ニコルソン法ともに、(1)の厳密解とよく一致することが分かる。

 

また、Δt=0.2Δx=0.4とし、陽解法で解いた結果は次の通り。

 

Explicit2.png

 

このとき、

  

となるので、この条件で陽解法を用いて解こうとすると、伝播誤差のために、計算を進めるたびに、数値解の振動が激しくなり、厳密解からの逸脱が大きくなることが分かる。

 


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

nice! 0

コメント 0

コメントを書く

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