CBD：循環参照を使って微分方程式を解く！(HAND MADE MAID)

ここ最近、循環参照ゴリ押し気味ですみませんが、マイブームなので生暖かい目で見守ってやってくれれば＾＾；


まずは、速度の1乗に比例した空気抵抗を受けながら落下する物体(質点)の運動をシミュレートしてみましょう

運動方程式はこうでしたね。

F=ma

F：力、m：質量、a：加速度

加速度aを速度vの時間微分としてみますと
a=dv/dtなので、運動方程式は

F=mdv/dt

になりますね。

では力Fはどうなるかといいますと、今回は速度の1乗に比例した空気抵抗を受けながら重力落下するという条件なので、
速度の向きを鉛直下向きに取って

F=mg-rv

となりますね。g：重力加速度、r：空気抵抗

ですから、運動方程式は結局

mdv/dt=mg-rv

となります。
しかしながら、これではエクセルさんが困ってしまうので、微分を差分に翻訳してあげましょう。
dv→Δvなので、v2-v1
dt→Δt
つまりdv/dt→(v2-v1)/dtとなるわけです。

これを運動方程式に当てはめますと

m(v2-v1)/dt=mg-rv2

になりますね。

これをv2についての式にまとめなおしますと

v2=(v1+gdt)/(1+rdt/m)

になります。

これを、図の速度vのところに入れてやります。

v=IF(sw,(v1+gdt)/(1+rdt/m),v0)

sw：初期化スイッチが0(オフ)のときはvの初期値v0を維持し
1(オン)になったら循環参照を開始して

vのセルに(v+gdt)/(1+rdt/m)

を代入する、という意味です。ここでvは自分自身のセルなので循環参照しています。
循環参照を開始する際に、ツール→オプション→計算方法→反復計算にチェックを入れ
反復計算回数を1にしておきましょう。変化の最大値はデフォルトの0.001で構いません。

dt=0.5、m=100/g(mg=100)、r=1、初期速度v0=0としておき
反復計算開始スイッチswを一旦0にしてから1にして、delでも押しっぱなしで再計算し続けると

速度vが0から終端速度vlast=mg/r=100に漸近していくのがわかりますね、わかります＾＾




＝＝＝＝＝＝＝＝＝
では次に、2階の微分方程式に挑戦してみましょう。
力が距離に比例する振り子を例に取ります。

運動方程式はmdv/dt=Fでしたね。
今回このFは、距離xに比例してxとは逆向きに作用する力ですので
F=-kxと置くことにします。

ここで、距離という変数が出てきたので、速度と距離の関係を使って距離を定義すると
速度は距離の時間微分なので

dx/dt=v

でしたね。

この2本の連立微分方程式

dx/dt=v
mdv/dt=-kx

を差分の形で解くことが
2階の微分方程式d^2x/dt^2=-kxを数値解析で解くことに相当します。

dx/dtも先ほどと同様に微分から差分の形にすると

(x2-x1)/dt=v2

の形になるので、x2の式にすると

x2=v2dt+x1

になります。

運動方程式F=maのほうも差分化しますと


dv/dt=-kx2/m

が

m(v2-v1)/dt=-kx2

になるので、v2の式にすると

v2=-kx2dt/m+v1


になります。

下の図の時間刻み幅dt=0.1、質量m=1、比例係数k=1、初期速度v0=1、初期変位x0=1
として、

変位xと速度vにそれぞれ

x=if(sw、vdt+x、初期値x0)
v=if(sw、-kxdt/m+v、初期値v0)

を意味する循環参照を施してやってdelボタン押しっぱなしとかで再計算し続けますと、
速度と変位がいつまでも往復する様子が描けると思います。


ちなみに、k=m=x0=v0=1にそろえて
v-x(速度：縦軸、変位：横軸)をグラフにして、
先ほどの空気抵抗のとき同様、循環参照初期化スイッチ(sw)を一旦オフ(0)ってからオン(1)にしますと

半径約実効値1.4(√(1^2+1^2)=√2)の円の上を延々と点が回り続けるのがわかりますね、わかります＾＾





宇宙ｷﾀｰｰｰ！！

にほんブログ村