# Runge-Kutta 2nd order of differential equations

Could I improve this somehow (not numerically, just C-wise)? Do I need to allocate the array since I am returning it and not just using it inside the function?

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

/* Approximates a solution to a differential equation on the form:
y'(t) + ay(t) = x(t)
y(0) = b
*/
double* runge_kutta_2nd_order(double stepSize, double a, double b, double (*x) (double), double upto)
{
int resultSize = ((int) (upto / stepSize)) + 2;
double yt = b;
double time;
double k1,k2,ystar1,ystar2;
int index = 1;

//double *results = (double*) malloc(resultSize * (sizeof(double)));
double *results = malloc(resultSize * sizeof(*results));
if(results == NULL)
{
printf("\nCould not allocate memory. Exiting program.");
exit(0);
}

results[0] = b;

for(time = 0; time < upto; time += stepSize) //<=
{
k1 = x(time) - a * yt;
ystar1 = yt + stepSize * k1;
k2 = x(time + stepSize) - a * ystar1;
ystar2 = yt + (k1 + k2) / 2 * stepSize;
yt = ystar2;
results[index] = ystar2;
index++;
}
return results;
}

void free_results(double **r)
{
free(*r);
*r = NULL;
}

double insignal(double t)
{
return exp(t/2)*(sin(5*t) - 10*cos(5*t));
}

int main(void)
{
int i;
double *res = runge_kutta_2nd_order(0.01,-1,0,&insignal,10);

printf("\nRunge Kutta 2nd order approximation of the differential equation:");
printf("\ny'(t) - y(t) = e^(t/2) * (sin(5t) - 10cos(5t))");
printf("\ny(0) = 0");
printf("\n0 <= t <= 10");

for(i=0; i<1001; i+=100){
printf("\ni = %lf => y = ", 0.01*i);
printf("%lf", res[i]);
}
printf("\n");

free_results(&res);

return 0;
}

• Is this a homework?
– jsalonen
Dec 4 '11 at 20:36
• I can guarantee you that your program correctly does what it does. Dec 4 '11 at 20:38
• no. the homework was implementing range-kutta in mathematica in mathematica and it is done. now im implementing it in C just to learn.
– Smith Johnson
Dec 4 '11 at 20:38
• I would consider splitting this one huge question into several, more specific subquestions. As the question stands, it may not be very helpful for other users.
– jsalonen
Dec 4 '11 at 20:40
• Im new with C. Just wondering if everything with pointers and arrays is correct. If I am allocating correctly and freeing correctly.
– Smith Johnson
Dec 4 '11 at 20:42

It is more usual to take a pointer to an already-allocated array, rather than allocate an array within the function itself. This allows the calling function to choose the best kind of allocation.

When doing this you'd likely also take the number of indexes to fill as an input, deriving the step size from that:

/* resultsize must be >= 2 */
void runge_kutta_2nd_order(double results[], size_t resultsize, double a, double b, double (*x) (double), double upto)
{
const double stepSize = upto / (resultsize - 1);
double yt = b;
size_t index;

results[0] = b;

for (index = 1; index < resultsize; index++)
{
const double time = index * stepSize;
const double k1 = x(time) - a * yt;
const double ystar1 = yt + stepSize * k1;
const double k2 = x(time + stepSize) - a * ystar1;
const double ystar2 = yt + (k1 + k2) / 2 * stepSize;
yt = ystar2;
results[index] = ystar2;
}
}


So, for example, in this case you could use automatic allocation in main():

int main(void)
{
int i;
double res[1001];

runge_kutta_2nd_order(res, 1001, -1, 0, &insignal, 10);

printf("\nRunge Kutta 2nd order approximation of the differential equation:");
printf("\ny'(t) - y(t) = e^(t/2) * (sin(5t) - 10cos(5t))");
printf("\ny(0) = 0");
printf("\n0 <= t <= 10");

for(i=0; i<1001; i+=100){
printf("\ni = %f => y = ", 0.01*i);
printf("%f", res[i]);
}
printf("\n");

return 0;
}