. 運動方程式の数値解法.. ニュートン方程式の近似速度は, 位置座標 の時間微分で, d と定義されます. これを成分で書くと, d d li li とかけます. 本来は が の極限をとらなければいけませんが, 有限の小さな値とすると 秒後の位置座標は速度を用いて, と近似できます. 同様にして, 加速度は, 速度 の時間微分で, d と定義されます. これを成分で書くと, d d li li とかけます. 本来は が の極限をとらなければいけませんが, 有限の小さな値とすると上の式から 秒後の速度は加速度を用いて, 4 と近似できます. ニュートン方程式は, で与えられます. 次元の場合に成分で書くと, 5 となります. 式 5 と近似式,4 を用いると物体の運動は,
6 を計算してやれば, 求まることが分かります... オイラー法時刻を 刻みに分割して とすると,6 式は, 7 と書き直すことができます. これにより物体の運動をシミュレートすることが可能となります. 微分方程式をこのように, g d g g 8 と近似して解く方法をオイラー法と言います. 以下にオイラー法のスキームを図示します. 時刻位置速度加速度力.. 物体の投射をオイラー法で解くここでは, オイラー法を使って斜方投射を解析するエクセルファイルを作成します. 斜方投射のニュートン方程式は, g で与えられ, 初速度を, 水平面からの投射角度を とすると初期条件は, si, cos, で与えられる.
... 初期値と定数の設定 まず, 右図のように初期値を入力します. 英数字は 全て半角英数で入力してください. それぞれの意味は, 時間刻みΔ 物体の質量 = での 座標 = での 座標 初速度の大きさ he 物体を投げ上げる時の水平線と の角度 初速度の 方向成分 図 初速度の 方向成分 です. 次ぎに入力した値に名前をつけます.B セルを選択して, プルダウン メニューの [ 挿 入 ] [ 名前 ] [ 定義 ] とクリックします. 図 のダイアログが開くので [OK] をクリックしま す. これで,. に という名前が付きます. 以下, 同様にして,B4 セルから B8 セル までそれぞれ,,,,, he と名前をつけていきます. 図 B9 セル,B セルにそれぞれ以下のように入力します. B9 =*coshe*pi8 B =*sihe*pi8 4 と同様にして,B9 セル,B セルにそれぞれ,, と名前をつけます.
... 計算式の入力 図 のように A セルから J セルまで文字を入力します. 図 A セルから A6 セルまで,,,,..., 5 と数字を入力します. そのために,A と A4 にそれぞれ, とを入力して,A と A4 の両方を選択し, 選択された枠の右下の をドラッグして A6 セルまで持っていくと簡単に入力できます. 図 4 によって, 時刻を計算します. そのために,B セルに =A* と入力します.の時と同様に B セルが選択された状態で, 枠の右下の をドラッグして B6 セルまで持っていくと時刻が入力されます. 4 を入力します. 軸方向には力は働かないので,I セルに を入力します.の時と同様に I セルが選択された状態で, 枠の右下の をドラッグして I6 セルまで持っていくと がコピーされます. 5 を入力します. g なので,J セルに =-**9.8 を入力します.の時 と同様に J セルが選択された状態で, 枠の右下の をドラッグして J6 セルまで持っていくとこの式がコピーされます. 6 を入力します. なので,G セルに =I と入力します.の時と同様に G セルが選択された状態で, 枠の右下の をドラッグして G6 セルまで持っていくとこの式がコピーされます. 7 を入力します. なので,H セルに =J と入力します.の時 と同様に H セルが選択された状態で, 枠の右下の をドラッグして H6 セルまで持っていくとこの式がコピーされます. 8 初期値を入力します. 各セルに以下のように入力します. C = D = E = = 9 いよいよ計算式を入力します. 各セルに以下のように入力します. 行目は入れた数式 の意味を示します. 4
C4 =C+E* * D4 =D+* * E4 =E+G* * 4 =+H* * 9 でセル C4,D4,E4,4 に入力した数式を の要領でセル C6,D6,E6, 6 までコピーします. ここまでで表は, 図 5 の様になっているはずです. 以上で, 数値的に差分方程式の解が得 られます. 図 5... 結果のグラフ化ここでは結果をグラフにします. セル C-C6 と D-D6 を選択します. プルダウンメニューの [ 挿入 ] [ グラフ ] とクリックします. 5
図 6 図 6のグラフウィザードダイアログが開くので,[ 散布図 ] を選んで,[ 次へ ] をクリックします. その後 回ダイアログが開きますが, とりあえず [ 次へ ] をクリックし最後のダイアログで [ 完了 ] をクリックします. 4 すると図 7のように自由落下の様子がグラフになります. -..4.6.8 - - 系列 -4-5 -6 図 7 6
..4. 異なる初期条件の計算 と he を変化さ せると, 異なる初期条件の時の解が得られます. 例えば図 8は,=5,he=45 の時の解です..5 -.5 4 系列 - -.5 図 8..5. 理論値との比較 斜方投射の問題の正確な答えを我々は知っています. 時刻 = に, から初速度, で物体を投げたときの時刻 での物体の座標は, g で与えられます. そこで, 数値解析による値を正確な理論値と比較して見ます. 各セルに以下のように入力します. L M L M _ec _ec =+*B =+*B-.5*9.8*B^ セル L を L6 まで,M を M6 までそれぞれコピーをする. こ れで正確な理論値が計算されます. 図 8の数値解と比較するために一緒にプロットします. そのために, 図 8のグラフを選択して, プルダウンメニューの [ グラフ ] [ データの追図 8 7
加 ] とクリックします. 図 9の [ データの追加 ] ダイアログが開くので, セル L-L6,M-M6 を選択し, [OK] をクリックします. 4 図 9の [ 形式を選択して貼り付け ] ダイアログが開くので,[ 先頭行を項目列として使用する ] にチェックを入れて,[OK] をクリックします. 図 9 5 図 の様に理論値と数値解析 値が同時にプロットされます..5 -.5 4 系列 系列 - -.5 図 問題 上の斜方投射の解析を, 速度に比例した空気抵抗がある場合に拡張しなさい. ニ ュートン方程式は, g となる. を適当に変えて結果をプロットしなさい. 問題 上の例を参考にして, 単振動を解析するエクセルファイルを作りなさい. 初期条件として,,,,he= とします. 運動方程式は, k ですから, 変わるのは初期条件と力の部分だけです. ただし,k= とします. 変位 を時刻 の関数としてプロットしなさい. また, 理論値と比較しなさい. * 間違いの元になるので, 前の課題と違う新しいエクセルファイルを作ること. 前の課題のファイルをコピーして使うと便利である. 8
.4. 蛙飛び法 問題 を解くと右図のような結果が得 られます. 単振動では,- から の間を 振動するはずなので明らかに間違って います. これはオイラー法による数値解 析の誤差のためです. オイラー法は簡単ですが, 誤差の大きい 方法です. ここでは, もう少し誤差を小 さくする方法を考えましょう. 誤差を小 さくするためには, 時間刻み を小さくする より良いアルゴリズムを考える.5.5.5 -.5 5 5 - -.5 - -.5 - という 通りの方法が考えられます. ここでは, 計算量を変えずにほんの少し式を変更す るだけで劇的に誤差を小さくできる蛙飛び法と呼ばれる方法を紹介します. オイラー法の 近似が良くない理由は, 時刻 での微分を と を結ぶ線分の傾きで近似していること にあります 下図. これを, と を結ぶ線分の傾きで近似すればかなり 良くなることが期待できます 下図 b. b 具体的には, 加速度と位置座標を評価する時間を速度を評価する時間と だけずらして, 9
9 とします. 以下に蛙飛び法のスキームを図示します. 時刻位置速度加速度力 蛙飛び法がオイラー法よりも精度の高い計算法であることを見ておこう. オイラー法はテイラー展開, ''' 6 '' ' で 次まで採った近似であり, 刻み幅 の 乗の誤差を持つ. 一方, 蛙飛び法は, ''' 6 '' ' ''' 6 '' ' の片々を引いた式から得られる近似式, ''' 4 で 次まで採った近似であり, 誤差は刻み幅 の 乗になる. 例えば. とすると 倍精度の良い計算となるのである. 言い方を変えると, 同じ精度の結果を得るのに 分の の時間で計算できることになる.
.5. 物体の投射を蛙飛び法で解く. で作成したエクセルファイルを編集して 蛙飛び法 のエクセルファイルを作ります. そのため, 前のファイルをコピーして別のファイルを作ります. 蛙飛び法では, 座標と速度の間に だけ時間差があります. そのため, 各セルに以下の値が入るように変更します. + + と は仕方ないのでオイラー法で,.5.5 により計算することにします.
.5.. エクセルファイルの編集 上記を実現するために, 各セルを以下のように変更します. 列目は数式の意味を示します. E =+.5*G*.5 E4 =E+G4* =+.5*H*.5 4 =+H4* セル E4 を E6 までコピーします. また, セル 4 を 6 までコピーします. 今度は右図のように, 理論値と数値解析値がほぼ一致しています. をもっと大きくしてもよく一致していることがわかります..5 -.5 - -.5 系列 系列 4 問題 問題 の単振動を蛙飛び法を用いて解きなさい. 問題 4 万有引力を受けた物体の運動を解きなさい. 質量 と M の物体に働く万有引力の M 大きさはG です. 質量 M の物体は原点にあり,M は に比べて十分大いため質量 M r の物体は動かないと仮定する. このとき質量 の物体に働く万有引力は原点に向きます. 従って, 質量 の物体に働く万有引力を成分で書くと, GM GM となります. 簡単のため, GM, とし, 初期条件,,,,he= としてエクセルで解きないさい. 問題 5 問題 4で初期条件,,,he を色々と変えて計算してみなさい.