[2階常微分方程式(線形系)] [数値解析]
[2階常微分方程式(線形系)]
保存系と並んでなんとかなるのが、線形系です。線形系なら非保存であってもなんとかなります。いっぽう保存系なら、非線形であってもなんとかなります。悩ましいなぁ~(^^;)。
この前は線形系かつ定数係数という場合をネコ先生にアップして頂きましたが、線形系なら非定数係数であってもなんとかなるはずです。やる事は定数係数の時と、基本的に同じだからです。
ただし自分はシンプレクティック法大好きなので、ここではシンプレクティック法の精神に則った準シンプレクティック法を使います。参考URLを最後にあげます。
ネコ先生も書いておられるように、この問題(1)は、1階の連立微分方程式にひき直せます。y'=dy/dt=pとして、
です。
ここでシンプレクティック法の精神を発揮して、等速運動させてから等加速度運動させる、に注意しながら(3)を前進差分に持ち込むと以下になります。τはtのステップ幅です。
です。
なので(5)は、
となります。
(6)最右辺の係数行列を、とすれば、
になります。特にj=0とすれば、
ここでy(0)=y0,p(0)=p0とすれば、yn=y(nτ),pn=p(nτ)です。よってt=0からnτ後のyn,pnを知りたいのであれば、(9)のA1A2・・・Anを計算すれば良い事になります。各Ajは(7)です。そんな事、数値的にはなんぼでもできます。以下はExcelの画面のコピーです。意味は余りにも明らかだと思うので、特に説明しません(^^)。
Excelで計算した結果、n=20,τ=0.05,nτ=1の時、
・・・に、なりました。従って(9)より、
なので、境界条件(2)を代入すると、
です。整理すると、
となります。これより、p0=1.391881057となりました。
y(0)=y0=1,p(0)=p0=1.391881057と初期値が決まったので、(4)を使って「準シンプレクティック、行きまぁ~す!」(^^)。
まず(4)の下段をpj+1について解き、(4)上段と連立させます。
これで準備完了です。y0=0,p0=1.391881057として(10)を使い、(y0,p0)→(y1,p1)→・・・→(y20,p20)と追跡して行けば、Excelにより(^^)図-1です。
図-1
理論解は、もちろんネコ先生のを使わせてもらいました。理論解は、どうやって出したんですか?(^^;)。この質問をしたくて、この記事を書いたようなもんです(※)。
ちなみに準シンプレクティック法は正式な用語でもなんでもなく、また確立された方法でもありません。でも、けっこういける気がしてます(^^)。
[参考URL]
http://www.sunagonet.co.jp/wp/wp-content/uploads/2018/01/c774c6cbdbfc66af0583295744d978fa.pdf
(執筆:ddt¹さん)
(※) 第15回 2階線形非同次微分方程式の解法2
http://nekodamashi-math.blog.so-net.ne.jp/2017-08-24
の問題3の(1)に出ております(^^)
2階線形(非同次)微分方程式の標準化(問題2)という計算のテクニックを使って求めています。まぁ、定数変化法の一種ですね。
コメント 0