Archive: rk4-test.tar.gz



// Simple (not optimized) implementation of Runge-Kutta 4
// Needs linking with math library (GCC option -lm)
// Copyright: (c) Petr Peringer, FIT TU Brno, 2012

#include <stdio.h>
#include <math.h>               // sin() - used only for error computation

// System: circle test
#define SIZE 2                  // state vector size of the system
// "Dynamic section" = continuous system description
// Input:  t = time of step start,
//         y = state vector = integrator outputs
// Output: yin = vector of derivatives = integrators input values
void f(double t, double *y, double *yin)
{
    yin[0] = y[1];
    yin[1] = -y[0];
}

// Runge-Kutta 4-th order method step
// Parameters: t = start time,
//             h = step size,
//             y = state vector
// Output:  new state in y
// Warning: it does not update time
void RK4step(double t, double h, double *y)
{
    double yin[SIZE];           // integrator input = derivative
    double ystart[SIZE];        // initial state
    // 4 stages results:
    double k1[SIZE];
    double k2[SIZE];
    double k3[SIZE];
    double k4[SIZE];

    int i;
    for (i = 0; i < SIZE; i++)  // save initial value
        ystart[i] = y[i];

    f(t, y, yin);               // yin = f(t, y(t))

    for (i = 0; i < SIZE; i++) {
        k1[i] = h * yin[i];     // k1 = h * f(t, y(t))
        y[i] = ystart[i] + k1[i] / 2;
    }

    f(t + h / 2, y, yin);       // yin = f(t+h/2, y(t)+k1/2)

    for (i = 0; i < SIZE; i++) {
        k2[i] = h * yin[i];     // k2 = h * f(t+h/2, y(t)+k1/2)
        y[i] = ystart[i] + k2[i] / 2;
    }

    f(t + h / 2, y, yin);       // yin = f(t+h/2, y(t)+k2/2)

    for (i = 0; i < SIZE; i++) {
        k3[i] = h * yin[i];     // k3 = h * f(t+h/2, y(t)+k2/2)
        y[i] = ystart[i] + k3[i];
    }

    f(t + h, y, yin);           // yin = f(t+h, y(t)+k3)

    for (i = 0; i < SIZE; i++) {
        k4[i] = h * yin[i];     // k4 = h * f(t+h, y(t)+k3)
        // Result:
        // y(t+h) = y(t) + k1/6 + k2/3 + k3/3 + k4/6;
        y[i] = ystart[i] + k1[i] / 6 + k2[i] / 3 + k3[i] / 3 + k4[i] / 6;
    }
    // Return: y = new state at time t+h
}

// Experiment control
int main()
{
    double y[SIZE] = { 0, 1 };
    double stepsize = 0.1;

    double t = 0;
    while (t < 20) {
        printf("%8.3f % -11g % -11g  % e\n", t, y[0], y[1], y[0] - sin(t));
        RK4step(t, stepsize, y);        // make step
        t += stepsize;          // advance the simulation time
    }
    printf("%8.3f % -11g % -11g  % e\n", t, y[0], y[1], y[0] - sin(t));
    return 0;
}