C code for simple pendulum doesn't give expected results

Question

I am trying to write a code for a network of FitzHugh NAgumo neurons in C, using RK4 integrator. As that wasn't working, I decided to try something simpler, a simple pendulum system. I am using separate functions - one to return the differential array, one for the integrator, and a few functions to add and multiply numbers to each element of an array.

Here is the code that returns all '0' values:

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

int N = 2;

float h = 0.001;

struct t_y_couple{
    float t;
    float* y;
};


float* array_add(int len_array_in,float array_in[2], float array_add[2]);
struct t_y_couple integrator_rk4(float dt,float t, float* p1);

float* oscnetwork_opt(float t, float *y);
float* array_mul(int len_array_in,float array_in[2], float num);


int main(void){
    /* initializations*/
    struct t_y_couple t_y;


    int i,iter;

//  time span for which to run simulation
    int tspan = 20;


//  total number of time iterations = tspan*step_size   
    int tot_time = (int) ceil(tspan/h);

//  Time array
    float T[tot_time];

//  pointer definitions
    float *p, *q;


//  vector to hold values for each differential variable for all time iterations
    float Y[tot_time][2];
//  2*N+sum_adj = total number of coupled differential equations to solve   

//  initial conditions vector for time = 0

    Y[0][0] = 0;
    Y[0][1] = 0.1;
//  set the time array
    T[0] = 0;

//  This loop calls the RK4 code
    for (i=0;i<tot_time-1;i++){
        p = &Y[i][0];   // current time
        q = &Y[i+1][0]; // next time step
//      printf("\n\n");
//      for (j=0;j<2*N+sum_adj;j++)
//          printf("value passed in from main function: %f ,%d ",*(p+j),j);
//      printf("\n");

//      call the RK4 integrator with current time value, and current 
//      values of voltage           
        t_y = integrator_rk4(h,T[i],p);  

//      Return the time output of integrator into the next iteration of time
        T[i+1] = t_y.t; 

//      copy the output of the integrator into the next iteration of voltage        
        q = memcpy(q, t_y.y, (2) * sizeof(float));
//      q = memcpy(q, p, (2*N+sum_adj) * sizeof(float) );

        printf("%f ",T[i+1]);
        for (iter = 0;iter<N;iter++)
            printf("%f ",*(p+iter));
        printf("\n");
    }


    return 0;
}

struct t_y_couple integrator_rk4(float dt,float t, float y[2])
{   
//  initialize all the pointers
    float *y1,*y2,*y3, *yout;
    float tout,dt_half;
    float *k1,*k2,*k3,*k4;
//  initialize iterator
    int i;

//  printf("\n");
    struct t_y_couple ty1;
    tout = t+dt;
    dt_half = 0.5*dt;
    float addition[2];

//  return the differential array into k1
    k1 = oscnetwork_opt(t,y);
//  multiply the array k1 by dt_half    
    k1 = array_mul(2,k1,dt_half);
//  add k1 to each element of the array y passed in
    y1 = array_add(2,y,k1);

//  for (i=0;i<2*N+sum_adj;i++)
//      printf("k1: %f y: %f y1: %f\n",*(k1+i), *(y+i), *(y1+i));

//  do the same thing 3 times
    k2 = oscnetwork_opt(t+dt_half,y1);
    k2 = array_mul(2,k2,dt_half);
    y2 = array_add(2,y,k2); 
//  for (i=0;i<2*N+sum_adj;i++)
//      printf("k2: %f y: %f y2: %f\n",*(k2+i), *(y+i), *(y2+i));

    k3 = oscnetwork_opt(t+dt_half,y2);
    k3 = array_mul(2,k3,dt);
    y3 = array_add(2,y,k3);
//  for (i=0;i<2*N+sum_adj;i++)
//      printf("k3: %f y: %f y3: %f\n",*(k3+i), *(y+i), *(y3+i));

