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

繰り返しの回数：n＝ (n×Δt 秒後に終了)
時間の刻み幅：Δt＝ [s] (n×Δt 秒後に終了)
振り子のひもの長さ：l＝ [m] 参考). 重力加速度g＝9.8[m/s]
角度の初期条件：θ0＝ [deg] 孤度法([rad / ラジアン])ではなく、
度数法([deg / °])で入力すること。
参考). 角速度の初期条件：ω0＝0[rad/s]
データ出力形式： CSV(半角コンマ区切り)
SSV(半角スペース区切り)
データ出力欄：


解説

単振り子の方程式は、二階微分方程式：

　 d2θ 　
l
＋gsinθ＝0
　 dt2 　
で表されるが、
dθ 　
＝ω
dt 　
と置くと、
dω 　 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θ ≈ θ」という近似は、不適当であることが示唆される。



ソースコード（新版）

simplependulum.html

<form><!-- getElementByIdを使用する為、name属性を指定しない -->
<table border="5" style="margin-left:auto;margin-right:auto;text-align:center;">
  <tr>
    <td style="text-align:right;">繰り返しの回数：<i>n</i>＝</td>
    <td style="text-align:left;"><input type="text" id="n" value="10000"></td>
    <td style="text-align:left;">(<i>n</i>&times;<i>&Delta;t</i> 秒後に終了)</td>
  </tr>
  <tr>
    <td style="text-align:right;">時間の刻み幅：<i>&Delta;t</i>＝</td>
    <td style="text-align:left;"><input type="text" id="dt" value="0.001">[s]</td>
    <td style="text-align:left;">(<i>n</i>&times;<i>&Delta;t</i> 秒後に終了)</td>
  </tr>
  <tr>
    <td style="text-align:right;">振り子のひもの長さ：<i>l</i>＝</td>
    <td style="text-align:left;"><input type="text" id="l" value="9.8">[m]</td>
    <td style="text-align:left;">参考). 重力加速度<i>g</i>＝9.8[m/s]</td>
  </tr>
  <tr>
    <td style="text-align:right;">角度の初期条件：<i>&theta;</i><sub>0</sub>＝</td>
    <td style="text-align:left;"><input type="text" id="theta0" value="45">[deg]</td>
    <td style="text-align:left;">孤度法([rad / ラジアン])ではなく、<br>
度数法([deg / °])で入力すること。<br>
参考). 角速度の初期条件：<i>&omega;</i><sub>0</sub>＝0[rad/s]</td>
  </tr>
  <tr>
    <td style="text-align:right;">データ出力形式：</td>
    <td style="text-align:left;" colspan="2">
      <input type="radio" name="separator" id="csv">CSV(半角コンマ区切り)<br>
      <input type="radio" name="separator" id="ssv" checked>SSV(半角スペース区切り)<br>
      <input type="button" value="計算する" onClick="spcalc()" style="width:200px;height:40px;"></td>
  </tr>
  <tr>
    <td style="text-align:right;">データ出力欄：</td>
    <td style="text-align:left;" colspan="2"><textarea id="data" cols="50" rows="20"></textarea></td>
  </tr>
</table>
</form>

simplependulum.js

