HTML5・二重振り子のシミュレーション

初期条件:
時間の刻み幅:Δt [s]  
質点1の質量:m1 [kg] 質点1の色:
質点2の質量:m2 [kg] 質点2の色:
振り子の棒1の長さ:l1 [m]  
振り子の棒2の長さ:l2 [m]  
角度1の初期条件:θ10 [deg] 孤度法([rad / ラジアン])ではなく、
度数法([deg / °])で入力すること。
参考). 角速度1の初期条件:
ω10=0[rad/s]
角度2の初期条件:θ20 [deg] 孤度法([rad / ラジアン])ではなく、
度数法([deg / °])で入力すること。
参考). 角速度2の初期条件:
ω20=0[rad/s]
お持ち帰りコーナー: データ生成:
「PNG画像を保存する」ボタンを押すと、以下に画像が生成される。
繰り返しの回数:n [s] (n×Δt 秒後に終了)
入力データ形式:



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

各テキストフィールドに数値を入力後、「x座標・y座標を計算する」ボタンを押すと、テキストエリアに
(時刻,質点1のx座標,質点1のy座標,質点2のx座標,質点2のy座標)の形式で、データが生成されるので、
そこで、「全て選択」して、テキストエディタを起動して、ファイルにコピーペーストし、名前を付けて保存する。
ここでは、Cドライブ以下のtempフォルダに「doublependulumxy.dat」という名前で保存する。

「角度・角速度を計算する」ボタンを押すと、テキストエリアに
(時刻,質点1の角度,質点1の角速度,質点2の角度,質点2の角速度)の形式で、データが生成されるので、
そこで、「全て選択」して、テキストエディタを起動して、ファイルにコピーペーストし、名前を付けて保存する。
ここでは、Cドライブ以下のtempフォルダに「doublependulumphase.dat」という名前で保存する。

「エネルギーを計算する」ボタンを押すと、テキストエリアに
(時刻,質点1の運動エネルギーK1,質点1のポテンシャルエネルギーU1,
質点2の運動エネルギーK2,質点2のポテンシャルエネルギーU2)の形式で、データが生成されるので、
そこで、「全て選択」して、テキストエディタを起動して、ファイルにコピーペーストし、名前を付けて保存する。
ここでは、Cドライブ以下のtempフォルダに「doublependulumenergy.dat」という名前で保存する。

これをグラフ作成ツール「wgnuplot」を用いて、グラフ化するには、
wgnuplotを起動して、例えば、以下のコマンドを入力する。

質点のx座標及びy座標の軌跡(2:3が質点1の軌跡、4:5が質点2の軌跡)

set terminal png
plot 'c:\temp\doublependulumxy.dat' using 2:3 with lines
set output 'c:\temp\doublependulumxy.png'
replot 'c:\temp\doublependulumxy.dat' using 4:5 with lines
exit

質点の角度の時間変化(1:2が質点1の角度の時間変化、1:4が質点2の角度の時間変化)

set terminal png
plot 'c:\temp\doublependulumphase.dat' using 1:2 with lines
set output 'c:\temp\doublependulumphase.png'
replot 'c:\temp\doublependulumphase.dat' using 1:4 with lines
exit

各エネルギーの時間変化(1:2がK1の時間変化、1:3がU1の時間変化、1:4がK2の時間変化、1:5がU2の時間変化)

set terminal png
plot 'c:\temp\doublependulumenergy.dat' using 1:2 with lines
replot 'c:\temp\doublependulumenergy.dat' using 1:3 with lines
replot 'c:\temp\doublependulumenergy.dat' using 1:4 with lines
set output 'c:\temp\doublependulumenergy.png'
replot 'c:\temp\doublependulumenergy.dat' using 1:5 with lines
exit

また、「運動の様子を描画する」ボタンを押すと、
振り子の棒1、棒2及び、質点1、質点2の運動の軌跡が、キャンバスに描画される。
ブラウザ、或いは、そのバージョンがHTML5に対応している必要があることに注意。
質点1及び、質点2は、描画する色を選択することができる。
また、その質量に比例して、球の半径が大きくなるように描いた。
動画として保存したい場合には、Bandicam等のソフトによる、動画キャプチャをオススメする。




