javascriptnumerical-methodsoderunge-kutta

Runge Kutta problems in JS


I'm attempting a Runge-Kutta implementation for a mass on a spring in Javascript and visualizing it with D3. The purpose is to compare it to Forward Euler and comment on the differences. My FE works fine and plots fine, but the Runge-Kutta is shooting off in a negative direction and never wrapping around.

Here's a plunkr with the vis and the code, but I'll add the JS (only for the ODE solvers) too.

// *** 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);

Solution

  • This has some very basic conceptual problems. For a coupled system you have to evaluate all operations simultaneously. That is, in y'(t)=f(y(t)) the function y(t) is vector valued, f has vectors as inputs and vectors as outputs. The Euler method can then be summarized as

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

    allowing flexibility in how the components of f are evaluated. RK4 follows then a similar scheme, the slopes k0,...,k3 are all values of the function f at various modified points.

    The Euler step certainly is not a part of the RK4 step and should also not be mixed up with the system function of the ODE.

    So you should use something in the direction of

    function EulerStep(x,v,h,odefnX,odefnV) {
        var kx = odefnX(x,v);
        var kv = odefnV(x,v);
        return [ x+h*kx, v+h*kv ];
    }
    
    function RK4Step(x,v,h,odefnX,odefnV) {
        var kx0 = odefnX(x,v);
        var kv0 = odefnV(x,v);
    
        var kx1 = odefnX(x+0.5*h*kx0,v+0.5*h*kv0);
        var kv1 = odefnV(x+0.5*h*kx0,v+0.5*h*kv0);
    
        var kx2 = odefnX(x+0.5*h*kx1,v+0.5*h*kv1);
        var kv2 = odefnV(x+0.5*h*kx1,v+0.5*h*kv1);
    
        var kx3 = odefnX(x+    h*kx2,v+    h*kv2);
        var kv3 = odefnV(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, odefnX, odefnV) {
        // 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, odefnX, odefnV); 
        }
        return values;
    }
    
    
    
    // ----------------------------------
    // Example usage
    
    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;
    }
    
    RK4(0.0, 1.0 , 0.01, 100, odefuncX, odefuncV);
    

    See forked Plunker

    // *** Chart variables *** //
    
    var margin = {top: 20, right: 20, bottom: 30, left: 50},
            width = 500 - margin.left - margin.right,
            height = 500 - margin.top - margin.bottom;
    
    var x = d3.scale.linear()
        .range([0, width]);
    
    var y = d3.scale.linear()
        .range([height, 0]);
    
    var xAxis = d3.svg.axis()
        .scale(x)
        .orient("bottom");
    
    var yAxis = d3.svg.axis()
        .scale(y)
        .orient("left");
    
    var line = d3.svg.line()
        .x(function(d) { return x(d.X); })
        .y(function(d) { return y(d.V); });
    
    var odeChart = d3.select("#chart-container").append("svg")
        .attr("width", width + margin.left + margin.right)
        .attr("height", height + margin.top + margin.bottom)
      .append("g")
        .attr("transform", "translate(" + margin.left + "," + margin.top + ")");
    
    // *** Functions for ODE Solvers *** //
    
    function odefunc(y) {
      var x=y[0]; var v=y[1];
        var k = 1; var m = 0.5; var g = 0;
        return [v, (-k/m)*x + g];
    }
    
    function EulerStep(y,h,odefn) {
        var k = odefn(y);
        return [ y[0]+h*k[0], y[1]+h*k[1] ];
    }
    
    function RK4Step(y,h, odefn) {
        var k0 = odefn(y);
     
        var k1 = odefn([y[0]+0.5*h*k0[0],y[1]+0.5*h*k0[1]]);
    
        var k2 = odefn([y[0]+0.5*h*k1[0],y[1]+0.5*h*k1[1]]);
    
        var k3 = odefn([y[0]+    h*k2[0],y[1]+    h*k2[1]]);
    
        return [ y[0]+h/6*(k0[0]+2*(k1[0]+k2[0])+k3[0]),
                 y[1]+h/6*(k0[1]+2*(k1[1]+k2[1])+k3[1]) ];
    }
    
    // FORWARD EULER
    function forwardEuler (x, v, h, n, odefn)
    {
        // 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+1);
    
        // Initial conditions
        values[0] = [x, v];
    
        for (i = 1; i <= n; ++i)
        {
            values[i] = EulerStep(values[i-1], h,odefn);
        }
    
        return values;
    }
    
    // 4TH ORDER RUNGE-KUTTA 
    function RK4 (x, v, h, n, odefn)
    {
        // Initialize an array to hold the values
        var values = new Array(n+1);
    
        // Initial conditions
        values[0] = [x, v];
    
        for (i = 1; i <= n; ++i)
        {
            values[i] = RK4Step(values[i-1], h, odefn);
        }
    
        return values;
    }
    
    // *** Setting up the data *** //
    
    var rkValues = RK4(1, 0, 0.1, 100, odefunc);
    var feValues = forwardEuler(1, 0, 0.025, 400, odefunc);
    
    var rkData = rkValues.map(function(d) {
      return {
         X: +d[0],
         V: +d[1]
      };});
    
    var feData = feValues.map(function(d) {
      return {
         X: +d[0],
         V: +d[1]
      };});
    
    x.domain(d3.extent(feData, function(d) { return d.X; }));
    y.domain(d3.extent(feData, function(d) { return d.V; }));
    
    odeChart.append("g")
      .attr("class", "x axis")
      .attr("transform", "translate(0," + height + ")")
      .call(xAxis);
    
    odeChart.append("g")
      .attr("class", "y axis")
      .call(yAxis)
    
    odeChart.append("path")
      .datum(rkData)
      .attr("class", "rkLine")
      .attr("d", line);
    
    odeChart.append("path")
      .datum(feData)
      .attr("class", "feLine")
      .attr("d", line);
    body {
        font: 10px sans-serif;
      }
    
      #chart-container {
        width: 600px;
        margin: 0 auto;
      }
    
      .axis path,
      .axis line {
        fill: none;
        stroke: #000;
        shape-rendering: crispEdges;
      }
    
      .x.axis path {
        display: none;
      }
    
      .rkLine {
        fill: none;
        stroke: steelblue;
        stroke-width: 1.5px;
      }
    
      .feLine {
        fill: none;
        stroke: red;
        stroke-width: 1.5px;
      }
    <script src="https://cdnjs.cloudflare.com/ajax/libs/d3/3.5.17/d3.min.js"></script>
    <div id="chart-container"> <h1>Runge Kutta Order 4</h1 ></div>