C Код для численного интегрирования работает на одном компьютере, но взрывается другим

Question

Я написал код для простого маятника с численным интегрированием с использованием метода rk4. Here's an image of expected result. Он работает на моем ноутбуке, работает Ubuntu 14.04, 64 бит (в результате он дает синусоидальную волну), но не работает на моем ПК, который запускает Debian 8, а также 64 бит. Here's an image of the wrong plot. Любая причина, почему это происходит?

Вот код:

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

Answer 1

У вас есть проблемы жизни с вашей переменной yout в integrator_rk4(). Вы назначаете адрес yout в ty1.y, но используете его вне этой функции. Это неопределенное поведение.

быстрое решение:

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 should not cast the return of malloc.