C Code for numerical integration works on one computer but blows up on another

Question

I have written a code for a simple pendulum with numerical integration using rk4 method. Here's an image of expected result. It works on my laptop, running Ubuntu 14.04, 64 bit, (it gives a sine wave as the result), but doesn't work on my PC, which runs Debian 8 and is also 64 bit. Here's an image of the wrong plot. Any reason why this would be happening?

Here's the code:

#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;
};


struct t_y_couple integrator_rk4(float dt, float t, float *p1);

void oscnetwork_opt(float t, float *y, float *dydt);

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] = 0;
  Y[0][1] = 3.14;
  //  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++)

    //      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));

    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[2], y2[2], y3[2], yout[2];
  float tout, dt_half;
  float k1[2], k2[2], k3[2], k4[2];
  //  initialize iterator
  int i;
  struct t_y_couple ty1;
  tout = t + dt;
  dt_half = 0.5 * dt;
  float addition[2];

  //  return the differential array into k1
  oscnetwork_opt(t, y, k1);
  //  multiply the array k1 by dt_half
  for (i = 0; i < 2; i++)
    y1[i] = y[i] + (k1[i]) * dt_half;
  //  add k1 to each element of the array y

  //  do the same thing 3 times
  oscnetwork_opt(t + dt_half, y1, k2);
  for (i = 0; i < 2; i++)
    y2[i] = y[i] + (k2[i]) * dt_half;

  oscnetwork_opt(t + dt_half, y2, k3);
  for (i = 0; i < 2; i++)
    y3[i] = y[i] + (k3[i]) * dt_half;
  oscnetwork_opt(tout, y3, k4);

  //  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;

  }
  //  add this to the original array
  for (i = 0; i < 2; i++)
    yout[i] = y[i] + addition[i];
  //  return a struct with the current time and the updated voltage array
  ty1.t = tout;
  ty1.y = yout;

  return ty1;
}

// function to return the vector with coupled differential variables for each
// time iteration
void oscnetwork_opt(float t, float y[2], float *dydt) {
  int i;
  dydt[0] = y[1];
  dydt[1] = -(1) * sin(y[0]);
}

We need to see the *smallest* bit of code wrapped in an `int main()` that exhibits this problem in order to help you. — Bathsheba, Feb 22 '17 at 11:08
This looks like undefined behavior. Maybe some variable which is not initialized and you get away with it on one computer and not on the other one. Whithout seeing your actual program we can't say much more. — Jabberwocky, Feb 22 '17 at 11:08
You should post some kind of code snippet because your question is too vague. — Sitram, Feb 22 '17 at 11:09
You don't seem to allocate any memory for the `y` member anywhere. — Lundin, Feb 22 '17 at 12:09
@Lundin , y is an input from the main function - I've passed a pointer to the array Y, by pointing p to the i'th row of Y in `p = &Y[i][0];` — Samyukta Ramnath, Feb 22 '17 at 12:18
Overall this program is very hard to follow. Try to keep things simple and use intuitive variable names. I suspect you have a large amount of subtle bugs lurking in this code. — Lundin, Feb 22 '17 at 12:42
`array_mul()` and `array_add()` are not used you should remove them and remove useless comment in your code like `// printf("final result : addition %f ", addition[i]);`. — Stargateur, Feb 22 '17 at 12:43

score 3 · Accepted Answer · edited May 23 '17 at 11:46

You have a problem of lifetime with your variable yout in integrator_rk4(). You assign address of yout to ty1.y but you use it outside this function. This is undefined behavior.

quick fix:

struct t_y_couple {
  float t;
  float y[2];
};

struct t_y_couple integrator_rk4(float dt, float t, float y[2]) {
  float y1[2], y2[2], y3[2], yout[2];

  // ...

  ty1.t = tout;
  ty1.y[0] = yout[0];
  ty1.y[1] = yout[1];

  return ty1;
}

You have a lot of useless allocation and you made "spaghetti code" with your global variable. You should not cast the return of malloc.

C Code for numerical integration works on one computer but blows up on another

1 Answers1