今日のアニソン、けものフレンズのみんなで『恋』ダンス

今日のアニソンは、けものフレンズのみんなで『恋』ダンスです。

「この曲はアニソンじゃない」なんて、かたいことは言わず言わず、この動画を愉しんで欲しいにゃ。
さらに、けものフレンズによるこの動画を♪

［2階常微分方程式（保存系）］

［2階常微分方程式（保存系）］

 

　　

　(1)は保存系の微分方程式です。保存系だと非線形であってもエネルギー積分があるので、それなりに出来る事があります(^^)。

 

　(1)の両辺にy'をかけ、tで積分すればエネルギー積分です。

　　

から、合成関数の微分公式の逆を使って、

　　

となり、エネルギー積分を得ます。ここにＣは任意定数。初期値問題ならばt＝0でのy'とyがわかっているのでＣを計算できますが、今は境界条件が(2)なのでＣを直接計算できません。それでもy(0)＝0を(3)に代入すれば、

　　

となり、Ｃが初期速度のかわりになってるのがわかります。Ｃを決める条件はもちろん(2)のy(1)＝1です。従ってＣを含んだ形で問題を解いた後でないと、Ｃを決定できないのもわかります。Ｃを決定できないという事は、未完の解答という事です(^^;)。

 

　ところが1次元の保存系の形式解は、必ずあるんですよ。(3)が変数分離形である事はすぐわかります。(3)をy'について解けば、

　　

ですから、

　　

です。(6)をtで積分すれば、再び合成関数の微分公式の逆を使い、

です。

　±1の積分は明らかに±t＋Ｄですから（Ｄも任意定数）、(6)，(7)より、

　　

が得られます。これにy(0)＝0とy(1)＝1を代入すれば、

　　

となるので、(9)，(10)の連立方程式からＣ，Ｄを決定して(8)に戻れば良い訳です。よって(8)中辺の具体的な形がわかればOKですが、たぶん駄目です。(8)中辺の積分は基本的に、

　　

がわかればいいのですが恐らく無理ですよね、ネコ先生？。岩波公式集にも載ってなかったもの(^^;)。

 

　そういう訳で定性的な検討に移行します。まず(3)の、という事は(5)のグラフを描いてみましょう。(5)の√内正より、y≦(3Ｃ)^1/3でなければならない事がわかります。その事に注意してExcelなんかでグラフを描けば図-1になります。y(t)がy(1)＝1に達するためには、1/3≦Ｃが必要です。

 

 

　問題の条件は、時間の経過tの増加とともに、y(t)はy(0)＝0 → y(1)＝1と変化するというものでしたが、図-1に示したように、y＝0 → 1となる位相空間での軌道はいくつもあります（図-1の赤矢印）。しかしそれは変です。何故ならy(0)＝0とy(1)＝1は、t＝0においてyとy'を与える普通の初期条件のかわりをしてるはずだからです。初期条件を与えれば、常微分方程式の解は一意に存在します。それが常微分方程式の解の存在定理です（これを本気で証明しようという人は、ほとんど見た事ありませんが(^^;)）。

 

　何かを見落としてます。答えを言うと、時間経過です。つまり図-1の3つの赤矢印のy＝0での時間t₀とy＝1での時間t₁との差は、t₁－t₀＝1と限らないという事です。初期条件を与える場合は同時刻でyとy'を与えるので、ふだんは経過時間なんか気にしませんが(^^;)。

　逆に言えばt₁－t₀＝1であるという制約条件が、位相空間上に無数に存在する、y＝0 → 1となる軌道から正しい軌道を選択してくれます。じっさい図-1の3つの軌道は、y(t₀)＝1ではありますが、y'(t₀)はそれぞれ違います。違った初期条件のもとでは位相空間での軌道は必ず違う。それが常微分方程式の解の存在定理の言ってる事です。だから違った3つの軌道になっていいです。と同時に、y＝1となる時刻t₁はそれぞれ違うしかない事にもなります。何故なら図-1の曲線は、微分方程式(1)そのものだからです。