二重振り子のシミュレーションの作り方

二重振り子のシミュレーションを作る際に必要なのは、以下の三点である。

解析力学
ラグランジアンから、ラグランジュの運動方程式を導出する際に必要。
詳細は、二重振り子の運動方程式を参照。
数値計算
4次のルンゲ・クッタ法を用いて、連立微分方程式を解く際に必要。
GUIプログラミング
描画する際に必要。以前は、Javaアプレットが主流だったが、
今後は、HTML5+JavaScriptが主流になっていくと思われる。




ソースコード

dpsimulation.html

<form><!-- getElementByIdを使用する為、name属性を指定しない -->
<table border="1" style="margin-left:auto;margin-right:auto;text-align:center;">
  <tr>
    <td style="text-align:center;" rowspan="8">
      <canvas id="doublependulum" width="500" height="500"></canvas></td>
    <th style="text-align:center;" colspan="3">初期条件:</th>
  </tr>
  <tr>
    <td style="text-align:right;">時間の刻み幅:<i>&Delta;t</i>=</td>
    <td style="text-align:left;">
      <input type="text" id="dt" value="0.01" style="text-align:right;">[s]</td>
    <td style="text-align:left;"> </td>
  </tr>
  <tr>
    <td style="text-align:right;">質点1の質量:<i>m</i><sub>1</sub>=</td>
    <td style="text-align:left;">
      <input type="text" id="m1" value="20" style="text-align:right;">[kg]</td>
    <td style="text-align:left;">質点1の色:
      <select id="m1color">
        <option value="#ff0000" selected>赤</option>
        <option value="#0000ff">青</option>
        <option value="#ffff00">黄</option>
        <option value="#008000">緑</option>
        <option value="#800080">紫</option>
      </select></td>
  </tr>
  <tr>
    <td style="text-align:right;">質点2の質量:<i>m</i><sub>2</sub>=</td>
    <td style="text-align:left;">
      <input type="text" id="m2" value="20" style="text-align:right;">[kg]</td>
    <td style="text-align:left;">質点2の色:
      <select id="m2color">
        <option value="#ff0000">赤</option>
        <option value="#0000ff" selected>青</option>
        <option value="#ffff00">黄</option>
        <option value="#008000">緑</option>
        <option value="#800080">紫</option>
      </select></td>
  </tr>
  <tr>
    <td style="text-align:right;">振り子の棒1の長さ:<i>l</i><sub>1</sub>=</td>
    <td style="text-align:left;">
      <input type="text" id="l1" value="100" style="text-align:right;">[m]</td>
    <td style="text-align:left;"> </td>
  </tr>
  <tr>
    <td style="text-align:right;">振り子の棒2の長さ:<i>l</i><sub>2</sub>=</td>
    <td style="text-align:left;">
      <input type="text" id="l2" value="100" style="text-align:right;">[m]</td>
    <td style="text-align:left;"> </td>
  </tr>
  <tr>
    <td style="text-align:right;">角度1の初期条件:<i>&theta;</i><sub>10</sub>=</td>
    <td style="text-align:left;">
      <input type="text" id="theta1" value="145" style="text-align:right;">[deg]</td>
    <td style="text-align:left;">孤度法([rad / ラジアン])ではなく、<br>
度数法([deg / °])で入力すること。<br>
参考). 角速度1の初期条件:<br><i>&omega;</i><sub>10</sub>=0[rad/s]</td></td>
  </tr>
  <tr>
    <td style="text-align:right;">角度2の初期条件:<i>&theta;</i><sub>20</sub>=</td>
    <td style="text-align:left;">
      <input type="text" id="theta2" value="145" style="text-align:right;">[deg]</td>
    <td style="text-align:left;">孤度法([rad / ラジアン])ではなく、<br>
