今日のアニソン、「鹿楓堂よついろ日和」から『桜色クリシェ』

今日のアニソンは、アニメ「鹿楓堂よついろ日和」から『桜色クリシェ』です。

このアニメは人気がないのだろうか、YouTubeにはこれしかなかったケロ(>_<)

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

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

 

 

　ここでは次の形の微分方程式を考えます。

　　

　まぁ～典型的な変数分離形ですね。で(1)に対して、4次のルンゲクッタ法をご紹介申しあげます。「ルンゲクッタ 4次」でググるとごまんとヒットしますから、それらからの引き写しと思って下さいな(^^;)。自分はここ、

　　https://www.sit.ac.jp/user/konishi/JPN/L_Support/SupportPDF/Runge-KuttaMethod.pdf

を愛用しとります。

　　

というyの更新手続きが、4次のルンゲクッタ法です。Δy_jをyの増分と言います。(2)の最後の式でyの増分Δy_jが計算されますが、k₀～k₃は(1)の場合、以下です。

　　

　なんとも単純じゃねぇ～ですか(^^)。しかもExcelにすぐ載ります。しかもしかも非常に高精度である事が知られています。例えば、

　　

 

を解いてみます。初期条件はy₀＝y(0)＝1として、x＝0～1の範囲で数値積分しy(1)の値を求めます。Δxはx＝0～1を10分割し0.1とします。y(1)＝y₁₀になります。理論解はy(1)＝eです。

 

 

　ちなみに、

　　e ＝2.71828182845905

　　y₁₀＝2.71827974413517

となり、|e－y₁₀|/e×100＝0.00008 % ～ 1/100万の相対誤差しかありません。このように4次のルンゲクッタ法さえ使っておけば、たいがいの微分方程式はなんなく解けます。

 

　本題に入ります。(1)でy₀を、たまたまdy/dx＝f(y₀)＝0に選んだらどうなるか？という話です。ルンゲクッタ法の実質計算部(3)に、この条件を適用します。

 

 

 

　従ってk₀＝k₁＝k₂＝k₃＝0。(2)に戻ってΔy₀＝0。なのでy₁＝y₀＋Δy₀＝y₀。xだけがx₁＝x₀＋Δx＝Δxと更新されます（x₀＝0）。yの値は変わりません。x₁＝Δxにおいて同じ事をやれば、yの値は変わってないので、y₂＝y₁＝y₀，x₂＝2Δxです。よってx＝0，Δx，2Δx，3Δx，4Δx，・・・において、y₀＝y₁＝y₂＝y₃＝y₄＝・・・となり、ルンゲクッタ法は定数関数y(x)＝y₀を計算します。これがネコ先生の仰る、蟻地獄です(^^;)。

 

　この結果は正しいでしょうか？。(1)のタイプの微分方程式は、dy/dx＝0の点に達した瞬間に、以後は定数関数y(x)＝y₀にしかならないのでしょうか？。微分方程式とは、直前の値が現在の値を決めるものだとすれば、これ以外の解はあり得ません。しかしちょっと（というかかなり）不味いんじゃないか？と思える例を、次にあげます。

 

　単振動もじつは、(1)の形で表せます。次式は単振動のエネルギー積分です。

　　

　(5)をdy/dxについて解いて、±の－の方を取ると、

　　

と(1)の形で表せます。(6)に初期条件y(0)＝y₀＝1を与えてやれば、y(x)＝cos(x)が計算されると考えがちですが、y＝1ではdy/dx＝0なので、さっきの話から定数関数y(x)＝1にしかなりません。しかもこれは(6)の解の一つです。y(x)＝1を(6)に代入すれば、常にdy/dx＝0なので。

　では4次のルンゲクッタ法がおかしいのか？。おかしくないのです。4次と言うだけあって4次のルンゲクッタは、4次関数までは厳密に計算できます。定数関数なんて屁の河童です。

　ところが初期条件として、y(0.01)＝y₀＝0.999950000416665なんかを与えてやると、突如さっきと同じように1/100万の精度でy(x)＝cos(x)を再現します。

 

　この蟻地獄は、(6)を2階微分方程式に直すことによって抜け出せます。(6)をxで微分して、微分した結果に(6)を考慮すると、

　　

です。(7)をルンゲクッタに持ち込むには、dy/dx＝pと書くことにして、

　　

とすれば十分です。(8)の一段目をxで微分し二段目に代入すれば、(7)に戻れます。この場合のルンゲクッタ法の定式化は、yとpの相互作用に注意しつつ一段目と二段目のそれぞれに(2)，(3)を適用すれば得られます。

　　

　　

 

 

　初期条件y(0)＝y₀＝1，p(0)＝p₀＝0、Δx＝2π/20＝π/10で(9)，(10)の計算を実行したのものが、図-2です。(8)による解は、きれいにy(x)＝cos(x)にのっています。しかも(7)に対して(8)を考え、そのルンゲクッタ法の定式化を行っただけで、自然にそうなりました。その原因は明らかですよね？。

　(8)ではp＝dy/dxがyで駆動されてるからです。それは(7)の形をみれば一目瞭然です。dy/dxの傾きまで考慮してるので、結局dy/dxは自律的に0から変化します。

 

　ここで再び、ネコ先生の仰る哲学問答が発生します。(6)を(7)の形で考えるという事は、単振動の運動方程式を与えたのと同じです。しかし微分方程式とは幾何学的関係を与えたものに過ぎなかったのではないか？。(6)のどこに、これは力学問題であると書かれているのか？。(6)～(10)の過程でいえる事は、たんに(6)はy(x)＝1で分岐解を持つ、ただそれだけではないのか？。(7)，(8)を選ぶべきだという積極的理由がどこにある？。

 

 

　・・・それがあるんです。それは最初に、(6)は単振動の微分方程式だと言い切ったからです(^^)。単振動は(7)で定義されます。

 

　以上を念頭に、

　　

に戻ります。ネコ先生と同じように、

　　http://nekodamashi-math.blog.so-net.ne.jp/2018-05-25-2

y₀＝y(-1)＝-1を初期値に取り、Δx＝2/80＝0.025とします。まず、数値解がdy/dx＝0から先はy(x)＝定数になって駄目なんですか？（蓮舫さんじゃないけど(^^)）が、基本的なところだと思えます。だって微分方程式は幾何学的関係を与えるだけのものなんだから、これはこれで正しく解けてるじゃないかと。

 

　・・・う～ん(^^;)。微分方程式が色んな場面に応用される事を考えると、たぶん一般的には芳しくないんですよ。y(x)＝0という定数解は現象として全然面白くない。むしろ「応用上は」x＝0でのもう一方の分岐解y＝x³を期待して、計算をするほうが遥かに多いと思うんですよね。いずれにしろdy/dx＝0の点で分岐解を全て探せないというのは、数学的に言ってまずいです。

 

（執筆：ddt³さん）