2017-02-22 2 views
-5

rk4 방법을 사용한 수치 적분과 간단한 진자에 대한 코드를 작성했습니다. Here's an image of expected result. 내 노트북에서 우분투 14.04, 64 비트 (결과적으로 사인파가 발생 함)가 작동하지만 데비안 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.