高階微分方程式の数値解法 - Sλの日記 --うつろな日々--

一般的に微分方程式っていうのはモデルから出発した場合かなり荒い形をしていて、それが計算機で直接扱えるようなものではなかったりします。が、それをどうにか頑張って、近似したりして線形化したりして、Strum-Liouvilleっぽい形にすることが出来ます。変数返還としまくれば。
そうすると、結局4階の方程式は、
   $a_0y^{(4)}+a_1y_{xxx}+a_2y_{xx}+a_3y_x+a_4y=0$
みたいな形になります。非同次でもまあ良いんですけどね。これがn解になっても理論的には数値的に解くことが可能です。
で、一方で、微分方程式の差分解法は、方程式が、
   $y_x=f(x,y)$
と与えられていた場合、
   $y_{n+1}=y_n+f(x,y)dx+O(dx^2)$
但し、dxは有限みたいになります。
で、これだと一次精度なので、これを四次精度まで上げると、
   $k_1=f(x,y)dx, k_2=f(x+dx/2, y+k1/2)dx, k_3=f(x+dx, y+k_2/2)dx, k_4=f(x+dx,y+k_3)dx$
   $y_{n+1}=y_n+1/6(k_1+2k_2+2k_3+k4)$
になります。これがRunge-Kutta法ですね。
っつーことで、どういうことを言いたいかと言うと、数値的に解くには、前処理として、方程式をy'=f(x,y)みたいな奇麗な形にしておかないといけないってことですね。っつーことで、高階の方程式もあたかも一階の方程式として扱えるように、行列方程式にします。
一番最初の方程式を、それぞれの解の導関数を元とするベクトルに、
   $y_i=(y, y_x, y_{xx}, y_{xxx}, y_{xxxx})^t$
みたいにして、この行列を使って上の方程式を大人しい行列方程式、
   $y'_i=b_{ij}y_j$
になります。まあここでEinsteinの総和規約を使ってるんですがね。まあそれはスルーする感じで。
っつーことで、高階方程式を解く場合は上の行列方程式を実装して解く訳です。そして一次精度のEuler法なんて使い物にならないんで、Runge-Kuttaを使う訳ですが、そうすると、Runge-Kuttaのうちf(x,y)の部分は行列方程式のうち、b_{ijy_j}になって、それに基づいてRunge-Kuttaのk₁だとか、k₂を作って解いてくことになるのですね。
これには昨日の夜家に帰るときに歩きながら思いつきましたよ。やっぱ学校帰りの自転車とか歩いてる時間ってすげーなって思いました。大体すげー着想はこういうときに思いつくことが多いんですよね。