JS 中的 Runge Kutta 问题

Question

我正在尝试使用 Javascript 对弹簧上的质量进行 Runge-Kutta 实现，并使用 D3 将其可视化。目的是将其与 Forward Euler 进行比较并评论差异。我的 FE 工作正常，情节也很好，但是 Runge-Kutta 正在向负面方向射击，并且从不环绕。

这是一个 plunkr 的 vis 和代码，但我也会添加 JS（仅适用于 ODE 求解器）。

// *** Functions for ODE Solvers *** //

function FEx (x, v, h)
{
    return x + h*v;
}

function FEv (x, v, h)
{
    var k = 1; var m = 0.5; var g = 0;

    return v + h*( (-k/m)*x + g );
}

function RKx (x, v, h)
{
    var k1 = FEx(x, v, h);
    var k2 = FEx(x+h/2*k1, v+h/2, h);
    var k3 = FEx(x+h/2*k2, v+h/2, h);
    var k4 = FEx(x+h*k3, v+h, h);

    return x + h/6*(k1 + 2*k2 + 2*k3 + k4);
}

function RKy (x, v, h)
{
    var k1 = FEv(x, v, h);
    var k2 = FEv(x+h/2, v+h/2*k1, h);
    var k3 = FEv(x+h/2, v+h/2*k2, h);
    var k4 = FEv(x+h, v+h*k3, h);

    return v + h/6*(k1 + 2*k2 + 2*k3 + k4);
}

// FORWARD EULER
function forewardEuler (x, v, h, n)
{
    // Initialize an array to hold the values
    // JS doesn't really support multi-dimensional arrays
    // so this is a "jagged" nested array
    var values = new Array(n);
    for(i = 0; i < values.length; i++)
        values[i] = new Array(2);

    // Initial conditions
    values[0] = [x, v];

    for (i = 1; i < n; ++i)
    {
        values[i][0] = FEx(values[i-1][0], values[i-1][1], h);
        values[i][1] = FEv(values[i-1][0], values[i-1][1], h);
    }

    return values;
}

// 4TH ORDER RUNGE-KUTTA 
function RK4 (x, v, h, n)
{
    // Initialize an array to hold the values
    var values = new Array(n);
    for(i = 0; i < values.length; i++)
        values[i] = new Array(2);

    // Initial conditions
    values[0] = [x, v];

    for (i = 1; i < n; ++i)
    {
        values[i][0] = RKx(values[i-1][0], values[i-1][1], h);
        values[i][1] = RKy(values[i-1][0], values[i-1][1], h);
    }

    return values;
}

// *** Setting up the data *** //

var rkValues = RK4(1, 0, 0.1, 100);
var feValues = forewardEuler(1, 0, 0.1, 100);

Answer 1

这有一些非常基本的概念问题。对于耦合系统，您必须同时评估所有操作。也就是说，在y'(t)=f(y(t)) 中，函数y(t) 是向量值，f 将向量作为输入，将向量作为输出。欧拉方法可以总结为

k = f(y[i]);
y[i+1] = y[i] + h*k;

允许灵活地评估f 的组件。然后RK4遵循类似的方案，斜率k0,...,k3都是函数f在各个修改点的值。

Euler 步当然不是 RK4 步的一部分，也不应该与 ODE 的系统函数混淆。

所以你应该使用方向的东西

function odefuncX(x,v) {return v;}

function odefuncV(x,v) { 
    var k = 1; var m = 0.5; var g = 0;
    return (-k/m)*x + g;
}

function EulerStep(x,v,h) {
    var kx = odefuncX(x,v);
    var kv = odefuncV(x,v);
    return [ x+h*kx, v+h*kv ];
}

function RK4Step(x,v,h) {
    var kx0 = odefuncX(x,v);
    var kv0 = odefuncV(x,v);

    var kx1 = odefuncX(x+0.5*h*kx0,v+0.5*h*kv0);
    var kv1 = odefuncV(x+0.5*h*kx0,v+0.5*h*kv0);

    var kx2 = odefuncX(x+0.5*h*kx1,v+0.5*h*kv1);
    var kv2 = odefuncV(x+0.5*h*kx1,v+0.5*h*kv1);

    var kx3 = odefuncX(x+    h*kx2,v+    h*kv2);
    var kv3 = odefuncV(x+    h*kx2,v+    h*kv2);

    return [ x+h/6*(kx0+2*(kx1+kx2)+kx3),
             v+h/6*(kv0+2*(kv1+kv2)+kv3) ];
}

// 4TH ORDER RUNGE-KUTTA 
function RK4 (x, v, h, n) {
    // Initialize an array to hold the values
    var values = new Array(n);

    // Initial conditions
    values[0] = [x, v];
    for (i = 1; i < n; ++i) {
        values[i] = RK4Step(values[i-1][0], values[i-1][1], h); 
    }
    return values;
}

见forked Plunker