プログラミング

2重(N重)振り子の数値シミュレーション – Javascriptで計算から描画まで

投稿日:

2重振り子を数値シミュレーションをJavascriptでやってみます。Javascriptでやる利点は計算後の結果表示アニメーションまで容易に行える事だと言えます。
2重振り子の解法に関する記事はWeb上に多数ありますが、本稿の特徴は

 ・2重ではなくN重で定式化して解く
 ・連立一次方程式をライブラリ等を使用せずLDL分解をして解く
 ・コードを全部公開

となっています。

N重振り子の定式化と解法

これはWeb上の情報を使わせてもらいましょう。
こちらの記事が参考になります。

https://note.com/sciencecafe_mc2/n/ne5077e5a62f1

上記の内容に1点だけ修正を加えます。
最終式導出過程で両辺をli で除算していますが、これは本稿では行わないです。理由は行列の対称性を保つためです。対称行列ならば連立一次方程式を解く際にLU分解ではなく、LDL分解を使用する事ができます。LDL分解のほうがコードが単純になると思われます。

また、ステップ幅を大きめに取りたい場合は誤差を抑えるためにホイン法やルンゲクッタ法を用いると良いと思われます。本稿では単純に1回の計算で増加量を求めています。

描画の工夫

全軌跡と現在の振り子の動きを表示する事を考えます。Canvasを2つ用意して、1つ目のcanvasに軌跡を、2つ目のcanvasに現在の振り子や状態量のテキストを表示させます。軌跡はステップごとに書き足していけば良いです。現在の振り子などは描画の前にclearRectにより前回の表示内容を全て消します。

アイデア次第でより良いアニメーションが作れると思います。自由にコードを改変してみて下さい。

ソースコード

mPendulum.html

<!DOCTYPE html>
<html lang="en">
<head>
  <meta charset="UTF-8">
  <meta name="viewport" content="width=device-width, initial-scale=1.0">
  <meta http-equiv="X-UA-Compatible" content="ie=edge">
  <title>Document</title>
  <script defer src="mPendulum.js"></script>
  <style>
    body {
      margin: 0;
      padding: 0;
      width: 100vw;
      height: 100vh;
      display: flex; 
      justify-content: center; 
      align-items: center;
    }

    canvas {
      position: absolute;
    }
  </style>
</head>

<body>
  
</body>
</html>


mPendulum.js

var mCanvasWidth = 768
var mCanvasHeight = 768
var Ox = mCanvasWidth/2
var Oy = mCanvasHeight/2

var canvas_bg = document.createElement("canvas");
canvas_bg.width = mCanvasWidth;
canvas_bg.height = mCanvasHeight;
document.body.append(canvas_bg)
var ctx_bg = canvas_bg.getContext('2d');
ctx_bg.fillStyle = "rgb(0,0,0)";
ctx_bg.fillRect(0, 0, mCanvasWidth, mCanvasHeight);

  // Grid
  ctx_bg.strokeStyle = "rgb(255,255,255)";
  ctx_bg.beginPath();
  ctx_bg.moveTo(0, mCanvasHeight/2);
  ctx_bg.lineTo(mCanvasWidth, mCanvasHeight/2);
  ctx_bg.stroke();

  ctx_bg.beginPath();
  ctx_bg.moveTo(mCanvasWidth/2, 0);
  ctx_bg.lineTo(mCanvasWidth/2, mCanvasHeight);
  ctx_bg.stroke();

var canvas_moving = document.createElement("canvas");
canvas_moving.width = mCanvasWidth;
canvas_moving.height = mCanvasHeight;
document.body.append(canvas_moving)
var ctx_m = canvas_moving.getContext('2d');

var mVideoFPS = 60;
var mFrameCount = 0

var g = 9.8
var dT = 1.0E-5;//1.0E-5;
var M = 1.0

var N_pendulum_max = 32;//Max
var N_pendulum = 2; //N_pendulum_max;//Max
var mLinkLength = 1.0;
var vec_theta = Array(N_pendulum_max).fill(Math.PI*0.95)
var vec_thetaDot = Array(N_pendulum_max).fill(0.0)
var vec_L = Array(N_pendulum_max).fill(mLinkLength)
var vec_pos = Array(N_pendulum_max*2).fill(0.0)
var vec_pre_pos = Array(N_pendulum_max*2).fill(0.0)


for(var i=0;i<N_pendulum_max;i++){
  //vec_L[i] = mLinkLength;
}//

var L_ = 0
for(var i=0;i<N_pendulum;i++){
  L_ += vec_L[i]
}//

var mScale = mCanvasWidth / (L_*2*1.05)  
console.log("mScale:"+mScale)


