微分方程式の境界条件についてのよもやま話 [ひとこと言わねば]
微分方程式の境界条件についてのよもやま話
たとえば、次のような２階の常微分方程式があるとする。
この微分方程式の境界条件が
と与えられるものをディリクレの境界条件（第１種の境界条件）、
と１階微分の形で与えられたものをノイマンの境界条件（第２種の境界条件）などという。
今日の数学の記事で書いたけれど、ノイマンの境界条件をうまくプログラムに入れるのは、結構、厄介。
単に面倒なだけではなく、数値解はノイマン条件に、結構、敏感に反応するので、とんでもない結果をもたらすことがあるにゃ。うまく入れないと下の図のように数値計算と厳密解が大きく食い違ったりする。
ノイマンの境界条件をプログラムにうまく入れることが如何に大変か、
この件についてddt³さんがきっと何か素敵な記事を書いてくれると思うにゃ。みんな、期待するにゃ。
微分方程式の境界条件は、ディリクレ、ノイマンの境界条件の他に、混合型の境界条件（第３種の境界条件）という更に厄介なものがある。
非定常１次元の熱伝導方程式
の場合は、
や
といった奴だにゃ。
（１）でも、結構、大変なのに、（２）なんかは涙がチョチョギ出るほど大変だと思うね。
（２）は、プログラムの部分的修正ではダメで、計算に使うアルゴリズム、使用する計算法などを全面的に見直し、プログラムを一から作り直さないといけないかもしれない。
そして、
新たに作ったプログラムによる計算結果が現実と大きく食い違うなんて可愛い話ではすまず、計算方法――計算に使うアルゴリズム、スキームなど――によっては、数値解が簡単に発散したり強く振動したりして、収束解が得られないなんて事態に遭遇するかもしれない(>_<)
⑨の定常解、つまり、
である、
の場合ですら、４次の代数方程式を解かなければ、この微分方程式は解くことができないんだから、（２）の条件を入れて解くことが如何に大変か容易に想像できるだろう。
差分法で解こうとする場合、最も簡単な方法でも、境界条件（２）から
という４次方程式が出てくる！！
単体の方程式（３）だけならば、ニュートン法や２分法を用いて、この近似解を求めることができるだろうが、
というn−1本の線形の連立方程式（４）に非線形の方程式（３）を同時に解かないといけないのだから、これはそんなに簡単な話ではない。
ひょっとしたら、（３）と（４）からなる連立方程式系は実数解をもたないかもしれない(^^ゞ
４次方程式（３）は上のグラフのような状態になっており、解けないかもしれない(^^ゞ
２種や３種の境界条件をうまく計算プログラムに入れるってのは、それほど、大変なんだケロ。
そして、多くの数値解析、数値計算の（入門的な？）本には、こういった条件の取り扱い方については、ほとんど書いていないのであった。
英語で書かれた「Computer Analisys」といったようなタイトルの１万円くらいの、結構、専門的な英語の本であっても、この手の話はほとんど出ていないと思うよ。
ってことは、ネムネコが、いま、書いている記事には、１万円以上の価値があるってこと(・・?
これらの記事を読んだヒト、一人から３００円ではなく千円くらいのお賽銭をもらってもバチは当たらないに違いない（笑）。
今日のアニソン、東方から『幻想河童行進曲』 [今日のアニソン]
ねこ騙し数学、海外のサイトから攻撃を受ける！！ [ひとこと言わねば]
差分法による２階常微分方程式の境界値問題の数値解法２ [数値解析]
差分法による２階常微分方程式の境界値問題の数値解法２
まず、いきなり、簡単な微分方程式の問題！！
問題 次の微分方程式を解け。
【解】
微分方程式
の特性方程式は
したがって、①の基本解は
なので、①の一般解は
（１） 境界条件が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=1における境界条件がy(1)=0とy'(1)=0と境界条件の種類が違い、（２）では
の値が与えられていないからだ。
したがって、y'(1)=0から
を求める方法をあらたに考えないといけない。
この最も簡単な解決方法は、x=1における境界条件y'(1)=0を、後退差分を用いて
と書き換え、（１）式に（２）式を新たに加えて解くというもの。
これで、未知数の数がn、連立方程式の数がnと一致し、解くことができるはずである。
この考え方に基づいてプログラムを書き換えて、n=10、h=0.1の条件で解いたものが次の通り。
解くことはできたが、厳密解と数値解が一致しないので、n=100、h=0.01の条件で再計算したものは次の通り。
これくらい計算格子を細かくとれば、それなりに良好な計算結果が得られるが、このような小手先の変更では、実用に足りないことがわかるだろう。
ちなみに、
問題の（１）のような境界条件をディリクレ条件、（２）のx=1のように１階の微分で境界条件が与えられるものをノイマン条件というが、ノイマン条件をプログラムにうまく入れるのは、結構、厄介なのだ。
そもそも、（１）式と（２）式は、必ずしも、整合していないし・・・。
参考までに、問題（２）のx=1における境界条件を
とディリクレ条件に直し、n=10、h=0.1で計算したものを以下に示す。
