Numerical integration is like magic

gistfile1.txt

#include <stdio.h>

void integrate_euler(float t0, float dt, float x0, float v0, float *x, float *v, float(*accelerate)(float,float,float))
{
    float a = accelerate(t0, x0, v0);
    float dx = v0 * dt + 0.5 * a * dt * dt;
    float dv = a * dt;

    *x = x0 + dx;
    *v = v0 + dv;
    //update();
}

void integrate_rkn(float t0, float dt, float x0, float v0, float *x, float *v, float(*accelerate)(float,float,float))
{
    float a0 = accelerate(t0, x0, v0);
    float k1 = 0.5 * dt * a0;
    float k = 0.5 * dt * (v0 + 0.5 * k1);
    // update();
    float a1 = accelerate(t0 + 0.5 * dt, x0 + k, v0 + k1);
    float k2 = 0.5 * dt * a1;
    float a2 = accelerate(t0 + 0.5 * dt, x0 + k, v0 + k2);
    float k3 = 0.5 * dt * a2;
    float l = dt * (v0 + k3);
    // update();
    float a3 = accelerate(t0 + dt, x0 + l, v0 + 2.0 * k3);
    float k4 = 0.5 * dt * a3;

    float dx = dt * (v0 + (k1 + k2 + k3)/3.0);
    float dv = (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 3.0;
    // update()

    *x = x0 + dx;
    *v = v0 + dv;
}

void integrate_2step(float t0, float dt, float x0, float v0, float *x, float *v, float(*accelerate)(float,float,float))
{
    // NOTE: this seems to be better than euler but still not very good
    float a0 = accelerate(t0, x0, v0);
    float k1 = 0.5 * dt * a0;
    float k = 0.5 * dt * (v0 + k1);
    // update();
    float a1 = accelerate(t0 + 0.5 * dt, x0 + k, v0 + k1);
    float k2 = 0.5 * dt * a1;

    //printf(
        //"a0: %2.2f\t k1: %2.2f\t k: %2.2f\t a1: %2.2f\t k2: %2.2f\n",
        //a0, k1, k,  a1, k2);

    float dx = dt * (v0 + (k1 + 3.0 * k2)/4.0);
    float dv = (k1 + 3.0 * k2) / 2.0;
    // update()

    *x = x0 + dx;
    *v = v0 + dv;
}

#define abs(x) ((x) < 0 ? -(x) : (x))

float acc(float t, float x, float v)
{
    (void)t; (void)x; (void)v;
    return -x;
    //return -x - 0.1 * v;
    //return 1;
}


int main()
{
    float dt = 1.0/120.0;
    float t0 = 0, t1 = 100 + dt / 2;

    float x0 = 1, v0 = 0;

    for(float t = t0; t <= t1; t += dt)
    {
        float x, v;
        integrate_rkn(t, dt, x0, v0, &x, &v, acc);
        printf("%04.4f\t%03.8f\t%03.8f\t%03.8f\n", t+dt, x, v, acc(t, x, v));
        x0 = x; v0 = v;
    }

    return 0;
}