C код для простого маятника не дает ожидаемых результатов

Question

Я пытаюсь написать код для сети нейронов FitzHugh NAgumo в C, используя интегратор RK4. Поскольку это не сработало, я решил попробовать что-то более простое, простая маятниковая система. Я использую отдельные функции - один для возврата дифференциального массива, один для интегратора и несколько функций для добавления и умножения чисел на каждый элемент массива.

Вот код, который возвращает значения всех «0»:

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

Вот код, который делает вернуть синусоиду для омеги переменной на моем ноутбуке, но не на моем компьютере (и 64 бит, ноутбук под управлением Ubuntu 14.04 и ПК под управлением 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; 
} 

Я подозреваю, что проблема заключается главным образом в том, как я передаю массивы через функцию - как это работает, когда у меня нет отдельной функции для запуска интегратора. Поблагодарите за помощь.

Answer 1

Ваш код имеет неопределенное поведение. По своей сути это связано с тем, что вы используете статические указатели в своих функциях умножения и добавления. Позвольте мне сконденсировать его в проблемную часть:

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. 

Сначала избавьтесь от этого статического электричества. Затем вы можете сделать одну из двух вещей:

1. Просто выделите память внутри функций и верните указатель на нее. Это будет вызывающими коды ответственности освободить его в конце, как это:

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

2. передать указатель в ваших функции для их заполнений, оставляя управление памятью полностью вызывающего код:

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