    k4 = oscnetwork_opt(tout,y3);
    k4 = array_mul(2,k4,dt);
    yout = array_add(2,y,k4);
//  for (i=0;i<2*N+sum_adj;i++)
//      printf("k4: %f y: %f y4: %f\n",*(k4+i), *(y+i), *(yout+i));     

//  Make the final additions with k1,k2,k3 and k4 according to the RK4 code
    for (i=0;i<2;i++){
        addition[i] = ((*(k1+i)) + (*(k2+i))*2 + (*(k3+i))*2 + (*(k4+i))) *dt/6;
//      printf("final result : addition %f  ", addition[i]);
    }
//  printf("\n");
//  add this to the original array
    yout = array_add(2,y,addition);

//  for (i=0;i<2*N+sum_adj;i++)
//      printf("yout: %f  ",*(yout+i));
//  printf("\n");
//  return a struct with the current time and the updated voltage array
    ty1.t = tout;
    ty1.y = yout;

    return ty1;
}



// function to add two arrays together, element by element
float* array_add(int len_array_in,float array_in[2], float array_sum[2]){
    int i;
    static float *array_out_add= NULL;
    if (array_out_add != 0) {
        array_out_add = (float*) realloc(array_out_add, sizeof(float) * (2));
    } else {
        array_out_add = (float*) malloc(sizeof(float) * (2));
    }

    for (i=0;i<len_array_in;i++){
        array_out_add[i] = array_in[i]+array_sum[i];
//      printf("before adding: %f, add amount: %f , after adding: %f, iteration: %d\n ", array_in[i], array_sum[i], array_out[i],i);
    }
    return array_out_add;
//  return 0;
}


// function to multiply each element of the array by some number
float* array_mul(int len_array_in,float array_in[2], float num){
    int i;
    static float *array_out_mul= NULL;
    if (array_out_mul != 0) {
        array_out_mul = (float*) realloc(array_out_mul, sizeof(float) * (2));
    } else {
        array_out_mul = (float*) malloc(sizeof(float) * (2));
    }
    for (i=0;i<len_array_in;i++){
        array_out_mul[i] =array_in[i]*num;
    }
    return array_out_mul;
//  return 0;
}


// function to return the vector with coupled differential variables for each time iteration
float* oscnetwork_opt(float t, float *y){

//  define and allocate memory for the differential vector
    static float* dydt = NULL;
    if (dydt != 0) {
        dydt = (float*) realloc(dydt, sizeof(float) * (2));
    } else {
        dydt = (float*) malloc(sizeof(float) * (2));
    }


    dydt[0] = y[1];
    dydt[1] = -(0.1*0.1)*sin(y[0]);
    return dydt;
}

Here is a code that does return a sine wave for the omega variable on my laptop but not on my PC (both 64 bit, laptop running Ubuntu 14.04 and PC running Debian 8).

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

int N = 2;

float h = 0.0001;

struct t_y_couple{
    float t;
    float* y;
};


float* array_add(int len_array_in,float array_in[2], float array_add[2]);
struct t_y_couple integrator_rk4(float dt,float t, float* p1);

float* oscnetwork_opt(float t, float *y);
float* array_mul(int len_array_in,float array_in[2], float num);


int main(void){
    /* initializations*/
    struct t_y_couple t_y;


    int i,iter,j;

//  time span for which to run simulation
    int tspan = 20;


//  total number of time iterations = tspan*step_size   
    int tot_time = (int) ceil(tspan/h);

//  Time array
    float T[tot_time];

//  pointer definitions
    float *p, *q;


//  vector to hold values for each differential variable for all time iterations
    float Y[tot_time][2];
//  N = total number of coupled differential equations to solve     

//  initial conditions vector for time = 0

    Y[0][0] = 3.14;
    Y[0][1] = 0;
//  set the time array
    T[0] = 0;

//  This loop calls the RK4 code
    for (i=0;i<tot_time-1;i++){
        p = &Y[i][0];   // current time
//      q = &Y[i+1][0]; // next time step
//      printf("\n\n");
//      for (j=0;j<N;j++)
//          printf("value passed in from main function: %f ,%d ",*(p+j),j);
//      printf("\n");

//      call the RK4 integrator with current time value, and current 
//      values of voltage           
        t_y = integrator_rk4(h,T[i],p);  

//      Return the time output of integrator into the next iteration of time
        T[i+1] = t_y.t; 

//      copy the output of the integrator into the next iteration of voltage        
//      q = memcpy(q, t_y.y, (2) * sizeof(float));
//      q = memcpy(q, p, (N) * sizeof(float) );

        printf("%f ",T[i+1]);
        for (iter = 0;iter<N;iter++){
            Y[i+1][iter] = t_y.y[iter];
            printf("%f ",Y[i+1][iter]);
        }
        printf("\n");
    }


    return 0;
}

