微分方程式を解くアルゴリズム (1/1)


a) 基本

基本は、「連続変数の微分方程式を近似的な差分方程式に変換して解く」です。これはすべてのアルゴリズムに共通して言えることです。 たとえばy'=1という微分方程式を考えます。これはすぐに解けて、y=x+cが解になっています。これを差分にして解く、というのは次のようなことです。まず考えるxの領域は、たとえば0≦x≦1だったとします。この領域をN個に分割して、幅1/Nの線分がN個ある状況を考えます。x=k/N (k=0,1,2,...)のところに線分の端が来るわけです。差分化とは、xの微分を、この線分の両端での差に置き換える操作です。つまり、y'を (y((k+1)/N)-y(k/N))/(1/N)のような差分で置き換えるわけです。これはNを十分大きくとると微分に戻りますが、有限のNだとその「近似」を与えます。いったんこのように近似すると、差分方程式は単なる漸化式なので、初めの数項(今の例だとy(0)の値)が分かっていれば、コンピュータは順番にy(k/N)を求めていってくれます。このように差分化によって漸化式に落とすことで、数値的に方程式を解くことが可能になります。

一般的な状況では、次のような状況になります。。 考えている微分方程式が閉区間 [a,b] で定義されているとして、

などとおいていけば、微分方程式が差分化されることになります。 ポイントは、

  • 出発点が必要。 (初期条件。上の例では
  • 細かいステップで逐次、解を計算していく。
           (ここの計算法の違いでいくつかのアルゴリズムに別れる)
  • 閉区間の終わりにきたら終了。

あたりでしょうか。Mathematica では差分方程式を RSolve で解く事ができますが、線形なものに限ります。数値計算で解きた いような微分方程式には非線形なものが多いため、このままではあま り有効ではありません。

また、差分化の定義は、線分の幅を小さくとった場合に微分に戻るもの、ということなので、 いろんなバリエーションがあります。一般に、一番素朴な差分化より、改良した差分化 を使ったアルゴリズムの方が精度よく解を求められます。以下ではこの改良した差分化の方法を紹介します。

b) Euler法

逐次解いて行く部分に Taylor展開を使うアルゴリズムです。 微分方程式がy'=f(x,y)の形だったとします。Euler法とは、この差分化を 以下のように行う方法のことです。

を採用します。最後の項は最も素朴な差分化では必要ないのですが、これがあることによって差分化による誤差を減らすことができます。しかし、この表式からわかるようにfの微分を計算する必要があり、高速計算には不向きなアルゴリズムです。

c) Runge-Kutta法(2次)

Runge-Kutta では、逐次解いて行く部分に

を使います。ただし、

です。通常は

(微小区間内で、両端で定義される値ではなく、中間での値を採用することに相当。)

および α=β=1 に選びます。したがって

を使って逐次、計算していきます。導関数評価に微分を計算しなくてすむ アルゴリズムになっています。

d) 4次のRunge-Kutta法

最もよく使われるアルゴリズムで、c)の自然な拡張になっています。

を使って逐次計算を行います。ただし

です。実際のRunge-Kutta法では、この他にも高速かつ精度を保つ改善が いろいろ行われています。 (例えば、関数が急激に変化するところではメッシュを細かく取り、 そうでないところでは大きめのメッシュを取って高速化する、など。 Adaptive Runge-Kutta法と呼ばれます。)

e) 多重ステップ法

いままで紹介した方法はすべて、逐次解いていく際に、次のステップの 解を直前のものから求めていました。 直前だけでなく、その前の数ステップ分も考慮して次の解を求める方法を 一般に多重ステップ法といいます。

例えば、4ステップ前のものを考慮した Adams-Bashforth 予測子/修正子法では、 1ステップあたり4回の導関数評価が必要だった4次の Runge-Kutta法 に較べ、 1回で済むために計算が高速化できます。 その反面、これらの4点が十分な精度を持った値でなければ、以降の計算が無意味 なものになってしまうという性格を持っています。 Runge-Kutta法に較べると解の安定性に問題があり、この方法を利用するときには 解くべき微分方程式の全体的な性質をよく把握しておくことが必要になります。

Mathematica で数値的に微分方程式を解く場合、NDSolve という関数が用意されています。 これはいくつかの多重ステップ法を問題に応じて自動的に使い分けています。 また、Runge-Kuttaで解くことも、オプションを指定することで可能になっています。


Home Top Back