差分法による２階常微分方程式の境界値問題の数値解法２

差分法による２階常微分方程式の境界値問題の数値解法２

まず、いきなり、簡単な微分方程式の問題！！

問題　次の微分方程式を解け。

【解】

微分方程式

の特性方程式は

したがって、①の基本解はなので、①の一般解は

（１）　境界条件がy(0)=1,y(1)=0なので、②より

これを解くと

したがって、

（２）　②より

境界条件はy(0)=1、y'(1)=0だから

これを解くと

よって、

（解答終）

こんな問題を解きたいわけではなく、問題の（２）を差分法を使って解くことが今回のテーマ。

いきなり、（２）を解くのは大変なので、問題の（１）について考えることにする。

閉区間[0,1]をn等分し、

さらに、

と書くことにする。

差分法を用いると、

と近似することができるので、微分方程式

の近似式は

 

になる。

未知数は、であり、連立方程式（１）の式の本数はn−1だから、連立方程式（１）を解くことによって、微分方程式を数値的に解くことができる。

プログラムを新たに作ってもいいけれど、過去に作ったより一般的な

を解くプログラムを再利用し、[0,1]をn=10、h=0.1と１０分割し計算した結果は次の通り。

n=10、h=0.1と粗い計算でも、微分方程式の解を正確に再現していることがわかるだろう。

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

計算に使用したプログラム



! 差分法を用いて y''+f(x)y'+g(x)y=h(x) の境界値問題（ディリクレ条件）を解くプログラム
! ただし、固有値問題は解けない!!
parameter (ns=100)
real x(0:ns), y(0:ns)

x = 0.; y =0;

n=10

a=0.; b=1. ! 境界のx座標
x(0)=a; x(n)=b
y(0)=1.; y(n)=0. ! 境界条件

call Solver(x,y,n)

write(6,*) ' i      x         y'
do i=0, n
    write(6,100) i, x(i), y(i), exact(x(i))
end do

100 format(i3,1x,f9.5,1x,f9.5,1x,f9.5)
end

function exact(x)
e=exp(1.)
exact=1/(1-e)*exp(-x)- e/(1-e)*exp(-2.*x)
end

function f(x)
f=3.
end

function g(x)
g=2.
end

function h(x)
h=0.
end

! ここより下はいじると危険
! ブラックボックスとして使うべし


subroutine Solver(x,y,n)
real a(n),b(n),c(n),d(n)
real x(0:n), y(0:n)

n1=n-1
dx=(x(n)-x(0))/n

do i=1,n1
    x(i)=x(0)+i*dx
end do

! 差分法によって得られる連立方程式の係数の計算
do i=1,n1
    a(i)=1-dx/2.*f(x(i))
    b(i)=-(2-dx*dx*g(x(i)))
    c(i)=1+dx/2.*f(x(i))
    d(i)=dx*dx*h(x(i))
end do
    d(1)=d(1)-a(1)*y(0) ! 境界条件
    d(n1)=d(n1)-c(n1)*y(n) ! 境界条件

call tdma(a,b,c,d,n1) ! TDMAで連立方程式を解く

do i=1,n1
    y(i)=d(i) ! 計算結果をセット
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
 

このプログラムを利用すれば、問題の（２）の微分方程式の境界値問題も簡単に解けそうに思うだろう。

しかし、そうは問屋が卸してくれない。

問題の（１）と（２）では、x=0の境界条件がy(1)=0とy'(1)=0と、境界条件の種類が違うため、（２）ではの値が与えられていないからだ。

したがって、を求める方法をあらたに考えないといけない。

この最も簡単な解決方法は、境界条件

を後退差分を用いて

と書き換えて、式（１）に式（２）を新たに加えて解くというもの。

これで未知数の数がn−1、方程式の数がn−1となり、解くことができるはずである。

この考え方に基づいてプログラムを書き換えて、n=10、h=0.1の条件で解いたものが次の通り。

厳密解と数値解が一致しないので、n=100、h=0.01の条件で計算したものは次の通り。

