So-net無料ブログ作成

放物型偏微分方程式の解法 ルンゲ・クッタ法 [ネコ騙し数学]

放物型偏微分方程式の解法 ルンゲ・クッタ法

 

次の熱伝導方程式

  

の右辺を差分法を用いて次のよう

  

と近似し、

とおくと、ルンゲ・クッタ法を用いて、(1)式の近似解を求めることができる。

すなわち、

  

 

例によって、(2)を用いて、次の問題を解くと、右のグラフのようになる。

 

問題

  

を、初期条件

  

境界条件

  

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

 

 

Runge-Kutta-graph-001.png

 

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

 

Δt=0.1Δx=0.4とすると、さらに良好な結果が得られ、厳密解との差はほとんど認められなくなる。

 

Runge-Kutt-Graph-003.png

 

なのですが、ルンゲ・クッタ法も陽解法の一種なので、

  

のとき数値解が振動することがあるので、ルンゲ・クッタ法を用いて熱伝導方程式を解くとき、

  

になるようにΔtΔxを取ったほうがよい。

このことを示すために、Δt=0.36Δx=0.5のとき、すなわち、

として計算すると、40ステップのt=14.4の前後で振動を始め、t=18ではさらに激しく振動する。

Runge-Kutta-graph-004.png

Runge-Kutta-graph-005.png

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

 

 

! ルンゲ・クッタ法による解法
parameter (nt=20, nx=10)
real t(0:nx, 0:nt)
real dk(4,0:nx)
real l,kappa

l=4.; kappa=0.5
dt=0.2
dx = l/nx
dx2=dx*dx
c=kappa/dx2

dk= 0
   
! 初期化
do i=0, nx
    do j=0,nt
        t(i,j)=0.
    end do
end do

! 初期条件
do i=0,nx
    x=i*dx
    t(i,0)=x*(l-x)
end do

! 境界条件
do j=0,nt
    t(0,j)=0.
    t(nx,j)=0.
end do

! Runge-Kutta Method
do j=1,nt
    do i=1,nx-1
        dk(1,i) = dt*f(t(i-1,j-1),t(i,j-1),t(i+1,j-1),c)
    end do
   
    do k=2,3
        do i=1,nx-1
            dk(k,i)=dt*f(t(i-1,j-1)+dk(k-1,i-1)/2,t(i,j-1)+dk(k-1,i)/2,t(i+1,j-1)+dk(k-1,i+1)/2,c)
        end do
    end do

    do i=1,nx-1
        dk(4,i)=dt*f(t(i-1,j-1)+dk(3,i-1),t(i,j-1)+dk(3,i),t(i+1,j-1)+dk(3,i+1),c)
    end do
   
    do i=1, nx-1
        t(i,j) = t(i,j-1)+(dk(1,i)+2.*dk(2,i)+2.*dk(3,i)+dk(4,i))/6
    end do
end do

write(6,*) '*** Resutl ***'
do j=0,nt
    write(6,100) j*dt,(t(i,j),i=0,nx)
end do

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

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

end

function f(x,y,z,c)
f=c*(x-2.*y+z)
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


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

nice! 0

コメント 0

コメントを書く

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

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