-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]);
}
我們需要查看包含在int main()中的*最小*位代碼,以顯示此問題以幫助您。 – Bathsheba
這看起來像未定義的行爲。也許有些變量沒有初始化,你可以在一臺計算機上而不是另一臺上運行。如果沒有看到你的實際節目,我們不能再多說了。 –
你應該發佈某種代碼片段,因爲你的問題太模糊。 – Sitram