度数法([deg / °])で入力すること。<br>
参考). 角速度2の初期条件:<br><i>&omega;</i><sub>20</sub>=0[rad/s]</td>
  </tr>
  <tr>
  <th style="text-align:center;">
    <input type="button" value="運動の様子を描画する" onClick="dpsimulation()"
      style="width:200px;height:40px;">
    <input type="button" value="PNG画像を保存する" onclick="saveasPNG()"
      style="width:200px;height:40px;"></th>
  <th style="text-align:center;" colspan="3">
    <input type="button" value="x座標・y座標を計算する" onClick="dpcalcxy()"
      style="width:200px;height:40px;">
    <input type="button" value="角度・角速度を計算する" onClick="dpcalcphase()"
      style="width:200px;height:40px;">
    <input type="button" value="エネルギーを計算する" onClick="dpcalcenergy()"
      style="width:200px;height:40px;"></th>
  </tr>
  <tr>
    <th style="text-align:center;">お持ち帰りコーナー:</th>
    <th style="text-align:center;" colspan="3">データ生成:</th>
  </tr>
  <tr>
    <td style="text-align:center;" rowspan="3">
「PNG画像を保存する」ボタンを押すと、以下に画像が生成される。<br>
      <div style="width:500px;height:500px;"><img id="newImg"></div></td>
    <td style="text-align:right;">繰り返しの回数:<i>n</i>=</td>
    <td style="text-align:left;">
      <input type="text" id="n" value="1000" style="text-align:right;">[s]</td>
    <td style="text-align:left;">(<i>n</i>&times;<i>&Delta;t</i> 秒後に終了)</td>
  </tr>
  <tr>
    <td style="text-align:right;">入力データ形式:</td>
    <td style="text-align:left;" colspan="2">
      <label for="csv"><input type="radio" name="separator" id="csv">CSV(半角コンマ区切り)</label>
      <label for="ssv"><input type="radio" name="separator" id="ssv" checked>SSV(半角スペース区切り)</label></td>
  </tr>
  <tr>
    <td style="text-align:center;" colspan="3">
      <textarea id="data" cols="70" rows="30"></textarea></td>
  </tr>
</table>
</form>

dpsimulation.js

// グローバル変数
var cvs; // canvas要素のオブジェクトを取得する為の変数
var cnt; // コンテキストを取得する為の変数

var w = 500; // canvasの横幅
var h = 500; // canvasの高さ
var hw = 250; // canvasの横幅の半値幅
var hh = 250; // canvasの高さの半値幅
var x1,x2,y1,y2; // 質点1及び質点2のx座標及びy座標を宣言
var x1draw,x2draw,y1draw,y2draw,m1draw,m2draw; // 描画用の変数を宣言

var bgcolor = "#c0c0c0"; // canvasの背景色
var m1color = "#ff0000"; // 質点1の色を表す変数を初期値と共に宣言
var m2color = "#0000ff"; // 質点2の色を表す変数を初期値と共に宣言

var g = 9.8; // 重力加速度
var t = 0, dt = 0.01; // 時間及び、時間の刻み幅を表す変数を初期値と共に宣言
var m1 = 20, m2 = 20; // 質点1及び、質点2の質量を表す変数を初期値と共に宣言
var l1 = 100, l2 = 100;
// 振り子の棒1の長さ及び、振り子の棒2の長さを表す変数を初期値と共に宣言
var theta1 = 145, theta2 = 145;
// 角度1の初期条件及び、角度2の初期条件を表す変数を初期値と共に宣言
var omega1 = 0,omega2 = 0;
// 角速度1及び、角速度2の初期条件は0とする

var mu,dth,cdth,sdth; // μ、Δθ、cosΔθ、sinΔθを宣言
var j11,j12,j13,j14,j21,j22,j23,j24;
var k11,k12,k13,k14,k21,k22,k23,k24;
// 4次のルンゲクッタ法で使用する変数を宣言

window.onload = function(){
  cvs = document.getElementById("doublependulum"); // canvas要素のオブジェクトを取得
  cnt = cvs.getContext("2d"); // コンテキストを取得

  setup();
  draw();
};

