このページを正しく読むためには MathML の表示のための数式フォントが必要です． Firefox, Netscape(6.0以降)など Moziila 系ブラウザであれば， 次の手順に従ってフォントをインストールしてみてください． → Firefox のインストールとその他の設定の手順
なお，IE の場合には MathPlayer というのも組み込まないとハングするようです．

落下運動のシミュレーション

Euler 法で運動方程式を解く

力は加速度を生み出す

質量が $m$ の物体に， 大きさ$ F $ の力 が働くと， その物体は 加速度運動を起こします．その加速度$ a $ は次の式で 表されます．

(1a) $ a = \frac{F}{m} $
あるいは次のようにも表されます．むしろこちらのほうが有名です．
(1b) $ F = ma $
これはニュートンの第2法則と呼ばれる最も重要な力学の法則です． 力学的な観点からは，力というのは物体に加速度を与えるものです．

ちなみに $ F $, $ m $, $ a $は それぞれ force, mass, acceleration の頭文字から採られています． acceleration はアクセルという言葉のもとですね．

次に加速度は，速度が時間的に増加減少する割合ですから， 速度$ v $ を時間$ t $ で微分したものとして 表されます．

(2) $a = \frac{dv}{dt} $
さらにいうと，速度は位置$ x $ を時間で微分したものです．
(3) $v = \frac{dx}{dt} $
どんな物体の運動も 以上の関係式から導き出すことができ， いつ，どこから，どんな速度でその物体が動き始めたかが分れば， その後の位置や速度は完全に予測できます．必要なことは 力 $ F $ がどのような性質を持っているかを知ることだけです．

以上からお分かりのように，物体の運動を解くという問題は， 物体に働く力がどんな性質を持っているかという点にポイントがあります． 力には日常的な場での重力，万有引力（大きなスケールで見た重力）， 電気的な力，摩擦力などなど，さまざまのものがあります．

ここまでの説明では，物体は大きさを持っていないことを前提にして います．大きさをもつ物体は回転したり変形したりできるので， その取り扱いはかなり複雑になります．しかし，原理的な部分では， 上のような運動方程式が基になっています．
さらに流動する気体や液体の 運動を扱うことになると，普通の力学と原理はほとんど共通なのですが， 数学的にかなり難しい取り扱いを要求されます．気象のシミュレーションが 大掛かりになるのは，そのような難しさを持っているからです．

2つの微分方程式

実際に物体の運動を解くには，まず式(1a)から加速度 $ a $ を求めます．質量はたいていの場合不変ですから，力の性質だけが 知られていれば，そこから加速度が求められます．もっとも，実際に値を出すの ではなく，式の形を求めるのです． 式(3)はいつでも変りません． そうしておいて， (2)と(3)の2つの微分方程式を連立させて解くのです． このやり方さえ分っていれば，ロケットでも太陽系でも，その運動を 予測することが出来るのです．運動のシミュレーションというと難しそうですが， 実は簡単なものなのですね．

Euler 法で運動方程式を解く

Euler 法で式(2),(3)の微分方程式を解くことを考えましょう． 加速度 $ a $ は，力の性質によって変化しますから， 位置と時間の関数になり，$ a(x,t) $ のように表されます． たとえば地球から遠い人工衛星は弱い重力で引っ張られますし， 綱引きのように引っ張る力が時々刻々変ることもあるというわけです．

なお 場合によっては，力が速度の関数にもなることもありますが，それは 必要なときに考えることにしましょう．これはたとえば空気抵抗のように， 速く運動するほど大きな抵抗力を受けるようなケースに当てはまります．

式(2),(3)を Euler 法の 形で書くと，それぞれ次の(4),(5)のようになります．

(4) $ v_1 = v_0 + a(x_0,t_0) \delta t $
(5) $ x_1 = x_0 + v_0 \delta t $