setInterval(async () => {

  // Calculation each angle
  // for(var i=0;i<N_pendulum;i++){
  //   vec_theta[i] += 0.01 * (i+1)
  //   //console.log("vec_theta["+i+"]:"+vec_theta[i])
  // }//

  var N = N_pendulum;

  var A = new Array(N)
  var B = new Array(N)
  var Lmat = new Array(N)
  for (var i=0;i<N;i++){
    A[i] = new Array(N).fill(0.0)
    B[i] = new Array(N).fill(0.0)
    Lmat[i] = new Array(N).fill(0.0)
  }

  var b = new Array(N).fill(0.0)
  var D = new Array(N).fill(0.0)
  var x = new Array(N).fill(0.0)
  var y = new Array(N).fill(0.0)

  for (var t=0;t<1000;t++)
  {
      for (var i=0;i<N;i++){
          b[i]=0.0;
          D[i]=0.0;
          x[i]=0.0;
          y[i]=0.0;
          for (var j=0;j<N;j++) {
              A[i][j] = 0.0;
              B[i][j] = 0.0;
              Lmat[i][j] = 0.0;
          }//
          Lmat[i][i] = 1.0;
      }//

      //Log.i("MyApp","A,B");
      for (var i=0;i<N;i++){
          for (var j=0;j<N;j++) {
              for (var k=Math.max(i,j);k<N;k++) {
                  A[i][j] += M;
              }//
              B[i][j] = A[i][j];
              //A[i][j] *= L * Math.cos( vec_theta[i] - vec_theta[j] );
              A[i][j] *= vec_L[i]*vec_L[j]*Math.cos( vec_theta[i] - vec_theta[j] );

              if (i==j){
                  B[i][j] *= g*vec_L[i]*Math.sin( vec_theta[i] );
              }//
              else {
                  //B[i][j] *= L*vec_thetaDot[j]*vec_thetaDot[j]*Math.sin( vec_theta[i] - vec_theta[j] );
                  B[i][j] *= vec_L[i]*vec_L[j]*vec_thetaDot[j]*vec_thetaDot[j]*Math.sin( vec_theta[i] - vec_theta[j] );
              }//

          }//
      }//
      //Log.i("MyApp","A,B End");

      for (var i=0;i<N;i++) {
          for (var j = 0; j < N; j++) {
              b[i]+= -B[i][j];
          }//
      }//


      D[0] = A[0][0];

      //Log.i("MyApp","D,Lmat");
      for (var k=1;k<N;k++) {
          //Log.i("MyApp","k="+k);
          for (var i = 0; i <= k-1; i++) {
              //Log.i("MyApp","i="+i);
              Lmat[k][i] = A[k][i];

              for (var j = 0; j <= i-1; j++) {
                  //Log.i("MyApp","j="+j);
                  Lmat[k][i] -= Lmat[k][j]*Lmat[i][j]*D[j];
              }//
              Lmat[k][i] /= D[i];
          }//

          D[k] = A[k][k];
          for (var i = 0; i <= k-1; i++) {
              D[k] -= Lmat[k][i]*Lmat[k][i]*D[i];
          }//

      }//
      //Log.i("MyApp","D,Lmat End");

      y[0] = b[0];
      for (var k=1;k<N;k++) {
          y[k] = b[k];
          for (var i = 0; i <= k-1; i++) {
              y[k] -= Lmat[k][i] * y[i];
          }//
      }//

      for (var k=0;k<N;k++) {
          y[k] /= D[k];
      }

      x[N-1] = y[N-1];
      for (var k=N-2;k>=0;k--) {
          x[k] = y[k];
          for (var i = k+1; i < N; i++) {
              x[k] -= Lmat[i][k]*x[i];
          }//
      }//

      for (var k=0;k<N;k++) {
          vec_thetaDot[k] += x[k]*dT;
          vec_theta[k] += vec_thetaDot[k]*dT;
      }

      //mStep += 1;

  }//for


  // Energy monitoring
  var K = new Array(N).fill(0.0)
  var U = new Array(N).fill(0.0)

  for (var i=0;i<N;i++){
      K[i] = 0.0;
      U[i] = 0.0;
      for (var j=0;j<=i;j++){
          for (var k=j;k<=i;k++){
              if (j==k){
                  K[i] += 0.5*vec_L[j]*vec_L[k]*vec_thetaDot[j]*vec_thetaDot[k]*Math.cos( vec_theta[j]-vec_theta[k]);
              }
              else{
                  K[i] += vec_L[j]*vec_L[k]*vec_thetaDot[j]*vec_thetaDot[k]*Math.cos( vec_theta[j]-vec_theta[k]);
              }
          }//k
          U[i] += vec_L[j] * Math.cos( vec_theta[j] );
      }//j
      K[i] *= M;
      U[i] *= M*g;
  }//i

  var E = 0.0;
  for (var i=0;i<N;i++){
      E += K[i] - U[i];
  }//

  

  // Update each position
  var offset_angle = - Math.PI /2 
  vec_pos[0] = vec_L[0] * Math.cos( vec_theta[0] + offset_angle )
  vec_pos[1] = vec_L[0] * Math.sin( vec_theta[0] + offset_angle)
  //console.log("vec_pos[0]:"+[vec_pos[0], vec_pos[1]])
  for(var i=1;i<N_pendulum;i++){
    vec_pos[i*2+0] = vec_pos[(i-1)*2+0] + vec_L[i] * Math.cos( vec_theta[i] + offset_angle )
    vec_pos[i*2+1] = vec_pos[(i-1)*2+1] + vec_L[i] * Math.sin( vec_theta[i] + offset_angle )
    //console.log("vec_pos["+i+"]:"+[vec_pos[i*2+0], vec_pos[i*2+1]])
  }//

  if(mFrameCount==0){
    for(var i=0;i<N_pendulum;i++){
      vec_pre_pos[i*2+0] = vec_pos[i*2+0]
      vec_pre_pos[i*2+1] = vec_pos[i*2+1]
    }
  }

  // Path
  //ctx_bg.strokeStyle = "rgb(255,0,0)";
  

  for(var i=0;i<N_pendulum;i++){
    ctx_bg.strokeStyle = "#18EBF9";
    if( vec_pre_pos[i*2+1] <= 0 ){
      ctx_bg.strokeStyle = "#fff100";
    }
  
    ctx_bg.beginPath();
    ctx_bg.moveTo(Ox+vec_pre_pos[i*2+0]*mScale, Oy-vec_pre_pos[i*2+1]*mScale);
    ctx_bg.lineTo(Ox+vec_pos[i*2+0]*mScale, Oy-vec_pos[i*2+1]*mScale);
    ctx_bg.stroke();

    vec_pre_pos[i*2+0] = vec_pos[i*2+0]
    vec_pre_pos[i*2+1] = vec_pos[i*2+1]
  }




  ctx_m.clearRect(0, 0, mCanvasWidth, mCanvasHeight)

  // Link
  //ctx_m.strokeStyle = "#18EBF9"; //"rgb(255,255,60)";
  ctx_m.strokeStyle = "rgb(50,50,255)" //"blue"
  ctx_m.beginPath();
  ctx_m.moveTo(Ox, Oy);
  ctx_m.lineTo(Ox+vec_pos[0]*mScale, Oy-vec_pos[1]*mScale);
  ctx_m.stroke();
  for(var i=1;i<N_pendulum;i++){    
    //ctx_m.strokeStyle = "#fff100" //"rgb(255,255,60)";
    ctx_m.strokeStyle = "yellow"
    ctx_m.beginPath();
    ctx_m.moveTo(Ox+vec_pos[(i-1)*2+0]*mScale, Oy-vec_pos[(i-1)*2+1]*mScale);
    ctx_m.lineTo(Ox+vec_pos[i*2+0]*mScale, Oy-vec_pos[i*2+1]*mScale);
    ctx_m.stroke();
  }
  
  // Mass
  for(var i=0;i<N_pendulum;i++){
    ctx_m.fillStyle = "rgb(50,50,255)" ;
    if(i==1){
      ctx_m.fillStyle = "yellow"
    }
    ctx_m.beginPath () ;
    ctx_m.arc( Ox+vec_pos[i*2+0]*mScale, Oy-vec_pos[i*2+1]*mScale, 10, 0, Math.PI*2, false ) ;
    ctx_m.fill() ;
  }

  // Energy
  ctx_m.font = "15px serif";
  ctx_m.fillStyle = 'rgb(255,0,0)'
  ctx_m.fillText(" E:"+ E.toExponential(3), 10, 20);  


  mFrameCount += 1

}, 1000/mVideoFPS)

