# 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;
}
```