Archive: rk4-test.tar.gz

//////////////////////////////////////////////////////////////////
// Simple (not optimized) implementation of Runge-Kutta 4
// used for solving of initial value problem (Circle test)
// Needs linking with math library (GCC option: -lm)
// Copyright: (c) Petr Peringer, FIT TU Brno, 2012, 2022

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

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

//////////////////////////////////////////////////////////////////
// Runge-Kutta 4-th order method
// Parameters: t = step 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 der[SIZE];           // integrator inputs = derivatives
    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, der);               // der = f(t, y(t))

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

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

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

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

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

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

    for (i = 0; i < SIZE; i++) {
        k4[i] = h * der[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
}

//////////////////////////////////////////////////////////////////
// Simulation experiment control
int main()
{
    double y[SIZE] = { 0, 1 };  // initial values:  y'(0)=1, y(0)=0
    double stepsize = 0.1;

    printf("# Circle test solution (RK4, step=%g)\n", stepsize);
    printf("#  time   y           y'           yerror\n");

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