これらを解く手順は次の通りです．

  1. $ a(x,t) $ はすでに知られている．
  2. 問題の設定から， 初期条件として $ x_0 $, $ v_0$ を与える．
  3. 式(4)を使って $ v_1 $ を求める．
  4. 式(5)を使って $ x_1 $ を求める．
  5. $v_1 \leftarrow v_0 $ , $x_1 \leftarrow x_0 $ の代入を行って，上の手続きを必要な回数繰り返す．

計算のための設定

それでは物体を地上から投げ上げて，それが落下してくる運動を Euler 法で解いてみましょう． 必要な設定は次の通りです．

ところで，上の設定では質量$ m $ が登場しますが， 実際には，質量は関係なくなってしまいます．なぜかというと， 物体に働く力（重力）は$ -m g $ で，一方物体の質量は$ m $ なので， 式(1a)から $ a = -g $ となるからです． さらに $ a $ は定数で，位置にも時間にも依存しません． したがって，式(4),(5)はさらに単純化されることになります．

実際の計算のためのパラメータとして，ここでは次のように定数を決めてみます． 

a = -9.8
vstart = 9.8
tstart = 0.0
xstart = 0.0
tend = 2.2
dt = 0.01

解析的な解

結果を吟味するために，解析的な厳密解も用意しましょう． 物体を上方に投げ上げる問題は，高校の物理の基本問題として 取り上げられ，その解は次の式で表されます．

(6) $ v = v_0 + a t $
(7) $ x = x_0 + v_0 t + \frac{1}{2} a t^2 $
GNUPLOT を使えば，これらは即座にプロットできます．

結果の比較

下に得られた結果を示します． 赤，緑はいずれも Euler 法によるもので，dt をそれぞれ 0.2 と 0.01 にとっ ています．青の解析解と比べると，それなりに似た挙動はとっているものの， やはり「ちゃんと曲がりきれない」という Euler 法の近似の限界があることがよく 分ります．

運動方程式を2次のLeap-Frog 法で解く

それでは，運動方程式

(2) $a = \frac{dv}{dt} $
(3) $v = \frac{dx}{dt} $
を Leap-Frog 法で解くことを考えましょう． ここで留意しなければならないのは，式(2)はニュートンの第2法則 $ F = ma $ に 直接結びついた ものであって， (3)はそれを受けて位置を計算するためのものに過ぎないということです． その観点から考えると，まずは速度 $ v$ が 加速度 $ a $ によって，時間とともに どう変っていくかという点を式(2)で追うのが筋です．

そしてその結果から式(3)の左辺を得ると， 今度は位置 $ x$ を計算することができ，もしも 力が位置に依存するのであれば（しばしばそうなります）， 結果的に $ a $ は $ x $ に依存することになって，式(2) の計算に組み込まれます。

それでは式(2)(3)を Leap-Frog 法の精神で差分方程式にすると どうなるでしょうか。 とりあえずその結果を示します。

(4) $v_1 = v_0 + a_0 \delta t $ 現在の位置と加速度から半歩先の速度を計算する
(5) $x_2 = x_0 + 2 v_1 \delta t $ 現在の位置と半歩先の速度から1歩先の位置を計算する
(6) $ v_2 = v_1 + a_1 \delta t$ 半歩先の位置と加速度から1歩先の速度を計算する

ここで，$x_0,v_0$などは「現在」の値， $x_1,v_1$などは「半歩先」の値， $x_2,v_2$などは「1歩先」の値と考えます．

Leap-Frog法で解く落下運動

