2017-02-22 122 views
-5

我已經使用rk4方法編寫了一個簡單擺錘的數值積分代碼。 Here's an image of expected result. 它運行在我的筆記本電腦上,運行Ubuntu 14.04,64位(它給出了一個正弦波),但不能在運行Debian 8的我的PC上運行,也是64位。 Here's an image of the wrong plot. 爲什麼會發生這種情況的任何原因?C數值積分代碼在一臺計算機上工作,但在另一臺計算機上爆炸

下面的代碼:

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

我們需要查看包含在int main()中的*最小*位代碼,以顯示此問題以幫助您。 – Bathsheba

+2

這看起來像未定義的行爲。也許有些變量沒有初始化,你可以在一臺計算機上而不是另一臺上運行。如果沒有看到你的實際節目,我們不能再多說了。 –

+1

你應該發佈某種代碼片段,因爲你的問題太模糊。 – Sitram

回答

3

你與你的變量youtintegrator_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

相關問題