これくらい計算格子を細かくとれば、それなりに良好な計算結果が得られるが、このような小手先の変更では、実用に足りないことがわかるだろう。

問題の（１）のような境界条件をディリクレ条件、（２）のx=1のときのように１階の微分で境界条件が与えられるものをノイマン条件というが、数値計算でノイマン条件を入れるのは、結構、厄介。

　――場合によっては、プログラムを全面的に書き換えることが迫られる！！――

（１）式と（２）式は、必ずしも、整合していないし(^^ゞ

参考までに、問題（２）のx=1における境界条件を

とディリクレ条件にし、n=10、h=0.1で計算したものを以下に示す。

定常１次元熱伝導方程式の数値計算

定常１次元熱伝導方程式の数値計算

熱伝導率λが温度Tなどによって変化する場合、非定常の１次元熱伝導方程式は次のようになる。

ここで、ρは密度、cは比熱、は単位時間単位体積あたりに発生（消滅）する熱量である。

したがって、熱の発生がない場合、１次元の定常熱伝導方程式は

となり、境界条件が与えられれば、この解を求めることができる。

たとえば、熱伝導率λが一定で、x=0における温度がT₀、x=lにおける温度がT₁の場合、

さらに、単位時間単位面積当たりに通過する熱量（熱流束）は

である。

熱伝導率が一定の場合、温度Tは直線的に変化するので、わざわざ数値計算をする必要はない。

そこで、熱伝導率λが温度によって変化する場合の定常１次元熱伝導方程式の差分法を用いた数値解法について考えてみる。

（２）は

となるので、が

の場合、差分法を用いて

と近似することが可能。

しかし、プログラムがすこし複雑になるので、今回、この計算法は採用しないことにする。

そこで、今回は、

という差分を用い、（２）式を

ここで、とは、それぞれ、との中点における熱伝導率。

それで、

 

そして、熱伝導率λは

として、これを解くプログラムを作ってみた。

ちなみに、この微分方程式の解は

である。

計算領域0≦x≦1をn=１０分割した計算結果は、次の通り。

青い直線は熱伝導率が一定の場合。

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

parameter (n=10)

real t(0:n), gam(0:n)

real a(n),b(n),c(n),d(n)

eps=1.e-6

dx =1./n

! 変数の初期化

t=0.; gam=1.

a=0.; b=0.; c=0.; d=0.

! 境界条件

t(n)=1.

do k=1, 10

! 熱伝導率Γを計算

do i=1, n

gam(i)=1+t(i)

end do

! 連立方程式の係数を計算

do i=1, n-1

dx2=dx*dx

aw=0.5*(gam(i-1)+gam(i))/dx2

ae=0.5*(gam(i)+gam(i+1))/dx2

ap=aw+ae

a(i)=-aw

b(i)=ap

c(i)=-ae

d(i)=0.

end do

a(1)=0.; d(1)=0.5*(gam(0)+gam(1))/dx2*t(0)

c(n-1)=0.; d(n-1)=0.5*(gam(n-1)+gam(n))/dx2*t(n)

! TDMAで連立方程式を解く

call tdma(a,b,c,d,n-1)

! 計算結果の更新と収束判定

err=0.0

do i=1,n-1

err=amax1(err,abs(1-t(i)/d(i)))

t(i)=d(i)

end do

if(err.lt.eps) exit

end do

do i=0,n

x=i*dx

write(*,100) x, t(i), -1+sqrt(1+3.*x)

end do

! ファイルへの出力

open(1, file='Netsu.dat', status='replace')

do i=0,n

x=i*dx

write(1,100) x, t(i),-1+sqrt(1+3.*x)

end do

close(1)

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

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

なお、このプログラムでは、とを

と、点における熱伝導率の算術平均を用いて計算している。

熱伝導率のこの計算の是非については議論が必要だろうが、計算結果は極めて良好なのだから、「終わりよければ全てよし」ということで(^^)

本問は、いちおう、Tの非線形の微分方程式なので、計算においては、

と、反復１回前の温度を用いて熱伝導率を計算し、線形化して解いている。

反復計算の収束条件は