時間発展問題の数値計算

単振動問題

○単振動

運動方程式は、質量 m の物体の位置座標を Z 、力を F とすれば、

  d2 Z
m ------ = F
  d t2

ここで、 F が基準位置 (簡単のため Z = 0 とする)からの変位 Z に比例する 復元力 -k Z (k > 0)であるような場合を考えよう。 たとえば、ばね定数 k のばね振り子を想定すればよい。すると、

  d2 Z       k
  ----- = - --- Z
  d t2       m

  d2 Z
  ------ = - ω2 Z                                     (1)
  d t2

これの解析解は、

Z = Z1 exp ( i ω t ) + Z2 exp ( -i ω t )

Z = Z3 cos ( ω t ) + Z4 sin ( ω t )

Z = Z0 cos ( ω t + φ0 )

○単振動の差分解

2階微分方程式を差分化するには、 2階差分(差分の差分)をとる方法も考えられるが、実際にプログラムの 形にする方法は簡単ではなさそうである(わたしには説明できない)。 ここでは、連立1階微分方程式に変形する方法を説明しよう。

速度 W も変数と考えれば、(1)式は次の2つの式に分けて書くことができる。

  d Z
  --- = W                                                (2a)
  d t

  d W
  --- = - ω2 Z                                         (2b)
  d t

減衰問題と同様に、前方差分で表現したプログラム例を、 forward.p, forward.f, forward.c として用意した。 ただし、その中では、- ω² を a とおいている。

注意: ディレクトリは別だが、ファイル名が減衰問題のときと同じである。 前と同じディレクトリにコピーする場合は、名前を変えないと 上書きするので、前の情報が失われる。

なお、この例文の「真の解」 ztrue を計算する部分は、

      tzero  =  0.0d0
      wzero  =  0.0d0

○課題 (自由課題を含む)

まず、単振動問題の前方差分のプログラム例を動かし、結果をグラフにしてみる。 deltat, nt(ステップ数)などを変更したほうがよいかもしれない。

減衰問題の例で示したように、誤差がΔtに応じてどう変わるか (理論どおりにΔtの1乗に比例するか)を調べてみるとよい。ただし、 計算の最終時刻でたまたま誤差が小さくなる可能性もあるので、 計算結果から誤差の代表的値を取り出すフィルターをくふうしてみよう。

同じ問題について、差分を、後方差分 (陰公式なので、式変形が必要) または、前方予測後方修正型の差分 (陽公式)のいずれかに変え、 前方差分の場合と、誤差の出かたを比較する。

余裕があれば、その他の差分法もためしてみてほしい。

それぞれの差分法の特徴(誤差の出かた)を、わかった範囲で、次の表のように まとめてみよう。