Javascriptで数値計算その3（Symplectic 8次で三体問題、8の字解）

三体問題の解で、8の字になるものがある。

http://www.scholarpedia.org/article/Three_body_problem

http://www.clas.kitasato-u.ac.jp/~fujiwara/nBody/rims/RIMS2003Abstract.pdf

以前はExcel VBAでシンプレクティック 8次で計算したけれど、今度はJavascript版。

nmax:

※Canvasが現れるまでOK押さないで。

プログラムソースはこちら。

function func100421() {
    var ctx = document.getElementById('c100421').getContext('2d');
   
    w = new Array(8);
    c = new Array(17);
    d = new Array(17);
    px = new Array(4);
    py = new Array(4);
    qx = new Array(4);
    qy = new Array(4);
    m = new Array(4);
   
    w[1] = 0.311790812418427;
    w[2] = -1.55946803821447;
    w[3] = -1.6789692825964;
    w[4] = 1.66335809963315;
    w[5] = -1.06458714789183;
    w[6] = 1.36934946416871;
    w[7] = 0.629030650210433;
    w[0] = 1. - 2. * (w[1] + w[2] + w[3] + w[4] + w[5] + w[6] + w[7]);
   
    d[1] = w[7];
    d[15] = w[7];
    d[2] = w[6];
    d[14] = w[6];
    d[3] = w[5];
    d[13] = w[5];
    d[4] = w[4];
    d[12] = w[4];
    d[5] = w[3];
    d[11] = w[3];
    d[6] = w[2];
    d[10] = w[2];
    d[7] = w[1];
    d[9] = w[1];
    d[8] = w[0];
    d[16] = 0.;

    c[1] = 0.5 * w[7];
    c[16] = 0.5 * w[7];
    c[2] = 0.5 * (w[7] + w[6]);
    c[15] = 0.5 * (w[7] + w[6]);
    c[3] = 0.5 * (w[6] + w[5]);
    c[14] = 0.5 * (w[6] + w[5]);
    c[4] = 0.5 * (w[5] + w[4]);
    c[13] = 0.5 * (w[5] + w[4]);
    c[5] = 0.5 * (w[4] + w[3]);
    c[12] = 0.5 * (w[4] + w[3]);
    c[6] = 0.5 * (w[3] + w[2]);
    c[11] = 0.5 * (w[3] + w[2]);
    c[7] = 0.5 * (w[2] + w[1]);
    c[10] = 0.5 * (w[2] + w[1]);
    c[8] = 0.5 * (w[1] + w[0]);
    c[9] = 0.5 * (w[1] + w[0]);
   
    n = 3;
   
    m[1] = 1.;
    m[2] = 1.;
    m[3] = 1.;
   
    qx[1] = 0.97000436;
    qy[1] = -0.24308753;
   
    px[1] = -0.5 * 0.93240737;
    py[1] = -0.5 * 0.86473146;
   
    qx[2] = -0.97000436;
    qy[2] = 0.24308753;
   
    px[2] = -0.5 * 0.93240737;
    py[2] = -0.5 * 0.86473146;
   
    qx[3] = 0.;
    qy[3] = 0.;
   
    px[3] = 0.93240737;      
    py[3] = 0.86473146;
    G = 1.;
    ctx.clearRect(0,0,300,300);
    ctx.strokeStyle ="rgb(0,0,0)";
    ctx.strokeRect(0,0,300,300);

    nmax = parseInt(document.Form100421.nvalue.value);

    G=1.0;

    t=0.;
    dt=0.01;

    ctx.strokeStyle ="rgb(255,0,0)";
    ctx.beginPath();   
    ctx.moveTo(150.+150.*qx[1]/2.,150.-150.*qy[1]/2.); 
   
    for (i=1;i<=nmax;i++) {

   
     for (k=1;k<=16;k++) {
         for (j=1;j<=3;j++) {
                qx[j] = qx[j] + dt * c[k] * px[j] / m[j];
                qy[j] = qy[j] + dt * c[k] * py[j] / m[j];
            }

         for (j=1;j<=3;j++) {
                fx = 0.;
                fy = 0.;

                for (ii=1;ii<=3;ii++) {
                   if (ii!=j) {

                     r = Math.sqrt((qx[ii] - qx[j]) * (qx[ii] - qx[j]) + (qy[ii] - qy[j]) * (qy[ii] - qy[j]));
                     fx = fx + G * m[ii] * m[j] * (qx[ii] - qx[j]) / (r*r*r);
                     fy = fy + G * m[ii] * m[j] * (qy[ii] - qy[j]) / (r*r*r);
                   }
                }
   
                px[j] = px[j] + dt * d[k] * fx;
                py[j] = py[j] + dt * d[k] * fy;
            }
        }

        t = t + dt;
        ctx.lineTo(150.+150.*qx[1]/2.,150.-150.*qy[1]/2.);   
  }
  ctx.stroke();

}