時間発展問題の数値計算

単振動問題

○単振動

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

  d2 Z
m ------ = F
  d t2
のような2階微分のはいった形をしている。

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

  d2 Z       k
  ----- = - --- Z
  d t2       m
となり、ここで ω2 = ( k / 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つの定数を決めるために、 初期条件として、 t = 0 で、位置 Z だけでなく、速度 d Z / d t も 与えてやる必要がある。

○単振動の差分解

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

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

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

  d W
  --- = - ω2 Z                                         (2b)
  d t
こうすると、見かけ上2階微分はないので、1階差分で表現できる。

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

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

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

      tzero  =  0.0d0
      wzero  =  0.0d0
であることを前提としており、これらを変えると正しくなくなる。 (もっと一般的に有効な表現にするにはどうしたらいいだろうか?)

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

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

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

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

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

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

減衰問題(指数関数)振動問題(三角関数)
前方差分____
後方差分____
[その他、ためしてみたもの]____

1996-06-29; 2000-06-28; 2002-07-11
増田 耕一 (MASUDA Kooiti)
[目次へ]