So-net無料ブログ作成

放物型偏微分方程式の解法3 クランク・ニコルソン法 [ネコ騙し数学]

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

 

微分方程式

  

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

  

となる。

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

  

(2)と(3)の平均を取ると、

  

になる。

  

とおくと、

  

になる。

ここで、

とおくと、

   

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

 

大胆に、(6)式の以外の項を落とすと、

  

となり、r>0のとき、

  

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

 

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

 

問題

  

を、初期条件

  

境界条件

  

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

 

 

ht-graph-004.png

 

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

 

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

 

ht-graph-010.png

ht-graph-011.png

ht-graph-012.png

 

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

 

! 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

 


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

nice! 0

コメント 0

コメントを書く

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

この広告は前回の更新から一定期間経過したブログに表示されています。更新すると自動で解除されます。