配列を使わずに常微分方程式を解く方法

最終更新日: 2000-11-17


常微分方程式の数値計算で、 データをすべて配列に保存し、 あとでまとめて出力するという方法を用いることがある。 プログラムは、たとえば次のようになる:


  反復 {
    p[j+1] = (A*p[j] + B*q[j])*dt
    q[j+1] = (C*p[j] + D*q[j])*dt
    t[j+1] = t[j] + dt
  }

  反復 {
    出力(t[j],p[j],q[j])
  }

しかし、ちょっと大がかりな計算をしようとすると、 これではすぐメモリが足りなくなってしまう。 配列なしで計算ができるなら、それに越したことはない。 それには、たとえば次のようにする:

/* 配列なしのプログラム: 例[1] */

  反復 {
    出力(t,p,q)
    pnew = (A*p + B*q)*dt
    qnew = (C*p + D*q)*dt
    p = pnew
    q = qnew
    t += dt
  }

最後にでてくる  t += dt  というのは 「tdtだけ増やす」という意味である。 ただし、FORTRANでは += という書き方がサポートされていないので
  t = t + dt
と書く必要がある(「tt+dtを代入せよ」という意味)。

例[1]の代わりに、次のように書いてもよい:

/* 配列なしのプログラム: 例[2] */

  反復 {
    出力(t,p,q)
    dp = (A*p + B*q)*dt
    dq = (C*p + D*q)*dt
    p += dp
    q += dq
    t += dt
  }

例[1]と例[2]のどちらでも好きなほうを用いてよい。

なお、上記の数値計算例では、 簡単化のためにEuler法によるプログラム例を示した。 実際には、 Runge-Kutta法などの 高精度の方法を用いるべきである。 課題として、2次あるいは4次のRunge-Kutta法のプログラムを、 配列を使わずに作成してみてほしい。


[このページの冒頭に戻る] [常微分方程式] [微分方程式]