それでは前に使った設定の下に,落下運動を Leap-Frog 法で解いてみることに しましょう． 実際には，半歩の幅を dt, 1歩の幅を dt2 ( =dt*2 : ループの 中で dt * 2 という計算を行わせるのは実行速度を遅くする ことになるので，あらかじめループの外でこのような値をもつ 変数を用意しています) とし， Euler 法のプログラムのループの中を書き換えるだけです． 得られた結果を下に示します．Euler 法だったら大きな誤差が出る ような時間刻みの設定を行ってみても， Leap-Frog 法は解析的な解によく一致した結果を与えて いることが分かります．これはたまたま 解が2次式で与えられるケースなので，厳密な解と Leap-Frog 法による解が一致する幸運に恵まれたのであって， いつでもこんなふうにぴったりというわけではありません． が，単純な Euler 法に比べて Leap-Frog 法が格段に よい精度の解を与えてくれることは確かです．

Java アプレットで運動方程式を解く

Java のプログラムを作成して実行するには，コンパイラのインストールが必要です。 インストールの方法はこちらです。

Euler 法で解いてみる

次のソースは，Euler 法でボールの運動を解くためのものです。 この段階では力は働いていないので，等速直線運動を見せるだけです。 ここからスタートしましょう。 次の3つのファイルをダウンロードしてください。

これらを利用するには次のようにします。なお行頭に '>' とあるのは， シェル画面からのコマンド入力を意味しています。

  1. 上記の3つのファイルを置いたディレクトリをカレントディレクトリとする。
  2. ソースをコンパイルして 拡張子が .classのファイルを作る。
    3. > javac Ballmotion01.java
  3. アプレットを走らせるためのHTML ファイルを appletviewer で読み込む。 
    4. > appletviewer ballmotion.html

    上のやり方の代わりに次のようにしても，このソースファイルの場合には 走るようになっています(ただしこれは裏技です)。
    > appletviewer Ballmotion01.java
    上の作業を繰り返して，プログラムを改善してください。
  4. 完成したらブラウザから上記の HTML ファイルを読み込むことで， アプレットが実行できます。

バネや振り子の振動をシミュレートする

バネの振動をシミュレートするには，次のフックの法則を使います。

$ F = -k x $
ここで $F, k , x$ はそれぞれ，力，バネ定数，伸びです。力の定義を これに変えるだけで，振動を起こさせることが出来ます。やってみてください。 ファイル名は SpringMotion01.javaとしてください。また バネ定数は kS とし，とりあえず 1.0 という値を与えてください。 このさい，力は $x$ 成分だけなので，いくらか単調なシミュレーションになります。

振り子も近似的にバネと同じような振動を見せます。 錘に働く重力と糸による運動方向の制限を考えることで， 位置と力の関係が得られますので，その関係を使って 運動方程式を解けばいいのです。

近似を改良する

これまでのところでは，運動がいまいちへんでした。これは近似がよくないためです。 今度は Leap-Frog 法を使って，より正確なものにしてみましょう。

惑星の運動をシミュレートする

次のソースは上のプログラムから作ったもので，太陽の周りを回る惑星の運動をシミュレーションで 見せるものです。ただ，Euler 法を使っているために， 一定の軌道を回らずにだんだんずれていってしまいます。これを Leap-Frog 法にしたら，改善されるかどうかやってみましょう。 改良版には， PlanetMotion_LF01.java と名前を付けてください。

PlanetMotion01.java

プログラムをバージョンアップする上での注意

Java のソースでは，ソースファイル名と，その中で与えられている 最外部のクラス名が一致していなければなりません。つまり， 上記のソースファイルの名前は 

Ballmotion01.java
となっているので，最も上のクラスの定義は次のようになっています。 
public class BallMotion01 extends Applet
  implements Runnable, ActionListener{
ついでにいうと，このクラス名はソースの他のところにも 現れることがしばしばあります。 実際，このファイルにも該当する箇所があるの ですべてを見つけて変更してください（エディタの機能を活用すること！）。

バージョンアップのためにソースファイル名を変更するときには，ファイルの 中のこのような箇所についても変更が必要です。間違えると わけの分からないことが起きたりしますので，注意してください。

ソースの PDF ファイル

印刷用に LaTeX で作って変換したものです。 → javasource.pdf

小波秀雄, 2006/06/07