Проблемы с Runge Kutta в JS

Question

Я пытаюсь реализовать реализацию Runge-Kutta для массы на пружине в Javascript и визуализировать ее с помощью D3. Цель состоит в том, чтобы сравнить его с Forward Euler и прокомментировать различия. Моя FE работает отлично, и сюжеты прекрасны, но Рунге-Кутта стреляет в отрицательном направлении и никогда не обматывает.

Вот plunkr с изображением и кодом, но я добавлю 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 в различных модифицированных точках.

Шаг Эйлера, конечно же, не является частью шага 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