［dy/dx＝0を通過する微分方程式の数値解］ その２

［dy/dx＝0を通過する微分方程式の数値解］　その２

 

 

　次にx＝0で解が止まらなかったとします。これは数値誤差によりy(0)＝ε≠0になったという事です。もしx＝0で解が止まるべきだという積極的理由があるなら、計算はx＝0で打ち切るべきです。じっさい今は(11)から、数値解はx＝0で止まっても良いとわかっています。しかし一般には、そんなこたぁ～事前にわかりません(^^;)。やってみて判断するのが普通です。だとすれば、y(0)＝ε≠0のεはルンゲクッタ法が許す精度内では妥当で、x＝0で解が止まらないのは正解かも知れないのです。もちろん疑ってかかるべきですが。

　ところがyの零点付近というのは、数値解法の生の誤差値が現れるところで、しかも相対誤差も意味をなさないので、εが非常に小さくても、それが本当に0かどうかの判定はかなり難しいのです。いずれにしろ今度もy＝x³にしか目を向けず、分岐解y(x)＝0を考えないというのは数学的に言ってまずいです。

 

　例を挙げますね。(11)は、ある力学問題を表していると仮定します。

　初期値y₀＝y(-1)＝-1，Δx＝2/80＝0.025で計算を開始したところ、y(0)＝ε≠0となりx＝0で解は止まらなかったとします。

　じっさいに計算すると、ε＝y(0)＝1.45×10^-5になります。Δx＝0.025の4次精度なので、yの精度はΔx³＝1.56×10^-5程度のはずです。これらを比べてみると、ほぼ同じです。

　これの意味は、4次のルンゲクッタ法では、ε程度の大きさは0と見分けがつかないという事です。従って正解がε＝0である可能性も非常に濃厚です。でも可能性は可能性です。

　いま我々は理論解を持っているので、初期値y(-1)＝-1を与えればy(0)＝0でなければならない事を知ってますが、普通はそうじゃありません。理論解を出せないからこそ数値計算します。そうするとy(0)＝0となるような初期値の与え方はy(-1-ε)であって、そのずれのためにy(0)＝ε≠0なのかも知れないからです。そうするとルンゲクッタ法は超正解を与えてる事になります。

 

　数学的には(x，y)＝(0，0)で解が分岐する可能性は否定できません。よって図-3のようになってしまった以上、y(0)＝0を初期値として(11)を数値積分した分岐解も探すべきです。ところがその計算は、蟻地獄です(^^;)。

 

　初期条件y(0)＝0から計算をスタートできるようにする、何らかのスターター手段が必要です。

 

　では物理的にはどうでしょう？。(11)は、ある力学問題を表しているのでした。分岐解のどちらかを選べるでしょうか？。力はyの2階微分に比例します。もしdy/dx＝0でもd²y/dx²≠0ならば、yはy(0)＝0から動くと言えます（前回は計算を間違えました(^^;)）。

　　

　(12)よりy＝0ではd²y/dx²＝0なので、速度0の時に力0なら動かないと結論できます。従ってy(0)＝0が真実なら、x＝0以後はy(x)＝0が正解です。しかし図-3は正しいのでしょうか？。その判断はできません。この段階は大学初年級レベルかな？(^^)。

 

　大学学部レベルになると、次のように言い出します。現実の物理系には減衰（摩擦）があるはずだ。よって減衰無しとして計算した結果がε＝y(0)＝1.45×10^-5と非常に小さいなら、そこで止まってしまうはずだと。何故なら(11)よりdy/dx＝0となるのは、y＝0の時のみだからy(0)＝0とみなすべきだ。よってx＝0以後はy(x)＝0が正解と。

 

 

　ところが院レベルになると、また趣向が変わります(^^)。

　摩擦があるから止まるって？。その結論は減衰定数の大きさに依存する。そんな事は、減衰項つきの微分方程式を解いてから言え。それよりも重要なことは、現実の物理系には必ず外乱がある事だ！。

　たとえx＝0以降の厳密解がy(x)＝0であったとしても(x，y)＝(0，0)は不安定停留点だ。それは位相空間で考えればわかる。

　外乱を認めるとx＝0でy＝dy/dx＝0に達しても、外乱によりy(0)は、微小に±δだけ、ほとんど必ずわけもなくジャンプすると考えられます。

　図-4に示すように、(11)のy(x)は常に単調増加です。もしジャンプ先が0＜δ＝y(0)なら、x＝0からy(x)は単調増加し続けます。よって図-3の赤点群は、定性的には正解という事になります。

　微小な外乱δは、微小という事以外は何もわかりません。従ってその理想化として、y(0)＝0を初期値とするy(x)＝0以外の数値解を計算する必要に迫られます。よってスターターが必要です。

 

　結局、数学科にいようと物理・工学系にいようと、スターターが必要になりました(^^;)。

 

　スターターとして自然なのは、件の微分方程式を高階微分方程式に持ち込むやり方です。単振動の時のように上手くいけば、解は自然にスタートしてくれます。じつは(11)に関してはdy³/dx³＝6です(^^)。3階微分まで考慮すれば、y＝0からでも数値解はスタートします。

 

　しかしこれには3つの問題があります。

　　1) 一般に高階微分であればある程、数値精度は悪くなり、計算手続きは面倒になり、計算量も増える。

　　2) 一般に何回微分すれば良いかわからない。

