単振り子のシミュレーション

繰り返しの回数:n (n×Δt 秒後に終了)
時間の刻み幅:Δt [s] (n×Δt 秒後に終了)
振り子のひもの長さ:l [m] 参考). 重力加速度g=9.8[m/s]
角度の初期条件:θ0 [deg] 孤度法([rad / ラジアン])ではなく、
度数法([deg / °])で入力すること。
参考). 角速度の初期条件:ω0=0[rad/s]



解説

単振り子の方程式は、二階微分方程式:
  d2θ  
l
gsinθ=0
  dt2  
で表されるが、
 

ω
dt  
と置くと、
  g  

=-
sinθ
dt   l  
となるので、両者の連立方程式の形で表すことが出来る。
従って、単振り子を数値計算によってシミュレーションする場合、
上述の式を4次のルンゲ・クッタ法(RK4)で解けばよい。




スクリプトプログラムの使い方

各テキストフィールドに数値を入力後、「計算する」ボタンを押すと、
新しいページにデータが生成されるので、そこで、「全て選択」して、
テキストエディタを起動して、ファイルにコピーペーストし、名前を付けて保存する。
ここでは、Cドライブ以下のtempフォルダに「simplependulum.dat」という名前で保存する。
これをグラフ作成ツール「gnuplot」を用いて、グラフ化するには、
gnuplotを起動して、以下のコマンドを入力する。

角度の時間変化(1:2がRK4による数値、1:4が近似後の数値)
set terminal png
plot 'c:\temp\simplependulum.dat' using 1:2 with lines
set output 'c:\temp\theta.png'
replot 'c:\temp\simplependulum.dat' using 1:4 with lines
exit

角速度の時間変化(1:3がRK4による数値、1:5が近似後の数値)
set terminal png
plot 'c:\temp\simplependulum.dat' using 1:3 with lines
set output 'c:\temp\omega.png'
replot 'c:\temp\simplependulum.dat' using 1:5 with lines
exit

位相空間での軌道(2:3がRK4による数値、4:5が近似後の数値)
set terminal png
plot 'c:\temp\simplependulum.dat' using 2:3 with lines
set output 'c:\temp\phasespace.png'
replot 'c:\temp\simplependulum.dat' using 4:5 with lines
exit




具体例

ここでは、一律にn=10000, Δt=0.001, l=9.8 に固定し、
角度の初期条件を5[deg], 30[deg], 45[deg], 60[deg], 90[deg]と変化させて、
横軸が時間t、縦軸が角度θで表された「角度の時間変化」、
横軸が時間t、縦軸が角速度ωで表された「角速度の時間変化」、
横軸が角度θ 、縦軸が角速度ωで表された、「位相空間での軌道」
をグラフ作成ツール「gnuplot」を用いて、グラフ化する。




角度の時間変化

θ0=5[deg]
θ0=30[deg]
θ0=45[deg]
θ0=60[deg]
θ0=90[deg]




角速度の時間変化

θ0=5[deg]
θ0=30[deg]
θ0=45[deg]
θ0=60[deg]
θ0=90[deg]




位相空間での軌道

θ0=5[deg]
θ0=30[deg]
θ0=45[deg]
θ0=60[deg]
θ0=90[deg]




考察

「角度の時間変化」、「角速度の時間変化」、「位相空間での軌道」
のいずれの結果を見ても、角度の初期条件が大きくなるに伴い、
4次のルンゲ・クッタ法(RK4)による数値と、
「sinθθ」と近似した場合の数値は、
次第にズレてくる。従って、角度の初期条件が大きい場合は、
「sinθθ」という近似は、不適当であることが示唆される。


Shadow Academy トップへ戻る

inserted by FC2 system