javascript - 翻译后的 JavaScript 中的舍入错误

Question

我翻译了一个旧的 Fortran 程序，该程序模拟阻尼驱动线性振荡器。这是 JS:

/*jslint browser: true, devel: true */
/*global VERLET */
/*global FORCE */
/*global $ */
function driven() {
  "use strict";
//MODULE shared
  //USE, INTRINSIC :: ISO_FORTRAN_ENV,dp=>REAL64!modern DOUBLE PRECISION
  //INTEGER :: i
  //INTEGER, PARAMETER :: i_max=5000
  //REAL(dp) :: x_read,v_read,const0,gamma,A_o,dt
  //REAL(dp), DIMENSION(:), ALLOCATABLE :: x,v,a,E,dE,time
  //REAL(dp), PARAMETER :: m=1.0
  var i_max,
    i_max_read,
    i,
    x_read,
    v_read,
    const0,
    const0_read,
    gamma,
    gamma_read,
    A_o,
    A_o_read,
    dt,
    dt_read,
    x,
    x_plot,
    v,
    v_plot,
    a,
    a_plot,
    time,
    f_osc,
    sim,
    result,
    nonlinear;
  //i_max = 5000;
  x = [];
  x_plot = [];//new combined array for Flot plotting
  v = [];
  v_plot = [];//new combined array for Flot plotting
  a = [];
  a_plot = [];//new combined array for Flot plotting
  time = [];
  f_osc = [];//unlike Fortran, this needs to be an array
//END MODULE shared

//PROGRAM nonlinear
  nonlinear = document.forms.nonlinear;
  //USE shared
  //IMPLICIT NONE
  //EXTERNAL VERLET
  //ALLOCATE(x(i_max + 1))
  //ALLOCATE(v(i_max + 1))
  //ALLOCATE(a(i_max + 1))
  //ALLOCATE(E(i_max + 1))
  //ALLOCATE(dE(i_max + 1))
  //ALLOCATE(time(i_max + 1))
  //PRINT *, ' '
  //PRINT *, 'Initial position of the mass? [m]'
  //PRINT *, ' '
  //READ *, x_read
  x_read = nonlinear.elements.x_read;
  //x(1) = x_read
  x[0] = parseFloat(x_read.value);
  //PRINT *, ' '
  //PRINT *, 'Initial velocity of the mass? [m/sec]'
  //PRINT *, ' '
  //READ *, v_read
  v_read = nonlinear.elements.v_read;
  //v(1) = v_read
  v[0] = parseFloat(v_read.value);
  //PRINT *, ' '
  //PRINT *, 'Value of k/m? [1/sec^2]'
  //PRINT *, ' '
  //READ *, const0
  const0_read = nonlinear.elements.const0_read;
  const0 = parseFloat(const0_read.value);
  //PRINT *, ' '
  //PRINT *, 'Value of the damping coefficient? [kg/sec]'
  //PRINT *, ' '
  //READ *, gamma
  gamma_read = nonlinear.elements.gamma_read;
  gamma = parseFloat(gamma_read.value);
  //PRINT *, ' '
  //PRINT *, 'Amplitude of the external force? [N]'
  //PRINT *, ' '
  //READ *, A_o
  A_o_read = nonlinear.elements.A_o_read;
  A_o = parseFloat(A_o_read.value);
  i_max_read = nonlinear.elements.i_max_read;
  i_max = parseFloat(i_max_read.value);
  //PRINT *, ' '
  //PRINT *, 'Time-step of the system? [sec]'
  //PRINT *, ' '
  //READ *, dt
  //PRINT *, ' '
  dt_read = nonlinear.elements.dt_read;
  dt = parseFloat(dt_read.value);
  sim = i_max * dt;
  result = sim.toFixed(1);
  document.getElementById('output').innerHTML = result;
  //time(1) = 0.0
  time[0] = 0.0;
  //i = 1
  i = 0;
  //DO
  do {
    //IF(i > i_max) EXIT
    //CALL VERLET
    VERLET(i, x, v, a, time, f_osc, dt, A_o, const0, gamma);
    x_plot[i] = [time[i], x[i]];
    v_plot[i] = [time[i], v[i]];
    a_plot[i] = [time[i], a[i]];
    //time(i + 1) = time(i) + dt
    time[i + 1] = time[i] + dt;
    //i = i + 1
    i = i + 1;
  //END DO
  } while (i < i_max);
  //OPEN(7, file='xt.dat', status='unknown')
  //WRITE(7,'(2f16.6)') (time(i),x(i),i=1,i_max)
  //CLOSE(7)
  //DEALLOCATE(x)
  //DEALLOCATE(v)
  //DEALLOCATE(a)
  //DEALLOCATE(E)
  //DEALLOCATE(dE)
  //DEALLOCATE(time)
  function doPlot(position) {//Flot
    $.plot("#placeholder", [//data
      {
        data: x_plot,
        label: "Position (m)",
        yaxis: 1,
        color: "red"
      },
      {
        data: v_plot,
        label: "Velocity (m/sec)",
        yaxis: 2,
        color: "green"
      },
      {
        data: a_plot,
        label: "Acceleration (m/sec/sec)",
        yaxis: 3,
        color: "blue"
      }
    ],//options
      {
        xaxis: { axisLabel: "Time (sec)" },
        yaxes: [
          { font: { color: "red" } },
          { font: { color: "green" } },
          { font: { color: "blue" } },
          { alignTicksWithAxis: position === "left" ? 1 : null }
        ],
        legend: {
          position: "nw",
          labelBoxBorderColor: null
        }
      }
      );
  }
  doPlot("left");
//END PROGRAM nonlinear
}