var spcalc = function(){
  var n = parseInt(document.getElementById("n").value); // 繰り返しの回数(n*dt秒後に終了)
  var dt = parseFloat(document.getElementById("dt").value); // 時間の刻み幅(n*dt秒後に終了)
  var g = 9.8; // 重力加速度
  var l = parseFloat(document.getElementById("l").value); // 振り子のひもの長さ
  var theta0 = parseFloat(document.getElementById("theta0").value); // 角度の初期条件

  var t,theta,omega,i; // 時刻、角度、角速度、ループ変数を宣言
  var j1,j2,j3,j4,k1,k2,k3,k4; // 4次のルンゲクッタ法で使用する変数を宣言

  var data = ""; // データを表す変数を空の初期値と共に宣言
  var separator = ""; // データの区切り文字

  // データの区切り文字の種類を選択して取得
  if(document.getElementById("csv").checked == true){ separator=","; } // 半角コンマ区切り
  if(document.getElementById("ssv").checked == true){ separator=" "; } // 半角スペース区切り

  theta0 = theta0*Math.PI/180; // 度数法から孤度法への変換
  t = 0; // 時刻を初期化
  theta = theta0; // 角度にその初期条件を代入
  omega = 0; // 角速度の初期条件は0とする
  thetaapprox = theta0*Math.cos(Math.sqrt(g/l)*t); // 近似後の角度
  omegaapprox = -theta0*Math.sqrt(g/l)*Math.sin(Math.sqrt(g/l)*t); // 近似後の角速度

  data = t + separator + Math.round(1000*(theta*180/Math.PI))/1000
           + separator + Math.round(1000*(omega*180/Math.PI))/1000
           + separator + Math.round(1000*(thetaapprox*180/Math.PI))/1000
           + separator + Math.round(1000*(omegaapprox*180/Math.PI))/1000;
  // 時刻t=0での、角度、角速度、近似後の角度、近似後の角速度
  // (孤度法から度数法へ変換して、小数点以下第3位未満を四捨五入)
  // をデータを表す変数dataに代入して、データを上書きする。

  for(i=0;i<n;++i){
    j1=dt*omega;
    k1=dt*(-g/l)*Math.sin(theta);
    j2=dt*(omega+k1/2);
    k2=dt*(-g/l)*Math.sin(theta+j1/2);
    j3=dt*(omega+k2/2);
    k3=dt*(-g/l)*Math.sin(theta+j2/2);
    j4=dt*(omega+k3);
    k4=dt*(-g/l)*Math.sin(theta+j3);

    t=Math.round(1000*(t+dt))/1000;
    theta=theta+(j1+2*j2+2*j3+j4)/6;
    omega=omega+(k1+2*k2+2*k3+k4)/6;
    thetaapprox=theta0*Math.cos(Math.sqrt(g/l)*t); // 近似後の角度
    omegaapprox=-theta0*Math.sqrt(g/l)*Math.sin(Math.sqrt(g/l)*t); // 近似後の角速度

    data += "\n" + t
           + separator + Math.round(1000*(theta*180/Math.PI))/1000
           + separator + Math.round(1000*(omega*180/Math.PI))/1000
           + separator + Math.round(1000*(thetaapprox*180/Math.PI))/1000
           + separator + Math.round(1000*(omegaapprox*180/Math.PI))/1000;
    // 時刻t=0での、角度、角速度、近似後の角度、近似後の角速度
    // (孤度法から度数法へ変換して、小数点以下第3位未満を四捨五入)
    // をデータを表す変数dataに代入して、データを上書きする。
  }

  document.getElementById("data").value = data; // テキストエリアへ出力する
};


ソースコード（旧版）

simplependulum.html

<form name="spform">
<table border="5" style="margin-left:auto;margin-right:auto;text-align:center;">
<tr>
<td style="text-align:right;">繰り返しの回数：<i>n</i>＝</td>
<td style="text-align:left;"><input type="text" name="n"></td>
<td style="text-align:left;">(<i>n</i>&times;<i>&Delta;t</i> 秒後に終了)</td>
</tr>
<tr>
<td style="text-align:right;">時間の刻み幅：<i>&Delta;t</i>＝</td>
<td style="text-align:left;"><input type="text" name="dt">[s]</td>
<td style="text-align:left;">(<i>n</i>&times;<i>&Delta;t</i> 秒後に終了)</td>
</tr>
<tr>
<td style="text-align:right;">振り子のひもの長さ：<i>l</i>＝</td>
<td style="text-align:left;"><input type="text" name="l">[m]</td>
<td style="text-align:left;">参考). 重力加速度<i>g</i>＝9.8[m/s]</td>
</tr>
<tr>
<td style="text-align:right;">角度の初期条件：<i>&theta;</i><sub>0</sub>＝</td>
<td style="text-align:left;"><input type="text" name="theta0">[deg]</td>
<td style="text-align:left;">孤度法([rad / ラジアン])ではなく、<br>
度数法([deg / °])で入力すること。<br>
参考). 角速度の初期条件：<i>&omega;</i><sub>0</sub>＝0[rad/s]</td>
</tr>
</table>
<input type="button" value="計算する" onClick="spcalc()">
</form>

