Runge-Kutta + Lorenz attractor in C and JavaScript

rk4.c

#include <stdlib.h>
#include "rk4.h"

static struct vector *a;
static struct vector *b;
static struct vector *c;
static struct vector *d;
static struct vector *tmp1;
static struct vector *tmp2;
static struct vector *tmp3;
static struct vector *tmp4;
static struct vector *tmp5;
static struct vector *tmp6;
static struct vector *tmp7;
static struct vector *tmp8;
static struct vector *tmp9;
static struct vector *tmp10;
static struct vector *tmp11;

struct vector *create_vector(unsigned int size) {
    struct vector *v = malloc(sizeof(struct vector));
    v->elts = malloc(size * sizeof(double));
    v->length = size;
    return v;
}

static struct vector *vec_add(struct vector *v1, struct vector *v2, struct vector *out) {
    unsigned int i = v1->length;
    while (i--) {
        out->elts[i] = v1->elts[i] + v2->elts[i];
    }
    return out;
}

static struct vector *vec_mul(double s, struct vector *v, struct vector *out) {
    unsigned int i = v->length;
    while (i--) {
        out->elts[i] = s * v->elts[i];
    }
    return out;
}

static struct vector *vec_eval(vf_t *fs, struct vector *vs, struct vector *out) {
    unsigned int i = vs->length;
    while (i--) {
        out->elts[i] = (*fs[i])(vs);
    }
    return out;
}

void rk4_step(double h, vf_t *fs, struct vector *ps, struct vector *out) {
    vec_eval(fs, ps, a);
    vec_eval(fs, vec_add(ps, vec_mul(h/2, a, tmp1), tmp2), b);
    vec_eval(fs, vec_add(ps, vec_mul(h/2, b, tmp3), tmp4), c);
    vec_eval(fs, vec_add(ps, vec_mul(h, c, tmp5), tmp6), d);
    vec_add(ps, vec_mul(h/6, vec_add(vec_add(a, d, tmp7), vec_mul(2, vec_add(b, c, tmp8), tmp9), tmp10), tmp11), out);
}

void rk4_init(unsigned int size) {
    a = create_vector(size);
    b = create_vector(size);
    c = create_vector(size);
    d = create_vector(size);
    tmp1 = create_vector(size);
    tmp2 = create_vector(size);
    tmp3 = create_vector(size);
    tmp4 = create_vector(size);
    tmp5 = create_vector(size);
    tmp6 = create_vector(size);
    tmp7 = create_vector(size);
    tmp8 = create_vector(size);
    tmp9 = create_vector(size);
    tmp10 = create_vector(size);
    tmp11 = create_vector(size);
}

rk4.h

struct vector {
    unsigned int length;
    double *elts;
};

typedef double (*vf_t) (const struct vector *);

void rk4_init(unsigned int size);
struct vector *create_vector(unsigned int size);
void rk4_step(double h, vf_t *fs, struct vector *ps, struct vector *out);

rk4.js

