RungeKutta method in C, pointers and arrays

Question

Now the program works, but is it correct? Sometimes it works in C anyway. Could I improve it somehow (not numerically, just C-wise)? I need to allocate the array as I do right? 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;
}

• nwellnhof: You shouldn't use a floating point comparison in the for loop condition. Because of rounding errors, it's essentially random how often the loop body is executed. Instead, pass resultSize to runge_kutta_2nd_order, compute stepSize from resultSize, and use index < resultSize as loop condition.

Answer 1

Generally, it is preferrable to put the \ns at the end of the corresponding printf instead of relying on them being present on the next one.

printf("\n"); // <-- you don't need to complicate your whole logic just because of this guy!
 // put him on a line of his own.
printf("\y'(t) - y(t) = e^(t/2) * (sin(5t) - 10cos(5t))\n");
printf("\y(0) = 0\n");
printf("0 <= t <= 10\n");
for(i=0; i<1001; i+=100){
 printf("i = %lf => y = %lf\n", 0.01*i, res[i]);
}

Also, as a general principle I don't like returning malloc-ed values from functions since doing so

1. forces you to remember to free the data even though the corresponding malloc is hidden away inside the other function and

2. couples your algorithm (the numeric calculations) with the memory management. It would be preferable if you were free to use any sort of memory instead of being forced into using malloc-free (as you currently are)

   One possibility is passing the result array as an argument to your function:

   //so you can malloc...
   double *result = malloc(/**/); 
   runge_kutta( /*...*/, result);
   free(result);
   //...or not
   double results[MAX_SIZE];
   runge_kutta( /*...*/, result);
   

   Do note that in the previous example you now need to be able to know the length of the result array beforehand. Most of the times this is trivial but in some cases you might need to package this in a subroutine.

Finally, since this is a numerically intensive problem and the list of results is not really needed (since I guess you only need to look at them one at a time and in sequential order) you might want to consider reorganizing things so that you can get away without any allocation for the results array.

For, example, we could have the Runge Kutta function call a callback on each result it gets instead of storing it to an array.

void runge_kutta_2nd_order(
 /*the original arguments and...*/,
 void (*onPoint)(int index, double x_i) /*the callback*/
)
{
 //...
 onPoint(index, ystar2);
 //...
}