var setup = function(){
  dt = parseFloat(document.getElementById("dt").value);
  // 時間の刻み幅を表す文字列をテキストフィールドから取得して、変数に代入
  m1 = parseFloat(document.getElementById("m1").value);
  m2 = parseFloat(document.getElementById("m2").value);
  // 質点1及び質点2の質量をテキストフィールドの文字列を実数化して取得
  m1color = document.getElementById("m1color").value
  // 質点1の色を表す文字列をテキストフィールドから取得して、変数に代入
  m2color = document.getElementById("m2color").value
  // 質点2の色を表す文字列をテキストフィールドから取得して、変数に代入
  l1 = parseFloat(document.getElementById("l1").value);
  l2 = parseFloat(document.getElementById("l2").value);
  // 振り子の棒1の長さ及び振り子の棒2の長さをテキストフィールドの文字列を実数化して取得
  theta1 = parseFloat(document.getElementById("theta1").value);
  theta2 = parseFloat(document.getElementById("theta2").value);
  // 角度1の初期条件及び角度2の初期条件をテキストフィールドの文字列を実数化して取得

  n = parseInt(document.getElementById("n").value);
  // 繰り返しの回数をテキストフィールドから取得して、変数に代入

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

  t=0; // 時刻を初期化
  mu=1+m1/m2;
  theta1=theta1*Math.PI/180;
  theta2=theta2*Math.PI/180;
  // 度数法から孤度法への変換
  omega1=0; // 角速度1の初期条件は0とする
  omega2=0; // 角速度2の初期条件は0とする
};

// 共通の処理をまとめた関数
var draw = function(){
  cnt.fillStyle = bgcolor; // canvasの背景色(銀色)を指定
  cnt.fillRect(0,0,500,500); // canvas全体を背景色の矩形で塗りつぶす
  cnt.fillStyle = "#000000"; // x軸及び、y軸の描画色(黒色)を指定
  drawLine(0,hh,w,hh); // x軸を描画
  drawLine(hw,0,hw,h); // y軸を描画

  x1=parseInt(l1*Math.sin(theta1));
  x2=parseInt(x1+l2*Math.sin(theta2));
  y1=parseInt(l1*Math.cos(theta1));
  y2=parseInt(y1+l2*Math.cos(theta2));

  x1draw = x1 + hw;
  x2draw = x2 + hw;
  y1draw = y1 + hh;
  y2draw = y2 + hh;
  m1draw = Math.ceil(m1 / 2);
  m2draw = Math.ceil(m2 / 2);

  drawLine(hw,hh,x1draw,y1draw); // 棒1を描画
  drawLine(x1draw,y1draw,x2draw,y2draw); // 棒2を描画

  // 質点1を描画
  cnt.beginPath();
  cnt.fillStyle = m1color;
  cnt.arc(x1draw,y1draw,m1draw,0,Math.PI*2,true);
  cnt.fill();

  // 質点2を描画
  cnt.beginPath();
  cnt.fillStyle = m2color;
  cnt.arc(x2draw,y2draw,m2draw,0,Math.PI*2,true);
  cnt.fill();
};

// 直線を引く関数
var drawLine = function(x0, y0, x1, y1){
  cvs = document.getElementById("doublependulum"); // canvas要素のオブジェクトを取得
  cnt = cvs.getContext("2d"); // コンテキストを取得

  cnt.beginPath();
  cnt.moveTo(x0,y0);
  cnt.lineTo(x1,y1);
  cnt.stroke();
};

