问题描述
我在用 C 编写的用于求解一维线性对流方程的程序有问题。基本上我已经初始化了两个数组。第一个数组 (u0_array) 是一个数组,数组元素在区间内设置为等于 2 或 0.5
我遇到的问题是代码末尾的嵌套 for 循环。此循环将计算下一个点所需的更新方程应用于数组中的每个元素。当我运行我的脚本并尝试打印结果时,对于循环的每次迭代,我得到的输出只是 -1.IND00。 (我正在遵循我也在下面附上的 Matlab 风格的伪代码)我对 C 很陌生,所以我的经验不足。我不知道为什么会这样。如果有人可以提出可能的解决方案,我将不胜感激。到目前为止,我已经在下面附上了我的代码和一些评论,以便您可以遵循我的思考过程。我还附上了我遵循的 Matlab 风格的伪代码。
# include <math.h>
# include <stdlib.h>
# include <stdio.h>
# include <time.h>
int main ( )
{
//eqyation to be solved 1D linear convection -- du/dt + c(du/dx) = 0
//initial conditions u(x,0) = u0(x)
//after discretisation using forward difference time and backward differemce space
//update equation becomes u[i] = u0[i] - c * dt/dx * (u0[i] - u[i - 1]);
int nx = 41; //num of grid points
int dx = 2 / (nx - 1); //magnitude of the spacing between grid points
int nt = 25;//nt is the number of timesteps
double dt = 0.25; //the amount of time each timestep covers
int c = 1; //assume wavespeed
//set up our initial conditions. The initial veLocity u_0
//is 2 across the interval of 0.5 <x < 1 and u_0 = 1 everywhere else.
//we will define an array of ones
double* u0_array = (double*)calloc(nx,sizeof(double));
for (int i = 0; i < nx; i++)
{
u0_array[i] = 1;
}
// set u = 2 between 0.5 and 1 as per initial conditions
//note 0.5/dx = 10,1/dx+1 = 21
for (int i = 10; i < 21; i++)
{
u0_array[i] = 2;
//printf("%f,",u0_array[i]);
}
//make a temporary array that allows us to store the solution
double* usol_array = (double*)calloc(nx,sizeof(double));
//apply numerical scheme forward difference in
//time an backward difference in space
for (int i = 0; i < nt; i++)
{
//first loop iterates through each time step
usol_array[i] = u0_array[i];
//printf("%f",usol_array[i]);
//MY CODE WORKS FINE AS I WANT UP TO THIS LOOP
//second array iterates through each grid point for that time step and applies
//the update equation
for (int j = 1; j < nx - 1; j++)
{
u0_array[j] = usol_array[j] - c * dt/dx * (usol_array[j] - usol_array[j - 1]);
printf("%f,u0_array[j]);
}
}
return EXIT_SUCCESS;
}
下面附上我所遵循的伪代码以供参考 1D linear convection pseudo Code (Matlab Style)
解决方法
使用 FP 数学代替整数除法。
这避免了后面的除以 0 和 -1.#IND00
// int dx = 2 / (nx - 1); quotient is 0.
double dx = 2.0 / (nx - 1);
OP 的代码与注释不匹配
// u[i] = u0[i] - c * dt/dx * (u0[i] - u[i - 1]
u0_array[j] = usol_array[j] - c * dt/dx * (usol_array[j] - usol_array[j - 1]);
更容易查看是否删除了多余的 _array
。
// v correct?
// u[i] = u0[i] - c * dt/dx * (u0[i] - u[i - 1]
u0[j] = usol[j] - c * dt/dx * (usol[j] - usol[j - 1]);
// ^^^^ Correct?
,
你可能想要的是 matlab 术语
for i = 1:nt
usol = u0
u0(2:nx) = usol(2:nx) - c*dt/dx*(usol(2:nx)-usol(1:nx-1))
end%for
这意味着您在空间维度上有 2 个内部循环,两个向量操作中的每一个都有一个。这两个单独的循环必须在 C 代码中显式
//apply numerical scheme forward difference in
//time an backward difference in space
for (int i = 0; i < nt; i++) {
//first loop iterates through each time step
for (int j = 0; j < nx; j++) {
usol_array[j] = u0_array[j];
//printf("%f",usol_array[i]);
}
//second loop iterates through each grid point for that time step
//and applies the update equation
for (int j = 1; j < nx; j++)
{
u0_array[j] = usol_array[j] - c * dt/dx * (usol_array[j] - usol_array[j - 1]);
printf("%f,",u0_array[j]);
}
}