1

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.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
Samyukta Ramnath
  • 371
  • 4
  • 18
  • You show two programs but only ask a question for the second. So what is the first program for? And what results do you get on the second program? – Werner Henze Feb 22 '17 at 08:40
  • give some idea about the result you are getting and what is the expected result. – Nish Feb 22 '17 at 08:48
  • I am actually asking about the first program. I wanted to highlight that the program with multiple modules does not give the expected result, while the code with one main function does give the expected result. The expected result being - values of theta and omega for different times, that form a sine wave when plotted against time. – Samyukta Ramnath Feb 22 '17 at 09:09

1 Answers1

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;
    }
}
StoryTeller - Unslander Monica
  • 165,132
  • 21
  • 377
  • 458
  • If I modify the code in the second way, then how do I, for example, return the differential array? For example, the function oscnetwork_opt sets an array with the coupled differential equations, which I return. If I write the following code: `void oscnetwork_opt(float y[2],float* dydt){ dydt[0] = y[1]; dydt[1] = -sin(y[0]); }` That wouldn't work either. – Samyukta Ramnath Feb 22 '17 at 09:30
  • @SamyuktaRamnath - In which way? I presented two options. – StoryTeller - Unslander Monica Feb 22 '17 at 09:32
  • Would I receive the results of my multiplication into a float* ptr in the calling function, if I tried the second method? In fact, trying this gives me a segmentation fault. – Samyukta Ramnath Feb 22 '17 at 09:37
  • @SamyuktaRamnath - You would do something like `float arr[2]; array_mul(2, array_in, num, arr);` The whole point is that you tell the function where to write, and decopule managing the memory from the calculation itself. – StoryTeller - Unslander Monica Feb 22 '17 at 09:43
  • Isn't that about the same as passing an array? – Samyukta Ramnath Feb 22 '17 at 11:02
  • 1
    @SamyuktaRamnath - No. You cannot pass arrays, since passing is always by value. This is an intricate point that novice programmers often fall over. A function can accept a pointer (to what is possibly the first element *of* an array), and arrays **decay** to pointers implicitly when you pass them. But you do not pass an array. `void foo(int a[2])` is just syntactic sugar for `void foo(int *a)`. There is no compiler check to make sure you pass an array of size 2 exactly. – StoryTeller - Unslander Monica Feb 22 '17 at 11:07