オレって実はスゴイ（笑） Part２ [ひとこと言わねば]
ネムネコが、昨夜、閃いた（思いついた）、４次のルンゲ・クッタ法を用いた２階常微分方程式の境界値問題の解法は、第一種の境界条件と第２種の境界条件の両方が使用されている場合にも、適用が可能なようだね。
h=0.1で、上の境界値問題を解いた場合の計算結果は次の通り。
参考までに、h=0.01とし、差分を用いた結果も以下に示す。
h=0.01とルンゲ・クッタ法の１０倍計算格子の間隔を細かくとっているものより、心もち、高精度で計算できていることがわかる。
ではあるが、
混合条件にまで拡張できるか、どうやって拡張したらよいかについてはまったくの白紙状態。
大体、ネムネコは、数値計算に関しては素人だにゃ。
化けネコではあるが、ネコがここまでできるってこと自体が既に驚異だにゃ。
そう思わないかい？
今日のアニソン、「名探偵ホームズ」から『グッバイ・スウィートハート』 [今日のアニソン]
差分法による２階常微分方程式の境界値問題の数値解法３ [数値解析]
差分法による２階常微分方程式の境界値問題の数値解法３
前回、次の微分方程式
のx=1における第２種の境界条件（ノイマンの境界条件）
を後退差分
と近似し計算したところ、うまくいかなかった。（下図参照）
これは、
（１）式を差分方程式にするときに用いた
であるのに対し、x=1における境界条件の近似で用いた
の誤差がO(h)と、（２）、（３）式に対して誤差が大きく、これが数値解の精度を悪化させたためである。
したがって、この問題を解決するためには、x=1における境界条件の近似を（４）より高精度にする必要がある。
そこで、
と[0,1]の外に仮想の点
を設け、さらに、微分方程式
がx=1でも成立すると仮定する。
すると、x=1における差分方程式は
になる。
x=1における境界条件に中心差分を用いると
これを（５）に代入すると、
これと、k=1〜n−1に関する差分方程式
を合わせると、未知数の個数は
のn個、対して、方程式の数もnとなるので、次の連立方程式を解くことによって、を全て求めることができる。
この方針に従ってプログラムを書き換え、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 ! nは分割数
a=0.; b=1. ! 境界のx座標
x(0)=a; x(n)=b
y(0)=1. ! 境界条件
call Solver(x,y,n)
! 計算結果の出力
write(*,*) ' i x y Exact'
do i=0, n, 1
write(*,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=2/(2-e)*exp(-x)- e/(2-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
dx2=dx*dx
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-dx2*g(x(i)))
c(i)=1+dx/2.*f(x(i))
d(i)=dx2*h(x(i))
end do
! 境界条件
d(1)=d(1)-a(1)*y(0) ! x=0 ディリクレ条件
! x=1 dy/dx=0 ノイマン条件
a(n)=2.
b(n)=(dx2*g(x(n))-2.)
c(n)=0.
d(n)=dx2*h(x(n))
call tdma(a,b,c,d,n) ! TDMAで連立方程式を解く
do i=1,n
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.3近傍、y=0付近での数値解と厳密解との差が多少気になるので、n=100、h=0.01で計算した結果を以下に示す。
h=0.1のとき、小数点２桁目に、h=0.01の場合は、小数点４桁目で、数値解と厳密解が異なっているが、数値計算の精度がこの程度だから、これはしょうがない。
オレって実はスゴイ（笑） [ひとこと言わねば]
ルンゲ・クッタ法（４次）で境界値問題
が簡単に解けてしまう(・・?
ルンゲ・クッタを３回使うだけで、何故か分からないけれど、簡単に解けてしまったケロよ。
この解法を思いついたネムネコがすごいってことか（笑）。
何でこんなに簡単に解けるだろう。不思議だな〜。