(function() {
    /* Simple 4th-order Runge Kutta method */

    /**
     * Add two vectors, put result in `out'
     * @param {Array.<number>} v1
     * @param {Array.<number>} v2
     * @return {Array.<number>}
     */
    function vec_add(v1, v2, out) {
        var i = v1.length;
        while (i--) {
            out[i] = v1[i] + v2[i];
        }
        return out;
    }

    /**
     * Multiply a vector by a scalar, put result in `out'
     * @param {number} s
     * @param {Array.<number>} v
     * @return {Array.<number>}
     */
    function vec_mul(s, v, out) {
        var i = v.length;
        while (i--) {
            out[i] = s * v[i];
        }
        return out;
    }

    /**
     * fs: Array of functions taking a vector and returning a double
     * vs: Array of vectors
     * For each function in fs, call it with the vs and put result
     * in corresponding element of `out'
     * @param {Array.<function(...[number]): number>} fs
     * @param {Array.<number>} vs
     * @return {Array.<number>}
     */
    function vec_eval(fs, vs, out) {
        var i = fs.length;
        while (i--) {
            out[i] = fs[i].call(null, vs);
        }
        return out;
    }

    var a, b, c, d;
    var tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, tmp10, tmp11;

    function init(size) {
        a = new Float64Array(size);
        b = new Float64Array(size);
        c = new Float64Array(size);
        d = new Float64Array(size);
        tmp1 = new Float64Array(size);
        tmp2 = new Float64Array(size);
        tmp3 = new Float64Array(size);
        tmp4 = new Float64Array(size);
        tmp5 = new Float64Array(size);
        tmp6 = new Float64Array(size);
        tmp7 = new Float64Array(size);
        tmp8 = new Float64Array(size);
        tmp9 = new Float64Array(size);
        tmp10 = new Float64Array(size);
        tmp11 = new Float64Array(size);
    }
    /**
     * Performs a single Runge-Kutta step
     * @param {number} h time step
     * @param {Array.<function(...[number]): number>} fs Array of time derivative functions
     * @param {Array.<number>} ps positions
     * @param {Array.<number>} out output array
     */
    function rk4(h, fs, ps, out) {
        vec_eval(fs, ps, a);
        vec_eval(fs, vec_add(ps, vec_mul(h/2, a, tmp1), tmp2), b);
        vec_eval(fs, vec_add(ps, vec_mul(h/2, b, tmp3), tmp4), c);
        vec_eval(fs, vec_add(ps, vec_mul(h, c, tmp5), tmp6), d);
        vec_add(ps, vec_mul(h/6, vec_add(vec_add(a, d, tmp7), vec_mul(2, vec_add(b, c, tmp8), tmp9), tmp10), tmp11), out);
    }
    exports.step = rk4;
    exports.init = init;
})();

strange.c

#include "rk4.h"
#include <sys/time.h>
#include <stdio.h>

static const unsigned int size = 3;
static const double sigma = 10;
static const double r = 28;
static const double b = 8.0/3.0;


static double dx(const struct vector *v) {
    return sigma * (v->elts[1] - v->elts[0]);
}

static double dy(const struct vector *v) {
    return r*v->elts[0] - v->elts[1] - v->elts[0]*v->elts[2];
}

static double dz(const struct vector *v) {
    return v->elts[0]*v->elts[1] - b*v->elts[2];
}


int main(void) {
    struct timeval tim;
    double t1, t2;
    struct vector *ps = create_vector(size);
    ps->elts[0] = 10;
    ps->elts[1] = -10;
    ps->elts[2] = 20;

    vf_t fs[3] = {&dx, &dy, &dz};
    int i;
    double h = 0.005;
    rk4_init(3);
    gettimeofday(&tim, NULL);
    t1 = tim.tv_sec+(tim.tv_usec/1000000.0);
    for (i = 0; i < 1000000; i++) {
        rk4_step(h, fs, ps, ps);
    }
    gettimeofday(&tim, NULL);
    t2 = tim.tv_sec+(tim.tv_usec/1000000.0);
    printf("%.6lf\n", t2-t1);
    return 0;
}

strange.js

(function() {
    var rk4 = require('./rk4');

    var sigma = 10;
    var r = 28;
    var b = 8/3;

    function dx(v) {
        return sigma * (v[1] - v[0]);
    }

    function dy(v) {
        return r*v[0] - v[1] - v[0]*v[2];
    }

    function dz(v) {
        return v[0]*v[1] - b*v[2];
    }

    /**
    * time derivative functions
    */
    var fs = [dx, dy, dz];

    /**
    * initial positions
    */
    var ps = new Float64Array([10, -10, 20]);
    var h = 0.005;

    rk4.init(ps.length);

    var start = (new Date()).getTime();
    var i = 1000000;
    while(i--) {
        rk4.step(h, fs, ps, ps);
    }
    var end = (new Date()).getTime();
    console.log(((end - start) / 1000));
})();