How do I write an RK4 (4th order Runge-Kutta) code in C which computes a system

Question

How do I write an RK4 (4th order Runge-Kutta) code in C which computes a system of three, first order, differential equations? I have started a code in C which computes a single first order differential equation I just need help converting that to solve this system:

S' = bSZ

Z' = bSZ + cR -aSZ

R' = dS + aSZ - cR

where a,b,c, and d are constants with values:

a = 0.005, b = 0.0095, c = 0.0001, d = 0.0001

The initial values of S, R, and Z are:

S(0) = 500

Z(0) = 0

R(0) = 0

The C code I have solves the differential equation:

Here is the C code:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

double rk4(double(*f)(double, double), double h, double t, double y)
{
   double   dy1 = h * f(t, y),
       dy2 = h * f(t + h / 2, y + dy1 / 2),
       dy3 = h * f(t + h / 2, y + dy2 / 2),
       dy4 = h * f(t + h, y + dy3);

   return y + (dy1 + 2 * dy2 + 2 * dy3 + dy4) / 6;
}

double rate(double t, double y)
{
   return t * sqrt(y);
}

int main(void)
{
   double *yout, t, y2;
   double t0 = 0, t1 = 10, h = .1;
   int i, n = 1 + (t1 - t0)/h;
   yout = malloc(sizeof(double) * n);

   for (yout[0] = 1, i = 1; i < n; i++)
       yout[i] = rk4(rate, h, t0 + h * (i - 1), yout[i-1]);

   printf("x yout rel. err. ------------ ");
   for (i = 0; i < n; i += 10) {
       t = t0 + h * i;
       y2 = pow(t * t / 4 + 1, 2);
       printf("%g %g %g ", t, yout[i], yout[i]/y2 - 1);
   }

   return 0;
}

Explanation / Answer

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

double rk4(double(*f)(double, double), double h, double t, double y)
{
   double   dy1 = h * f(t, y),
       dy2 = h * f(t + h / 2, y + dy1 / 2),
       dy3 = h * f(t + h / 2, y + dy2 / 2),
       dy4 = h * f(t + h, y + dy3);

   return y + (dy1 + 2 * dy2 + 2 * dy3 + dy4) / 6;
}

double rate(double t, double y)
{
   return t * sqrt(y);
}

int main(void)
{
   double *yout, t, y2;
   double t0 = 0, t1 = 10, h = .1;
   int i, n = 1 + (t1 - t0)/h;
   yout = malloc(sizeof(double) * n);

   for (yout[0] = 1, i = 1; i < n; i++)
       yout[i] = rk4(rate, h, t0 + h * (i - 1), yout[i-1]);

   printf("x yout rel. err. ------------ ");
   for (i = 0; i < n; i += 10) {
       t = t0 + h * i;
       y2 = pow(t * t / 4 + 1, 2);
       printf("%g %g %g ", t, yout[i], yout[i]/y2 - 1);
   }

   return 0;
}

Related Questions

Navigate

• Browse (All)
• Browse H
• Subjects

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.