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

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

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

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

【解】

微分方程式

の特性方程式は

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

（１）　境界条件が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で計算したものを以下に示す。