// 4次のルンゲ・クッタ法を用いて、角度・角速度の計算をする関数
var RungeKutta = function(){
  dth=theta1-theta2;
  cdth=Math.cos(dth);
  sdth=Math.sin(dth);
  j11=dt*omega1;
  j21=dt*omega2;
  k11=dt*(g*(Math.sin(theta2)*cdth-mu*Math.sin(theta1))
     -(l1*omega1*omega1*cdth+l2*omega2*omega2)*sdth)/(l1*(mu-cdth*cdth));
  k21=dt*(g*mu*(Math.sin(theta1)*cdth-Math.sin(theta2))
     -(mu*l1*omega1*omega1+l2*omega2*omega2*cdth)*sdth)/(l2*(mu-cdth*cdth));

  // omega1→(omega1+k11*0.5)
  // omega2→(omega2+k21*0.5)
  // theta1→(theta1+j11*0.5)
  // theta2→(theta2+j21*0.5)

  dth=(theta1+j11*0.5)-(theta2+j21*0.5);
  cdth=Math.cos(dth);
  sdth=Math.sin(dth);
  j12=dt*(omega1+k11*0.5);
  j22=dt*(omega2+k21*0.5);
  k12=dt*(g*(Math.sin(theta2+j21*0.5)*cdth-mu*Math.sin(theta1+j11*0.5))
     -(l1*(omega1+k11*0.5)*(omega1+k11*0.5)*cdth+l2*(omega2+k21*0.5)*(omega2+k21*0.5))*sdth)
     /(l1*(mu-cdth*cdth));
  k22=dt*(g*mu*(Math.sin(theta1+j11*0.5)*cdth-Math.sin(theta2+j21*0.5))
     -(mu*l1*(omega1+k11*0.5)*(omega1+k11*0.5)+l2*(omega2+k21*0.5)*(omega2+k21*0.5)*cdth)*sdth)
     /(l2*(mu-cdth*cdth));

  // (omega1+k11*0.5)→(omega1+k12*0.5)
  // (omega2+k21*0.5)→(omega2+k22*0.5)
  // (theta1+j11*0.5)→(theta1+j12*0.5)
  // (theta2+j21*0.5)→(theta2+j22*0.5)

  dth=(theta1+j12*0.5)-(theta2+j22*0.5);
  cdth=Math.cos(dth);
  sdth=Math.sin(dth);
  j13=dt*(omega1+k12*0.5);
  j23=dt*(omega2+k22*0.5);
  k13=dt*(g*(Math.sin(theta2+j22*0.5)*cdth-mu*Math.sin(theta1+j12*0.5))
     -(l1*(omega1+k12*0.5)*(omega1+k12*0.5)*cdth+l2*(omega2+k22*0.5)*(omega2+k22*0.5))*sdth)
     /(l1*(mu-cdth*cdth));
  k23=dt*(g*mu*(Math.sin(theta1+j12*0.5)*cdth-Math.sin(theta2+j22*0.5))
     -(mu*l1*(omega1+k12*0.5)*(omega1+k12*0.5)+l2*(omega2+k22*0.5)*(omega2+k22*0.5)*cdth)*sdth)
     /(l2*(mu-cdth*cdth));

  // (omega1+k12*0.5)→(omega1+k13)
  // (omega2+k22*0.5)→(omega2+k23)
  // (theta1+j12*0.5)→(theta1+j13)
  // (theta2+j22*0.5)→(theta2+j23)

  dth=(theta1+j13)-(theta2+j23);
  cdth=Math.cos(dth);
  sdth=Math.sin(dth);
  j14=dt*(omega1+k13);
  j24=dt*(omega2+k23);
  k14=dt*(g*(Math.sin(theta2+j23)*cdth-mu*Math.sin(theta1+j13))
     -(l1*(omega1+k13)*(omega1+k13)*cdth+l2*(omega2+k23)*(omega2+k23))*sdth)
     /(l1*(mu-cdth*cdth));
  k24=dt*(g*mu*(Math.sin(theta1+j13)*cdth-Math.sin(theta2+j23))
     -(mu*l1*(omega1+k13)*(omega1+k13)+l2*(omega2+k23)*(omega2+k23)*cdth)*sdth)
     /(l2*(mu-cdth*cdth));

  t=Math.round(1000*(t+dt))/1000;
  theta1=theta1+(j11+2*j12+2*j13+j14)/6;
  theta2=theta2+(j21+2*j22+2*j23+j24)/6;
  omega1=omega1+(k11+2*k12+2*k13+k14)/6;
  omega2=omega2+(k21+2*k22+2*k23+k24)/6;
};

