Javascriptで数値計算その2（Runge-Kutta8次のDormand-Prince法）

昨日、一応うちの環境ではルンゲクッタ4次でローレンツ方程式の計算ができた。そこでとうとうルンゲクッタ8次のDormand-Princeでローレンツ方程式を計算してみよう。

r:
σ:
b:
nmax:

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

プログラムはHTMLのソースを見ればわかるけど、一応ここにも記載。このままCとかにも移せるよ。

var r, sig, bb, n;

function func100420() {
    var ctx = document.getElementById('c100420').getContext('2d');
   
    kx = new Array(14);
    ky = new Array(14);
    kz = new Array(14);
    a = new Array(14);
    c = new Array(14);
    b = new Array(14);
    for (i=0; i <= 14; i++) {
       b[i] = new Array(14);
    }

    a[1] = 0;
    a[2] = 2. / 27.;
    a[3] = 1. / 9.;
    a[4] = 1. / 6.;
    a[5] = 5. / 12.;
    a[6] = 0.5;
    a[7] = 5. / 6.;
    a[8] = 1. / 6.;
    a[9] = 2. / 3.;
    a[10] = 1. / 3.;
    a[11] = 1.;
    a[12] = 0.;
    a[13] = 1.;

    for (j=1; j<=13; j++) {
        for (k=1; k <=13; k++) {
            b[j][k] = 0.;
        }
        kx[j] = 0.;
        ky[j] = 0.;
kz[j] = 0.;
        c[j]=0.;
    }

   
    b[2][1] = 2. / 27.;
   
    b[3][1] = 1. / 36.;
    b[3][2] = 1. / 12.;
   
    b[4][1] = 1. / 24.;
    b[4][3] = 1. / 8.;
   
    b[5][1] = 5. / 12.;
    b[5][3] = -25. / 16.;
    b[5][4] = 25. / 16.;
   
    b[6][1] = 1. / 20.;
    b[6][4] = 1. / 4.;
    b[6][5] = 1. / 5.;
   
    b[7][1] = -25. / 108.;
    b[7][4] = 125. / 108.;
    b[7][5] = -65. / 27.;
    b[7][6] = 125. / 54.;
   
    b[8][1] = 31. / 300.;
    b[8][5] = 61. / 225.;
    b[8][6] = -2. / 9.;
    b[8][7] = 13. / 900.;
   
    b[9][1] = 2.;
    b[9][4] = -53. / 6.;
    b[9][5] = 704. / 45.;
    b[9][6] = -107. / 9.;
    b[9][7] = 67. / 90.;
    b[9][8] = 3.;
   
    b[10][1] = -91. / 108.;
    b[10][4] = 23. / 108.;
    b[10][5] = -976. / 135.;
    b[10][6] = 311. / 54.;
    b[10][7] = -19. / 60.;
    b[10][8] = 17. / 6.;
    b[10][9] = -1. / 12.;
   
    b[11][1] = 2383. / 4100.;
    b[11][4] = -341. / 164.;
    b[11][5] = 4496. / 1025.;
    b[11][6] = -301. / 82.;
    b[11][7] = 2133. / 4100.;
    b[11][8] = 45. / 82.;
    b[11][9] = 45. / 164.;
    b[11][10] = 18. / 41.;
   
    b[12][1] = 3. / 205.;
    b[12][6] = -6. / 41.;
    b[12][7] = -3. / 205.;
    b[12][8] = -3. / 41.;
    b[12][9] = 3. / 41.;
    b[12][10] = 6. / 41.;
   
    b[13][1] = -1777. / 4100.;
    b[13][4] = -341. / 164.;
    b[13][5] = 4496. / 1025.;
    b[13][6] = -289. / 82.;
    b[13][7] = 2193. / 4100.;
    b[13][8] = 51. / 82.;
    b[13][9] = 33. / 164.;
    b[13][10] = 12. / 41.;
    b[13][12] = 1.;
   
    c[6] = 34. / 105.;
    c[7] = 9. / 35.;
    c[8] = 9. / 35.;
    c[9] = 9. / 280.;
    c[10] = 9. / 280.;
    c[12] = 41. / 840.;
    c[13] = 41. / 840.;

    ctx.clearRect(0,0,300,300);
    ctx.strokeStyle ="rgb(0,0,0)";
    ctx.strokeRect(0,0,300,300);

    r = parseFloat(document.Form100420.rvalue.value);
    sig = parseFloat(document.Form100420.svalue.value);
    bb = parseFloat(document.Form100420.bvalue.value);

    n = parseInt(document.Form100420.nvalue.value);

    x=1.e-7;
    y=1.e-7;
    z=1.e-7;
    t=0.;
    dt=0.01;

    ctx.strokeStyle ="rgb(255,0,0)";
    ctx.beginPath();   
    ctx.moveTo(150.+150.*x/30.,150.-150.*y/30.); 
   
    for (i=1;i<=n;i++) {
        kx[1] = fx100420(t + dt, x, y, z) * dt;
        ky[1] = fy100420(t + dt, x, y, z) * dt;
        kz[1] = fz100420(t + dt, x, y, z) * dt;
       
for (j=2;j<=13;j++) {

            xd = x
            yd = y
            zd = z
            for (k=1;k<=j-1;k++) {

                xd = xd + b[j][k] * kx[k];
                yd = yd + b[j][k] * ky[k];
                zd = zd + b[j][k] * kz[k];
            }
            kx[j] = fx100420(t + a[j] * dt, xd, yd, zd) * dt;
            ky[j] = fy100420(t + a[j] * dt, xd, yd, zd) * dt;
            kz[j] = fz100420(t + a[j] * dt, xd, yd, zd) * dt;
        }
       
        for (j=1;j<=13;j++) {
            x = x + c[j] * kx[j];
            y = y + c[j] * ky[j];
            z = z + c[j] * kz[j];
        }

        t = t + dt;
        ctx.lineTo(150.+150.*x/30.,150.-150.*y/30.);   
  }
  ctx.stroke();

}

function fx100420(t,x,y,z) {
     ans = -sig*(x-y);
return (ans);
}
function fy100420(t,x,y,z) {
     ans = -y-x*z+r*x;
return (ans);
}
function fz100420(t,x,y,z) {
     ans =x*y-bb*z;
return (ans);
}