　　3) 最悪の事態として何回微分してもdⁿy/dxⁿが、その点で0かもしれない。

　　（※ あまり知られていない、テーラー展開不能点(^^;)）

 

　例えば(6)を(7)へ持ち込んだとたん、計算量は倍に増えました。3階微分なら3倍です。

 

　もっと泥臭いスターターを考えます。泥臭いという事は、強い方法だという意味です。単振動(6)では、y(0)＝1でなく、y(0.01)＝0.999950000416665を初期値とするだけで非常に高精度の解が得られます。

　そういうものですぐ思い付けるのは、x＝0からΔx離れたx＝Δxにおけるy(Δx)を、逆向きルンゲクッタ法に従ってy(0)＝0となるように逆算する手です。(11)の場合なら、y(Δx)＝y₁として、

 

となりますけれど、上記の非線形方程式をy₁について解くのはまぁ～、あきらめた方が無難ですな(^^;)。ところでΔxはかなり小さいのでした。という事はy₁もyの零点付近の値なので、かなり小さいはずです。Δx×y₁はすごく小さいと予想できます。Δx×y₁について線形化しましょう。

 

　上記近似は、Δxが小さいほど正確なはずです。そこでΔx＝0.001とすると、y₁＝1.259422×10^-10となり、ほぼ倍精度実数の実用精度の限界です。この結果をもとにΔx＝-0.001とした逆向きルンゲクッタを行うと、y(0)＝-1.71×10^-9になります。y(-1)＝-1を初期値とした時の結果、y(0)＝1.45×10^-5よりもだいぶ0に近いです。y₁＝y(0.001)＝1.259422×10^-10からΔx＝0.001でx＝0.025までルンゲクッタを行うと、y(0.025)＝1.459363×10^-5です。以後これを初期値とし、Δx＝0.025でx＝1までルンゲクッタで計算します。

 

 

　そして図-5をいっしょうけんめい観察するのです。そうすると、ある疑念がふつふつと沸いてきますきます(^^)。y(x)はx＝0に関して反対称ではないのか？（！）。じっさい黒点のy(1)はy(1)＝0.996619です。もしy(1)＝1が正解なら誤差は1/1000のオーダーで、これはy(0.001)を逆算した時の計算誤差の影響と考えられます。

 

 

　図-6は図-5の赤点のx≦0側の絶対値と黒点を比較したものです。両者の比はほぼ常に1です。x＝0近傍で1からはずれるのは、両者の値が小さくなりすぎて0/0状態になり比が意味をなさないからだと判断できます。y(x)は、x＝0を境に反対称である可能性濃厚です。もしもそうなら、y(0)＝0でなければなりません。もともとの条件(11)が、y＝0に関して対称だからです。

 

　こうしていちおう、数値的にも次の結論が得られます。

　(11)のy(x)は、初期値y(-1)＝-1に対してy(0)＝0で数学的には解が分岐する。物理的には、どちらを選ぶかは状況による。

 

　でもあくまで「いちおう」です(^^;)。結論を確定させるには、もっと高精度の力業が・・・(^^)。

 

 

ネムネコのつぶやき

 

あくまで、微分方程式

　　

を数値的に解く場合に限っては、テーラー展開をそのまま利用し、難しいことは何もせずに、

　　

と逐次的に計算するほうがいいみたいですよ(^^ゞ

 

　　

だから、４次以上の高次微分は０。

 

昨日（６月２日）に気づいたのですが、

考えてみれば、微分方程式のyの解は３次関数なのだから、⑨式で正確に計算できなければおかしいんですよね。

 

ルンゲ・クッタ法で正確に計算できないのは、打切誤差や丸め誤差のために、x=0近傍での２次、３次と４次の高次微分を正確に計算できないためなんだと思います。

x=0近傍におけるルンゲ・クッタ法の数値解を元に、差分法を使って、x=0近辺での高次微分の値を計算すると、おそらく、そうなっているのではないでしょうか。

そして、これがx=0近傍での解の精度を悪化させ、x>0における計算精度を悪化させていると思います。

 

⑨式を用いて、表計算ソフトを用いて、Δx=0.1で計算した結果は次の通り。

 

 

そして、注目に値するのは、テーラー展開に基づくこの計算法だと、y=x³という解を返してくれることです。

オイラー法やルンゲ・クッタ法などの数値積分法だと、(x,y)=(0,0)にという正しい答を計算してしまうと、それ以降は蟻地獄状態に嵌り込んで、y=0という解になってしまいますから。

 

微分（級数解）と（数値）積分の違いが出ていて、非常に興味深い事例ではないでしょうか。

どちらも、テーラー展開を基礎にしているんですが・・・。

 

「４次のルンゲ・クッタ法の計算結果は正確である」という思い込み、経験則が、実は、「落とし穴」だったのかもしれませんね。

 

なお、電卓メーカーのCASIOさんに、高精度計算のための次のサイトがあります。

https://goo.gl/YTC6mU

 

ここで、Δx=0.004で計算してもらったところ、

 

 

 

十進数で１８桁くらい精度をもっているらしいですが、x=1のとき、1%くらいの誤差が出るみたいです。

 

 (x,y)=(1,1)を初期条件にし、x=−1まで逆向きに計算した結果と上の計算結果を足して２で割ると、おそらく、正しい計算結果になると思いますが（笑）。