simplependulum.js

function spcalc(){
  var n=parseInt(document.spform.n.value);
  /* 繰り返しの回数(n*dt秒後に終了) */
  var dt=parseFloat(document.spform.dt.value);
  /* 時間の刻み幅(n*dt秒後に終了) */
  var g=9.8; /* 重力加速度 */
  var l=parseFloat(document.spform.l.value);
  /* 振り子のひもの長さ */
  var theta0=parseFloat(document.spform.theta0.value);
  /* 角度の初期条件 */

  var t,theta,omega,i;
  /* 時刻、角度、角速度、ループ変数を宣言 */
  var j1,j2,j3,j4,k1,k2,k3,k4;
  /* 4次のルンゲクッタ法で使用する変数を宣言 */

  theta0=theta0*Math.PI/180;
  /* 度数法から孤度法への変換 */

  t=0; /* 時刻を初期化 */
  theta=theta0; /* 角度にその初期条件を代入 */
  omega=0; /* 角速度の初期条件は0とする */
  thetaapprox=theta0*Math.cos(Math.sqrt(g/l)*t);
  /* 近似後の角度 */
  omegaapprox=-theta0*Math.sqrt(g/l)*Math.sin(Math.sqrt(g/l)*t);
  /* 近似後の角速度 */

  document.writeln("<html>");
  document.writeln("<head>");
  document.writeln("<title>単振り子</title>");
  document.writeln("</head>");
  document.writeln("<body bgcolor='#ffd700' text='#000000'>");

  document.writeln("<pre>");
//  document.writeln(t," ",theta," ",omega);
  document.writeln(t," ",Math.round(1000*(theta*180/Math.PI))/1000,
    " ",Math.round(1000*(omega*180/Math.PI))/1000," ",
    Math.round(1000*(thetaapprox*180/Math.PI))/1000,
    " ",Math.round(1000*(omegaapprox*180/Math.PI))/1000);
  /* 時刻t=0での、角度と角速度(孤度法から度数法へ変換して出力)
     さらに、小数点以下第3位未満を四捨五入 */
  for(i=0;i<n;++i){
    j1=dt*omega;
    k1=dt*(-g/l)*Math.sin(theta);
    j2=dt*(omega+k1/2);
    k2=dt*(-g/l)*Math.sin(theta+j1/2);
    j3=dt*(omega+k2/2);
    k3=dt*(-g/l)*Math.sin(theta+j2/2);
    j4=dt*(omega+k3);
    k4=dt*(-g/l)*Math.sin(theta+j3);

    t=Math.round(1000*(t+dt))/1000;
    theta=theta+(j1+2*j2+2*j3+j4)/6;
    omega=omega+(k1+2*k2+2*k3+k4)/6;
    thetaapprox=theta0*Math.cos(Math.sqrt(g/l)*t);
    /* 近似後の角度 */
    omegaapprox=-theta0*Math.sqrt(g/l)*Math.sin(Math.sqrt(g/l)*t);
    /* 近似後の角速度 */

//    document.writeln(t," ",theta," ",omega);
  document.writeln(t," ",Math.round(1000*(theta*180/Math.PI))/1000,
    " ",Math.round(1000*(omega*180/Math.PI))/1000," ",
    Math.round(1000*(thetaapprox*180/Math.PI))/1000,
    " ",Math.round(1000*(omegaapprox*180/Math.PI))/1000);
  /* 時刻tでの、角度と角速度(孤度法から度数法へ変換して出力)
     さらに、小数点以下第3位未満を四捨五入 */
  }
  document.writeln("</pre>");

  document.writeln("</body>");
  document.write("</html>");
}


参考文献

  1. 「理系のためのJava入門」（工学社、2006年）


Shadow Academy トップへ戻る

inserted by FC2 system