2017-02-22 17 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数値積分のコードは1つのコンピュータで動作しますが、別のコンピュータで動作します

#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

関連する問題