メモ: わからないながらにPlaskで遊んでみた2

結果:

Plask:

http://www.plask.org/

参考:

web雑記帳: javascriptで4次Runge-Kutta法

重振り子

コード:

var plask = require('plask');

//
// refer to:
// http://yoshifuji.blogspot.jp/2011/07/javascript4runge-kutta.html
// http://www.ne.jp/asahi/tokyo/nkgw/gakusyu/rikigaku/2zyufuriko/zyufuriko_kaisetu/zyufuriko_kaisetu.html
function solver(x1,x2,x3,x4,t,dt){  
  var k11,k12,k13,k14,k21,k22,k23,k24,k31,k32,k33,k34,k41,k42,k43,k44;  
  
  k11 = dt * func1(x1,x2,x3,x4,t);  
  k12 = dt * func2(x1,x2,x3,x4,t);  
  k13 = dt * func3(x1,x2,x3,x4,t);  
  k14 = dt * func4(x1,x2,x3,x4,t);  
   
  k21 = dt * func1(x1 + k11/2, x2 + k12/2, x3 + k13/2, x4 + k14/2, t + dt/2);  
  k22 = dt * func2(x1 + k11/2, x2 + k12/2, x3 + k13/2, x4 + k14/2, t + dt/2);  
  k23 = dt * func3(x1 + k11/2, x2 + k12/2, x3 + k13/2, x4 + k14/2, t + dt/2);  
  k24 = dt * func4(x1 + k11/2, x2 + k12/2, x3 + k13/2, x4 + k14/2, t + dt/2);  

  k31 = dt * func1(x1 + k21/2, x2 + k22/2, x3 + k23/2, x4 + k24/2, t + dt/2);  
  k32 = dt * func2(x1 + k21/2, x2 + k22/2, x3 + k23/2, x4 + k24/2, t + dt/2);  
  k33 = dt * func3(x1 + k21/2, x2 + k22/2, x3 + k23/2, x4 + k24/2, t + dt/2);  
  k34 = dt * func4(x1 + k21/2, x2 + k22/2, x3 + k23/2, x4 + k24/2, t + dt/2);  

  k41 = dt * func1(x1 + k31, x2 + k32, x3 + k33/2, x4 + k34/2, t + dt);  
  k42 = dt * func2(x1 + k31, x2 + k32, x3 + k33/2, x4 + k34/2, t + dt);  
  k43 = dt * func3(x1 + k31, x2 + k32, x3 + k33/2, x4 + k34/2, t + dt);  
  k44 = dt * func4(x1 + k31, x2 + k32, x3 + k33/2, x4 + k34/2, t + dt);  

  var tmp= [];
  tmp[0] = x1 + (k11/6 + k21/3 + k31/3 + k41/6);  
  tmp[1] = x2 + (k12/6 + k22/3 + k32/3 + k42/6);  
  tmp[2] = x3 + (k13/6 + k23/3 + k33/3 + k43/6);  
  tmp[3] = x4 + (k14/6 + k24/3 + k34/3 + k44/6);

  return tmp;  
}
var m1= 1.;
var m2= 2.;
var l1= 9.8;
var l2= 9.8;
var  g= 9.8; //...
function func1(x1,x2,x3,x4,t){  
  return x2*0.99999;
}  
function func2(x1,x2,x3,x4,t){
  var c0= x1-x3;
  var cc= Math.cos(c0);
  var cs= Math.sin(c0);

  var c1= 1.0/(m1+m2*cs*cs);
  var c2= -(m1+m2)*Math.sin(x1)/l1*g;
  var c3=  m2*cc*Math.sin(x2)/l1*g;
  var c4= -m2*cc*cs*x2*x2;
  var c5= -m2*cs*(l2/l1)*x4*x4;
  return c1*(c2+c3+c4+c5); 
}  
function func3(x1,x2,x3,x4,t){  
  return x4*0.99999;
}
function func4(x1,x2,x3,x4,t){  
  var c0= x1-x3;
  var cc= Math.cos(c0);
  var cs= Math.sin(c0);
  
  var c1= 1.0/(m1+m2*cs*cs);
  var c2= (m1+m2)*cs*Math.cos(x1)/l2*g;
  var c3= (m1+m2)*cs*Math.sin(x2)*(l1/l2)*x2*x2;
  var c4= m2*cc*cs*x4*x4;
  return c1*(c2+c3+c4);
}  
//

var path= new PlaskRawMac.SkPath;

var dots1= [];
var dots2= [];
var dotsN= 90.;

var xs = [];
var t;
var dt;

var flag;
plask.simpleWindow({
  settings: {
    width: 400,
    height: 400
  },

  init: function() {
    var canvas = this.canvas, paint = this.paint;
    this.framerate(60);

    flag= false;
    this.on('keyDown', function(e) {
      switch (e.str) {
        case ' ': flag= true; break;
      }
    });
    this.on('keyUp', function(e) {
      switch (e.str) {
        case ' ': flag= false; break;
      }
    });

    for(var i=0; i<dotsN; i++) {
      dots1.push({x:0, y:0});
      dots2.push({x:0, y:0});
    }
    xs[0]= Math.PI/3;
    xs[1]= 0.1;
    xs[2]= Math.PI;
    xs[3]= -0.1;
    dt= 1./720.;
    t= 0.0;
  },

  draw: function() {
    var canvas = this.canvas, paint = this.paint;
    canvas.clear(0, 0, 0, 255);    

    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;
    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;
    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;
    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;

    var y= 200;
    var x= 200;
    var th1= xs[0];
    var th2= xs[2];
    y+= 10.*l1*Math.cos(th1);
    x+= 10.*l1*Math.sin(th1);
    dots1.push({x:x, y:y});
    dots1.shift();
    y+= 10.*l2*Math.cos(th2);
    x+= 10.*l2*Math.sin(th2);
    dots2.push({x:x, y:y});
    dots2.shift();

    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;
    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;
    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;
    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;

    y= 200;
    x= 200;
    var th1= xs[0];
    var th2= xs[2];
    y+= 10.*l1*Math.cos(th1);
    x+= 10.*l1*Math.sin(th1);
    dots1.push({x:x, y:y});
    dots1.shift();
    y+= 10.*l2*Math.cos(th2);
    x+= 10.*l2*Math.sin(th2);
    dots2.push({x:x, y:y});
    dots2.shift();


    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;
    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;
    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;
    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;

    path.reset();
    y= 200;
    x= 200;
    path.moveTo(x, y);
    var th1= xs[0];
    var th2= xs[2];
    y+= 10.*l1*Math.cos(th1);
    x+= 10.*l1*Math.sin(th1);
    path.lineTo(x, y);
    dots1.push({x:x, y:y});
    dots1.shift();
    y+= 10.*l2*Math.cos(th2);
    x+= 10.*l2*Math.sin(th2);
    path.lineTo(x, y);
    dots2.push({x:x, y:y});
    dots2.shift();


    if (flag) {
      paint.setStroke();
      paint.setColor(125,125,125,200);
      canvas.drawPath(paint, path);
    }

    paint.setFill();
    for(var i= 0; i<dotsN; i++) {
      var p1= dots1[i];
      paint.setColor(255, 125, 125, 255/dotsN*i);
      canvas.drawCircle(paint, p1.x, p1.y, i/dotsN*10.);
      var p2= dots2[i];
      paint.setColor(125, 125, 355, 255/dotsN*i);
      canvas.drawCircle(paint, p2.x, p2.y, i/dotsN*10.);
    }
  }
});