第二回-02 運動方程式を数値的に解いてみよう
運動方程式の例としておもりとばねによる単振動系を考え、シミュレーションを行ってみよう。
第一回のページでも紹介した、摩擦のない床におもりとばねが接続された系のシミュレーションをあらためて実行してみよう。
おもりと異なる側のばねの端は壁に固定されているものとする。
前回と同じファイルであるが、あらためてダウンロードしてみよう。ダウンロードの際に前回とは異なる名前を付けて保存しても良い。
ここで、「開発」タブ→「Visual Basic」ボタンをクリックし、プログラムを見てみよう。
「開発」タブがないという学生は、第一回のページをもう一度学びなおすこと。
プログラムの冒頭は以下のようになっている。
今回はここにこのシステムの運動方程式が記されていることに着目しよう。
前ページで運動方程式
は、x の微分 (速度) を表す新たな変数 y を導入することで、以下のように書き直すことができることを学んだ。
いま考えているシステムの運動方程式は
mx'' = - k x
なのであるから、これは以下のように書きなおすことができる。
x' = y
y' = - k x / m
この式が、Visual Basic 上では
dxdt(0) = x(1)
dxdt(1) = -k * x(0) / m
と書かれていることに着目して欲しい。
プログラム上の変数が以下の文字と対応していることが想像できるだろうか。
x(0) : x
x(1) : y
dxdt(0) : x'
dxdt(1) : y'
なお、プログラミングに興味のある学生向けに書いておけば、「(0)」、「(1)」 は Visual Basic における配列の添字である。
前回同様、デフォルトの状態でシミュレーションを行ってみよう。
「データ作成」ボタンを押すのだった。
このとき、プログラム中におもりの質量を定義する以下の行があることからわかるように、
おもりの質量を 1 [kg] としたときのシミュレーションが実行される。
Const m As Double = 1 ' 質量
今回は、おもりの質量を変えたときのシミュレーションも行ってみよう。
まずは準備として、おもりの質量が 1 [kg] のときのシミュレーション結果が消されないよう、
コピーして別の列に貼り付けよう。
下図では、A列とB列をコピーしてH列、I列に貼り付けている。
C列の「データ y」は今回は用いないのでコピーしなかった。
なお、貼り付けた後、分かりやすいようにI列の先頭を「m=1」と書き換えておこう。
次に、質量 m = 2 [kg] におけるシミュレーションを行ってみよう。
Visual Basic 上で以下の部分を「2」に書き換えれば良い。
そして、再び「データ作成」ボタンをクリックすると、m=2 に対するシミュレーション結果が得られる。
そこで下図のように、B 列のデータを先程の右隣、すなわち J 列に貼り付けよう。
さらに先程と同様に J列の先頭を「m=2」と書き換えておこう。
以上の準備のもと、H列～J列のデータをグラフ化してみよう。A列～C列ではないので注意すること。
前回同様、散布図 (直線) でグラフ化すると、下図のようになる。
m=2 の単振動が、m=1 のそれよりも遅くなっていることがわかるだろう。
この結果を元に課題を提出してもらう。
←動力学と静力学／解析解と数値解／課題→
Excel VBA で微分方程式をシミュレーションしように戻る