//SUBROUTINE VERLET()
function VERLET(i, x_temp, v_temp, a_temp, t_temp, f_osc_temp, dt, A_o, const0, gamma) {
  "use strict";
  //USE shared
  //IMPLICIT NONE
  //REAL(dp) :: x_temp,v_temp,t_temp,f_osc
  //EXTERNAL FORCE,ENERGY
  var m;
  m = 1.0;
  //x_temp = x(i)
  //v_temp = v(i)
  //t_temp = time(i)
  //CALL FORCE(x_temp,v_temp,t_temp,f_osc)
  FORCE(i, x_temp, v_temp, t_temp, f_osc_temp, A_o, const0, gamma);
  //x_temp = x(i)
  //v_temp = v(i)
  //CALL ENERGY!don't think I'm actually using this; no INTENT(OUT)
  //a(i) = f_osc/m
  a_temp[i] = f_osc_temp[i] / m;
  //x(i + 1) = x(i) + v(i)*dt + 0.5*a(i)*dt*dt
  x_temp[i + 1] = x_temp[i] + v_temp[i] * dt + 0.5 * a_temp[i] * dt * dt;
  //x_temp = x(i + 1)
  //v_temp = v(i)
  //t_temp = time(i + 1)
  //CALL FORCE(x_temp,v_temp,t_temp,f_osc)
  FORCE(i, x_temp, v_temp, t_temp, f_osc_temp, A_o, const0, gamma);
  //a(i + 1) = f_osc/m
  a_temp[i + 1] = f_osc_temp[i] / m;
  //v(i + 1) = v(i) + 0.5*(a(i + 1) + a(i))*dt
  v_temp[i + 1] = v_temp[i] + 0.5 * (a_temp[i + 1] + a_temp[i]) * dt;
  //x_temp = x(i + 1)
  //v_temp = v(i + 1)
  //CALL ENERGY!don't think I'm actually using this; no INTENT(OUT)
  return [x_temp, v_temp, a_temp, t_temp];
//END
}

//SUBROUTINE FORCE(xs,vs,ts,f_oscs)
function FORCE(i, xs, vs, ts, f_oscs, A_o, const0, gamma) {
  "use strict";
  //USE shared
  //IMPLICIT NONE
  //REAL(dp), INTENT(IN) :: xs,vs,ts
  //REAL(dp), INTENT(OUT) :: f_oscs
  //REAL(dp) :: f_t
  var f_t;
  //f_t = A_o*DCOS(2.0*ts)
  f_t = A_o * Math.cos(2.0 * ts[i]);
  //f_oscs = -const0*xs - gamma*vs + f_t
  f_oscs[i] = -const0 * xs[i] - gamma * vs[i] + f_t;
  return f_oscs;
//END
}

//SUBROUTINE ENERGY()
  //USE shared
  //IMPLICIT NONE
  //!REAL(dp), INTENT(OUT) ::
  //E(i) = 0.5*(v(i)**2 + const0*x(i)**2)
  //E(1) = 0.5*(v(1)**2 + const0*x(1)**2)
  //dE(i) = E(i) - E(1)
//END

原始 Fortran 源代码显示在注释中。比较原始 Fortran 程序的位置与时间的输出:

    0.000000        0.100000
    0.010000        0.129930
    0.020000        0.159693
    ...

...来自翻译后的 JavaScript:

    0.000000        0.100000
    0.010000        0.129930
    0.020000        0.159707
    ...
    

结果从那时开始变得更糟。 5000 步后，结果相差 0.01。我不习惯在 JS 中使用数组，所以错误可能出现在某个地方。这是我第一次遇到这样的事情。有人知道这里发生了什么吗？

<小时/>

更新 1

JS 本身:

/*jslint browser: true, devel: true */
/*global VERLET */
/*global FORCE */
/*global $ */
function driven() {
  "use strict";
  var i_max,
    i_max_read,
    i,
    x_read,
    v_read,
    const0,
    const0_read,
    gamma,
    gamma_read,
    A_o,
    A_o_read,
    dt,
    dt_read,
    x,
    x_plot,
    v,
    v_plot,
    a,
    a_plot,
    time,
    f_osc,
    sim,
    result,
    nonlinear;
  x = [];
  v = [];
  a = [];
  time = [];
  f_osc = [];//unlike Fortran, this needs to be an array
  time[0] = 0.0;
  i = 0;
  do {
    VERLET(i, x, v, a, time, f_osc, dt, A_o, const0, gamma);
    time[i + 1] = time[i] + dt;
    i = i + 1;
  } while (i < i_max);
}

function VERLET(i, x_temp, v_temp, a_temp, t_temp, f_osc_temp, dt, A_o, const0, gamma) {
  "use strict";
  var m;
  m = 1.0;
  FORCE(i, x_temp, v_temp, t_temp, f_osc_temp, A_o, const0, gamma);
  a_temp[i] = f_osc_temp[i] / m;
  x_temp[i + 1] = x_temp[i] + v_temp[i] * dt + 0.5 * a_temp[i] * dt * dt;
  FORCE(i, x_temp, v_temp, t_temp, f_osc_temp, A_o, const0, gamma);
  a_temp[i + 1] = f_osc_temp[i] / m;
  v_temp[i + 1] = v_temp[i] + 0.5 * (a_temp[i + 1] + a_temp[i]) * dt;
  return [x_temp, v_temp, a_temp, t_temp];
}

function FORCE(i, xs, vs, ts, f_oscs, A_o, const0, gamma) {
  "use strict";
  var f_t;
  f_t = A_o * Math.cos(2.0 * ts[i]);
  f_oscs[i] = -const0 * xs[i] - gamma * vs[i] + f_t;
  return f_oscs;
}

<小时/>

更新2

这是一个类似的 JS，没有表现出这种症状:

/*jslint browser: true, devel: true */
/*global $ */
function amp() {
  "use strict";
  var i_max,
    i,
    m,
    f_osc,
    theta,
    theta_plot,
    omega,
    omega_plot,
    a,
    e,
    e_plot,
    de,
    time,
    theta0,
    read_ratio,
    const0,
    T,
    result,
    omega_read,
    dt_read,
    dt,
    amplitude;
  i_max = 500;
  m = 1.0;
  f_osc = 0.0;
  theta = [];
  omega = [];
  a = [];
  e = [];
  de = [];
  time = [];
  e[0] = 0.5 * (Math.pow(omega[0], 2) + const0 * Math.pow(theta[0], 2));
  time[0] = 0.0;
  i = 0;
  do {
    f_osc = -const0 * Math.sin(theta[i]);
    a[i] = f_osc / m;
    e[i] = 0.5 * (Math.pow(omega[i], 2) + const0 * Math.pow(theta[i], 2));
    de[i] = e[i] - e[0];
    theta[i + 1] = theta[i] + omega[i] * dt + 0.5 * a[i] * dt * dt;
    f_osc = -const0 * Math.sin(theta[i + 1]);
    a[i + 1] = f_osc / m;
    omega[i + 1] = omega[i] + 0.5 * (a[i + 1] + a[i]) * dt;
    e[i] = 0.5 * (Math.pow(omega[i + 1], 2) + const0 * Math.pow(theta[i + 1], 2));
    de[i] = e[i] - e[0];
    time[i + 1] = time[i] + dt;
    i = i + 1;
  } while (i < i_max);
}

最大的区别是这个函数不包含我视为子例程的函数。我在这方面是否做了违法的事情？

Answer 1

第二次调用 FORCE() 所做的事情与 FORTRAN 对应的操作不同......

在 FORTRAN 中，你可以这样做:

a(i) = FORCE( x(i) , v(i) , time(i) ) / m
a(i+1) = FORCE( x(i) + v(i)*dt + 0.5*a(i)*dt*dt , v(i) , time(i+1) ) / m

在 JS 中，你做了另一件事:

a(i) = FORCE( x(i) , v(i) , time(i) ) / m
a(i+1) = FORCE( x(i) , v(i) , time(i) ) / m

您在处理索引方面遇到了真正的困难，因为 FORTRAN a(i+1) 的中间结果需要混合 x(i+1) 和 v(i)...
因此，仅传递 i 作为参数听起来不太合适
您可以设置 v[i+1]=v[i] 并尝试 FORCE(i+1,...) 但这听起来太棘手了。

要么坚持更接近的翻译，要么尝试完全重写。