-プログラミング

執筆者:

関連記事

Microsoft AzureでNode.jsゲームをデプロイ – web.configを忘れるな

Node.js+Three.jsでフォートナイトのような建築FPSゲームを自作しております。Paasサービスは手軽で太っ腹なHerokuを利用していますがサーバーがアメリカにあるため通信ラグに悩まされ …

Double (Multiple) pendulum simulation by Javascript

Source Code mPendulum.html <!DOCTYPE html> <html lang="en"> <head> <m …

Node.js+Three.js でTPSツクール

大人気ゲーム「フォートナイト」のようなTPSゲームを簡易的に作成できるサンプルプログラムを公開しました。MITライセンスとしていますので自由にコード改変して遊んでみてください。 目次 デモ動画ダウンロ …

Javascriptでまばたき検知 face-api.js landmarks

顔認識の技術の発展により、今日ではまばたきの検知まで容易に行えるようです。例えばPython、OpenCV、dlibを使用したものではhttps://www.pyimagesearch.com/201 …

Node.jsゲーム公開はHerokuが良い – 初めてのゲーム開発からデプロイまでの記録

ゲーム開発初心者が簡単なNode.jsオンラインゲームを作成して公開したい場合はHerokuで公開するのが良いです。Node.jsゲームの公開手段は様々ありますが、Herokuは無料プランで十分にテス …

スポンサーリンク