　しかし今まで導いた条件に、y＝0 → 1のための時間経過を教えてくれるようなものはあるんでしょうか？。　それが保存系の場合は「あるんです！」(^^)。それは式(8)そのものです。この議論を一般化すると、作用変数・角変数とかハミルトン・ヤコビ方程式の話になります。そういう訳で(8)中辺の具体的形がわかればいいんですけど、・・・けっきょくそれかい！(^^;)。

 

　再び出来ない話に帰着したので、再び定性的検討に移行しましょう(^^;)。今度は(8)中辺の被積分関数をg(y)として、g(y)のグラフを描いてみます。図-1の逆数のグラフを描くだけですよ。図-1から0≦y'だけで良いので、

 

 

です。(8)で任意定数Ｄは明らかに、時間の原点の設定を表しています。時間の原点の設定は任意にできるので、Ｄ＝0としましょう。そうすると(8)の積分は、図-2の曲線のy＝0～1における面積そのものになります。それがt₁－t₀＝1になる事が、制約条件です。1/y'をyで積分した結果の面積をＴとします。図-2の黄色いハッチ部は、y＝0～1における面積1を表します。

　図-2を眺めると、Ｃ＝0.5では恐らく1＜Ｔ，Ｃ＝0.7ではＴ≦1か1≦Ｔかは不確定、Ｃ＝1なら確実にＴ＜1です。またＴ(Ｃ)はＣに対して単調減少であろうという予想もできます。もしそれが正しいなら、二分法かなにかで数値的にＣを決定できます。Ｃが決定できれば(4)から、y'(0)もわかります。

 

　Ｔ(Ｃ)がＣに対して単調減少である事を示します。Ｄ＝0とすれば(8)から、

 

　　

　積分記号化の微分を使って、

　　

が得られます。図-2の条件下では被積分関数は正で、dＴ/dＣは負。yがy＝0から1にいたる経過時間Ｔ(Ｃ)は、Ｃに対して単調減少とわかります。こうなれば後は台形公式でも使って(11)の面積を計算し、Ｃを探すだけです。そのために最初は二分法を考えたんですが、(12)をみてたら縮小写像の方が便利そうだと気づきました。(12)をもう一回Ｃで微分すると、

　　

なのでd²Ｔ/dＣ²は正。よってy＝0～1の範囲でＴ(Ｃ)のグラフは、下に凸な単調減少です。という事はＳ＝Ｔ(Ｃ)－1＋ＣとＳ＝Ｃとの関係は、図-3のようになります。

 

 

　そうであれば例えば初期値Ｃ＝1から始めて、青点線矢印のように進んでいけば、Ｔ(Ｃ)－1＋Ｃ＝ＣすなわちＴ(Ｃ)＝1となるＣに収束します。それも相当に速く。手動でもできそうです(^^)。

 

　Excelでやってみた結果が以下です。

 

 

　単純に台形公式でＴの面積を計算し、結果を右の表にメモリながらＣの値を入れかえるの繰り返しです。ちなみに値コピーしないでＣの値を入れかえようとすると、「循環参照だぞっ！」ってExcelに怒られます(^^;)。

 

　最終結果は5回で収束しました。

打ち切り条件はＴ＝1が目標だったので、Ｔ－1＜10^-6です。いつもは関数の単調性さえわかれば、なんでもかんでも二分法に持ち込むんですけれど、きちんと状況を調べたうえで使用する縮小写像も捨てたもんじゃないな、と思いました。もし二分法でやったら、同等以上の精度に達するのに23回ほどかかります。

 

 

　Ｃ＝0.599942より(4)から、y'(0)＝1.095392です。これで普通に前進差分に持ち込めますが、じつはもう答えが出てます。図-4の計算表の最終段階は（n＝5）、表-2です。

 

表-2の面積計算は微分方程式(1)そのものです。積分はy＝0～1で行ったので、境界条件の半分はみたします。そしてΔＴの列は、yがある値に達するまでの経過時間tであり、y＝1でt＝1になるようにＣを決めたのでした。これはy＝y(t)の関数表じゃないですか！(^^)。

 

（執筆：ddt³さん）