var update = function(){
  RungeKutta();
  draw();
};

var dpsimulation = function(){
  setup();
  draw();
  setInterval("update()",dt*1000); // dt秒毎に実行
};

// PNG画像として保存
var saveasPNG = function(){
  var png = cvs.toDataURL("image/png");
  document.getElementById("newImg").src = png;
};

var data = ""; // データを表す変数を空の初期値と共に宣言
var separator = ""; // データの区切り文字
var i, n = 1000; // ループ変数と繰り返しの回数を初期値と共に宣言

// x座標・y座標を計算する関数
var dpcalcxy = function(){
  setup();
  x1=parseInt(l1*Math.sin(theta1));
  x2=parseInt(x1+l2*Math.sin(theta2));
  y1=parseInt(l1*Math.cos(theta1));
  y2=parseInt(y1+l2*Math.cos(theta2));

  data = t + separator + Math.round(1000*(x1))/1000
           + separator + Math.round(1000*(y1))/1000
           + separator + Math.round(1000*(x2))/1000
           + separator + Math.round(1000*(y2))/1000;
  // 時刻t=0での、質点1及び質点2のx座標及びy座標(小数点以下第3位未満を四捨五入)
  // をデータを表す変数dataに代入して、データを上書きする。

  for(i=0;i<n;++i){
    RungeKutta();
    x1=l1*Math.sin(theta1);
    x2=x1+l2*Math.sin(theta2);
    y1=l1*Math.cos(theta1);
    y2=y1+l2*Math.cos(theta2);

    data += "\n" + t + separator + Math.round(1000*(x1))/1000
                     + separator + Math.round(1000*(y1))/1000
                     + separator + Math.round(1000*(x2))/1000
                     + separator + Math.round(1000*(y2))/1000;
    // 時刻tでの、質点1及び質点2のx座標及びy座標(小数点以下第3位未満を四捨五入)
    // をデータを表す変数dataに代入して、データを上書きする。
  }

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

// 角度・角速度を計算する関数
var dpcalcphase = function(){
  setup();

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

  for(i=0;i<n;++i){
    RungeKutta();

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

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

// エネルギーを計算する関数
var dpcalcenergy = function(){
  setup();

  K1=m1*l1*l1*omega1*omega1/2;
  U1=-m1*g*l1*Math.cos(theta1);
  K2=m2*(l1*l1*omega1*omega1+l2*l2*omega2*omega2)/2
    +m2*l1*l2*omega1*omega2*Math.cos(theta1-theta2);
  U2=-m2*g*(l1*Math.cos(theta1)+l2*Math.cos(theta2));

  data = t + separator + Math.round(1000*(K1))/1000
           + separator + Math.round(1000*(U1))/1000
           + separator + Math.round(1000*(K2))/1000
           + separator + Math.round(1000*(U2))/1000;
  // 時刻t=0での、質点1及び質点2の運動エネルギーK及びポテンシャルエネルギーU
  // (小数点以下第3位未満を四捨五入)をデータを表す変数dataに代入して、データを上書きする。

  for(i=0;i<n;++i){
    RungeKutta();
    K1=m1*l1*l1*omega1*omega1/2;
    U1=-m1*g*l1*Math.cos(theta1);
    K2=m2*(l1*l1*omega1*omega1+l2*l2*omega2*omega2)/2
      +m2*l1*l2*omega1*omega2*Math.cos(theta1-theta2);
    U2=-m2*g*(l1*Math.cos(theta1)+l2*Math.cos(theta2));

    data += "\n" + t + separator + Math.round(1000*(K1))/1000
                     + separator + Math.round(1000*(U1))/1000
                     + separator + Math.round(1000*(K2))/1000
                     + separator + Math.round(1000*(U2))/1000;
    // 時刻tでの、質点1及び質点2の運動エネルギーK及びポテンシャルエネルギーU
    // (小数点以下第3位未満を四捨五入)をデータを表す変数dataに代入して、データを上書きする。
  }

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



Shadow Academy トップへ戻る

inserted by FC2 system