Excelな工学

2階微分方程式の境界値問題を解くマネゴトをしてみる

作成:2001-01-24

更新:2001-01-25


そもそも境界値問題って何? っていう疑問もあるかとは思いますが 追々説明するとして、まずは
============================================================
★1階微分方程式

Y'=f'(x)=6X+2
を x=0.0〜1.0の間で解いてみます。
f(0)=C=1
が与えられているとします。

Y(n)=f(X(n)) Y(n+1)=f'(X(n+1))
だとすると、微分を差分で近似すると

Y(n+1) - Y(n)
------------ = Y'(n)
X(n+1) - X(n)

ここで、ΔX=X(n+1)-X(n)は常に等間隔だとすると

Y(n+1) - Y(n)
------------ = Y'(n)
   ΔX

だから
Y(n+1) = Y(n) + ΔX*Y'(n) = Y(n) + Δx * (6*X(n) +2 )

Y(0)が与えられているので
Y(1)=Y(0) + ΔX*Y'(0) = Y(0) + Δx * (6*X(0) +2 )
と、順番に求まって行きます。

Excelで解くと

Y(n) = Y(n-1) + ΔX*Y'(n) = Y(n) + Δx * (6*X(n) +2 ) だから、X=0.0 の時にY=1.0としておいて、
Y(近似)のセルは、左隣のセル + ΔX * (6*Xのセル+2)
ΔX=0.1
 ってな感じ。

代数的に解けば Y=f(x)=3x^2+2x+C なので、
初期値としてf(0)=C=1ならば 解はY=3x^2+2x+1 (これがY(真値))

真値との間に誤差があるのは、微分を差分で近似しているからで、もっと精度よくやりたい人は「ルンゲクッタ法」とか、その辺を参考にして下さい。(ルンゲクッタに関しては、数値解析の本に詳細に書いてあります。)Excelで簡単にやるには、こういう線形近似(オイラー法)の方が手軽だと思う。

============================================================
★2階微分方程式


★初期値値問題

一般に微分方程式の数値解析は初期値問題と境界値問題に
分かれます。1階の時は初期値=境界値だったわけで。

 微分方程式:Y'' = 6 を解く。
2階微分方程式を差分で近似して
一般に

Y(n+1) - Y(n)
------------ = Y'(n)
   ΔX

Y'(n+1) - Y'(n)
------------ = Y''(n)
   ΔX


Y'(n+1) = Y'(n) + ΔX*Y''(n)
Y (n+1) = Y (n) + ΔX*Y' (n)

この漸化式が使えるのは Y(0)=1 , Y'(0)=2 とか が与えられている時で、例えば

重力加速度の式が与えられていて、初速度と初期位置が与えられている、
X''=G=-9.8m/s2 , X'(0)=V(0)=0 X(0)=0
っていう、高校物理の最初の問題ですな。

これはExcelだと、上の表を二つ重ねてしまえばOKっすね。

Y'(n+1) = Y'(n) + ΔX*Y''(n) Y'(0)=2 , Y(0)=1
だから、X=0.0 の時にY'=2.0 , Y=1.0 としておいて、
ΔX=0.1で、
Y'のセルは、左隣のセル + ΔX × (Xのセル×6)

 

Y (n+1) = Y (n) + ΔX*Y' (n)
Yのセルは、左隣のセル + ΔX *(Y'のセル)


★境界値問題

 

ところが、実際の問題では、1階微分の初期値Y'(0)が与えられているんじゃなくて 解析区間の両端、 Y(0) = 1 , Y(1) = 6 が与えられていて、この条件の下に Y'' = 6 を解く、なんて必要がある時がある。
具体的には、熱伝導方程式で、両端の温度が分かっているので、間の温度分布が 知りたいとか、流体方程式とか、なんかそういう場合ですな。

これを解く漸化式を考えてみると(実は漸化式にならないんだけど) 微分の差分近似を、もう一回繰り返せば2階微分になるので 下記の差分で考えてみます。

Y(n+1) - Y(n)     Y(n) - Y(n-1)
------------ - --------------
    ΔX         ΔX
-------------------------------- = f''(X)
          ΔX

∴Y(n) = { Y(n+1) + Y(n-1) - Δx^2*f''(X(n)) } / 2

これってどうやって求めればいいんでしょか?
n=1とするとΔx=0.1として
y(1) = { y(2) + y(0) - 0.01*6 } / 2
なので、 y(1)を求めるには、y(2)を求めなければならない。

 で、普通ならここで諦めてしまうのですが、無理やりExcelの計算式に入れてしまうと 「循環参照」とかいうエラーになります。


これは、自分自身を計算式に入れてしまう時に出てくるエラーです。 (例えば A1B1 = A1B1 + A2B1 とか )

ここで、「反復計算」のオプションをつけます。
ツール/オプション/反復計算 最大反復回数=1000 変化の最大値=0.0001

あとは収束するまでひたすらF9を叩けばよろし。 すると収束して、理論値通りの解が求まります。 (変化の最大値の値を小さくしながら。但し0にするといつまでたっても計算が終わらない)

「境界値問題を解くマネゴト」というのは、敢えて収束する様な程式と近似式を使ったからであり、この方法が常に収束するとは限りません。
境界値問題の収束計算と収束加速に関しては、それだけで学術分野を持つ奥の深いものなのですが、どうせ自分のPCなので、収束しなかったとしても誰に迷惑かける わけでもなし、CPU時間を無駄食いして1行1行マークシートを塗ってコーディングしたのが全部無駄になって、その挙句、計算機センターの係員からブツブツ言われる心配も無いですし(っていつの時代の話をしてるんだか(笑))
なんかExcelのセルが収束していくのって、見ていて気持ち良いっすね(笑)

 


return