struct t_y_couple integrator_rk4(float dt,float t, float y[2])
{   
//  initialize all the pointers
    float y1[2],y2[2],y3[2], yout[2];
    float tout,dt_half;
    float k1[2],k2[2],k3[2],k4[2];
//  initialize iterator
    int i;
//  for(i=0;i<N;i++)
//      printf("into integrator %f ",y[i]);
//  printf("\n");
    struct t_y_couple ty1;
    tout = t+dt;
    dt_half = 0.5*dt;
    float addition[2];

//  return the differential array into k1
    k1[0] = y[1];
    k1[1] = -(1)*sin(y[0]);


//  multiply the array k1 by dt_half
    for (i=0;i<N;i++){
        y1[i] = y[i] + k1[i]*dt_half;
//      printf("k1: %f y: %f y1: %f\n",*(k1+i), *(y+i), *(y1+i));       
    }

//  do the same thing 3 times
    k2[0] = y1[1];
    k2[1] = -(1)*sin(y1[0]);


    for (i=0;i<N;i++){
        y2[i] = y[i] + k2[i]*dt_half;
//      printf("k2: %f y: %f y2: %f\n",*(k2+i), *(y+i), *(y2+i));       
    }
    k3[0] = y2[1];
    k3[1] = -(1)*sin(y2[0]);


    for (i=0;i<N;i++){
        y3[i] = y[i] + k3[i]*dt;
//      printf("k3: %f y: %f y3: %f\n",*(k3+i), *(y+i), *(y3+i));
    }


    k4[0] = y3[1];
    k4[1] = -(1)*sin(y3[0]);


//  Make the final additions with k1,k2,k3 and k4 according to the RK4 code

    for (i=0;i<N;i++){

        addition[i] = (( k1[i]*dt) + (k2[i])*2 + (k3[i])*2 + (k4[i]))/6;
        yout[i] = y[i] + addition[i];
//      printf("y[%d]: %f  +  addition[%d]: %f  = yout[%d] :%f  ",i,y[i],i, addition[i], i, yout[i]);

//      printf("final result : addition %f  ", addition[i]);
    }
//  add this to the original array

//  printf("\n");
//  return a struct with the current time and the updated voltage array
    ty1.t = tout;
    ty1.y = yout;

    return ty1;
}

I suspect that the problem lies mostly in the way I am passing the arrays through the functions - as it works when I do not have a separate function to run the integrator. Would appreciate some help on this.

Answer 1

Your code has undefined behavior. At its core, it's because you use static pointers in your multiplication and addition functions. Allow me to condense it into the problematic part:

k2 = array_mul(2,k2,dt_half);
k3 = array_mul(2,k3,dt); // This either makes k2 point to freed memory, or overwrites the values in the location it points to.

addition[i] = ((*(k1+i)) + (*(k2+i))*2 + (*(k3+i))*2 + (*(k4+i))) *dt/6; // This uses all four pointers as though they must point to distinct memory locations.

First get rid of that static. Then you can do one of two things:

1. Simply allocate memory inside the functions and return a pointer to it. It will be the calling codes responsibility to free it at the end, like this:

   free(k1);
   free(k2);
   // etc
   

2. Pass a pointer into a your functions for them to populate, leaving the memory management entirely to the calling code:

   // function to multiply each element of the array by some number
   void array_mul(int len_array_in,float *array_in, float num, float *array_out){
       int i;
       for (i=0;i<len_array_in;i++){
         array_out[i] =array_in[i]*num;
       }
   }