I’m learning C by myself using real-world examples. Given my background (physics) I’m implementing a 4th order Runge-Kutta method. I have a function which performs a single RK4 step with arbitrary input ODE system and timestep. When calculating intermediate steps for RK4, I have two approaches, but I can’t decide which one is more efficient.

In the 1st approach, I’m allocating an array ‘k’ on the heap, passing it to the RK4 step function, and doing calculations. You can see the related functions below:

void rk4Step(double t, double dt, double* p, double* y, double* k, int n_eqns, ODE* f) {
    double y_orig_val;
    for (int i = 0; i < n_eqns; i++) {
        y_orig_val = y[i];
        k[0] = dt * (*f)(p, t, y, i);

        y[i] += 0.5 * k[0];
        k[1] = dt * (*f)(p, t + 0.5 * dt, y, i);
        y[i] = y_orig_val;

        y[i] += 0.5 * k[1];
        k[2] = dt * (*f)(p, t + 0.5 * dt, y, i);
        y[i] = y_orig_val;

        y[i] += k[2];
        k[3] = dt * (*f)(p, t + dt, y, i);
        y[i] = y_orig_val;

        y[i] += 1.0 / 3.0 * (1.0 / 2.0 * (k[0] + k[3]) + (k[1] + k[2]));
    }
}

void solveODE(double t, double dt, int n_eqns, int n_params, int n_steps, ODE* f) {

    double* y = malloc(n_eqns * sizeof(double));
    double* p = malloc(n_params * sizeof(double));
    double* k = malloc(4 * sizeof(double));
    init(p, y, n_eqns);

    FILE* sol = fopen("sol.dat", "w");

    for (int i = 0; i <= n_steps; i++) {
        fprintf(sol, "%f\t", t);
        for (int j = 0; j < n_eqns-1; j++) {
            fprintf(sol, "%f\t", y[j]);
        }
        fprintf(sol, "%f\n", y[n_eqns-1]);

        rk4Step(t, dt, p, y, k, n_eqns, f);
        t += dt;
    }

    free(y);
    free(p);
    free(k);
    fclose(sol);
}

In the other approach the variable ‘k’ is not allocated on the heap, but on the stack in the body of the rk4Step function. As you can see, the rk4Step function is used in a loop. This was the variable k would be allocated on the stack again and again in each iteration.

I know that working with variables on the stack is much faster, but is it true in this case too?

>Solution :

No need to allocate on the heap:

void solveODE(double t, double dt, int n_eqns, int n_params, int n_steps, ODE* f) {

    double y[n_eqns];
    double p[n_params];
    double k[4];
    init(p, y, n_eqns);

    FILE* sol = fopen("sol.dat", "w");

    for (int i = 0; i <= n_steps; i++) {
        fprintf(sol, "%f\t", t);
        for (int j = 0; j < n_eqns-1; j++) {
            fprintf(sol, "%f\t", y[j]);
        }
        fprintf(sol, "%f\n", y[n_eqns-1]);

        rk4Step(t, dt, &p, &y, &k, n_eqns, f);
        t += dt;
    }

 
    fclose(sol);
}