qq_27104447 2015-04-22 12:30 采纳率: 40%
浏览 5469

C运行后老是结果等于-1.#IND00是怎么回事

我的题目是一个流体力学sod激波管问题,运行后结果都是-1.#IND00,还请各位大神指教啊
我程序源代码是:
#include
#include
float dx=1.0/100,e=0.01,r=1.4;
float a_1=-2/(30*dx),a_2=-5*a_1+1/(12*dx),a_3=10*a_1-2/(3*dx),a_4=-10*a_1,a_5=5*a_1+2/(3*dx),a_6=-a_1-1/(12*dx);
float*q[3];float*s[3];float*v[3];float*w1[3];float*w2[3];
float pr[100],u[100],p[100],c[100];
int i,j;
int main()
{
float dt;
float f[3][100],f_1[3][100],f_2[3][100],L[3][100],L_1[3][100],L_2[3][100];
float y[3][100],y_1[3][100],y_2[3][100],U[3][100],U_1[3][100],U_2[3][100],U_11[3][100],U_12[3][100],U_21[3][100],U_22[3][100];
float*(TIME(float f_1[3][100]));
float*(time(float f_2[3][100]));
float*(fun(float A[3][100]));
float*(FUN1(float B[3][100]));
float*(FUN2(float C[3][100]));
int n;
char*d="E:\x2.txt";char*a="E:\sudu2.txt";char*b="E:\midu2.txt";char*g="E:\yaqiang2.txt";
FILE*fd;FILE*fa;FILE*fb;FILE*fg;
fd=fopen(d,"a");fa=fopen(a,"a");fb=fopen(b,"a");fg=fopen(g,"a");
dt=0.01;
for(j=0;j<=99;j++)
{
if((j*dx)<0.5)
{
u[j]=0;pr[j]=1;p[j]=1;
}
else
{
u[j]=0;pr[j]=0.125;p[j]=0.1;
}
c[j]=sqrt(r*p[j]/pr[j]);
y[0][j]=u[j];y[1][j]=u[j]-c[j];y[2][j]=u[j]+c[j];
U[0][j]=pr[j];U[1][j]=u[j]*pr[j];U[2][j]=p[j]/(r-1)+pr[j]*pow(u[j],2)/2;
}
w1=FUN1(&y[0]);*w2=FUN2(&y[0]);
for(i=0;i<=2;i++)
{
for(j=0;j<=99;j++)
{
f_1[i][j]=
(*(w1+i)+j);f_2[i][j]=*(*(w2+i)+j);
}
}
for(n=1;n<=14;n++)
{
q=TIME(&f_1[0]);*s=time(&f_2[0]);
for(i=0;i<=2;i++)
{
for(j=0;j<=99;j++)
{
L_1[i][j]=
(*(q+i)+j);L_2[i][j]=*(*(s+i)+j);
L[i][j]=L_1[i][j]+L_2[i][j];
U_1[i][j]=U[i][j]+dt*L[i][j];
}
}
v=fun(&U_1[0]);
for(i=0;i<=2;i++)
{
for(j=0;j<=99;j++)
{
y[i][j]=
(*(v+i)+j);
}
}
w1=FUN1(&y[0]);*w2=FUN2(&y[0]);
for(i=0;i<=2;i++)
{
for(j=0;j<=99;j++)
{
U_11[i][j]=
(*(w1+i)+j);U_12[i][j]=*(*(w2+i)+j);
}
}

q=TIME(&U_11[0]);*s=time(&U_12[0]);
for(i=0;i<=2;i++)
{
for(j=0;j<=99;j++)
{
L_1[i][j]=
(*(q+i)+j);L_2[i][j]=*(*(s+i)+j);
L[i][j]=L_1[i][j]+L_2[i][j];
U_2[i][j]=U[i][j]*3/4+(U_1[i][j]+dt*L[i][j])*1/4;
}
}
v=fun(&U_2[0]);
for(i=0;i<=2;i++)
{
for(j=0;j<=99;j++)
{
y[i][j]=
(*(v+i)+j);
}
}
w1=FUN1(&y[0]);*w2=FUN2(&y[0]);
for(i=0;i<=2;i++)
{
for(j=0;j<=99;j++)
{
U_21[i][j]=
(*(w1+i)+j);U_22[i][j]=*(*(w2+i)+j);
}
}

q=TIME(&U_21[0]);*s=time(&U_22[0]);
for(i=0;i<=2;i++)
{
for(j=0;j<=99;j++)
{
L_1[i][j]=
(*(q+i)+j);L_2[i][j]=*(*(s+i)+j);
L[i][j]=L_1[i][j]+L_2[i][j];
U[i][j]=U[i][j]*1/3+(U_2[i][j]+dt*L[i][j])*2/3;
}
}
}
for(j=0;j<=99;j++)
{
pr[j]=U[0][j];u[j]=U[1][j]/pr[j];p[j]=(U[2][j]-U[0][j]*U[1][j]/2)*(r-1);
printf("%f,%f,%f,%f\n",j*dx,u[j],pr[j],p[j]);
fprintf(fd,"%f\n",j*dx);
fprintf(fa,"%f\n",u[j]);
fprintf(fb,"%f\n",pr[j]);
fprintf(fg,"%f\n",p[j]);
}

return 0;

}
float*(TIME(float f_1[3][100]))
{
float L_1[3][100];
for(i=0;i<=2;i++)
{
L_1[i][0]=-(a_1*f_1[i][0-3+100]+a_2*f_1[i][0-2+100]+a_3*f_1[i][0-1+100]+a_4*f_1[i][0]+a_5*f_1[i][0+1]+a_6*f_1[i][0+2]);
L_1[i][1]=-(a_1*f_1[i][1-3+100]+a_2*f_1[i][1-2+100]+a_3*f_1[i][1-1]+a_4*f_1[i][1]+a_5*f_1[i][1+1]+a_6*f_1[i][1+2]);
L_1[i][2]=-(a_1*f_1[i][2-3+100]+a_2*f_1[i][2-2]+a_3*f_1[i][2-1]+a_4*f_1[i][2]+a_5*f_1[i][2+1]+a_6*f_1[i][2+2]);
for(j=3;j<=97;j++)
{
L_1[i][j]=-(a_1*f_1[i][j-3]+a_2*f_1[i][j-2]+a_3*f_1[i][j-1]+a_4*f_1[i][j]+a_5*f_1[i][j+1]+a_6*f_1[i][j+2]);
}
L_1[i][98]=-(a_1*f_1[i][98-3]+a_2*f_1[i][98-2]+a_3*f_1[i][98-1]+a_4*f_1[i][98]+a_5*f_1[i][98+1]+a_6*f_1[i][98+2-100]);
L_1[i][99]=-(a_1*f_1[i][99-3]+a_2*f_1[i][99-2]+a_3*f_1[i][99-1]+a_4*f_1[i][99]+a_5*f_1[i][99+1-100]+a_6*f_1[i][99+2-100]);
q[i]=L_1[i];
}
return q;
}
float
(time(float f_2[3][100]))
{
float L_2[3][100];
for(i=0;i<=2;i++)
{
L_2[i][0]=a_1*f_2[i][0+3]+a_2*f_2[i][0+2]+a_3*f_2[i][0+1]+a_4*f_2[i][0]+a_5*f_2[i][0-1+100]+a_6*f_2[i][0-2+100];
L_2[i][1]=a_1*f_2[i][1+3]+a_2*f_2[i][1+2]+a_3*f_2[i][1+1]+a_4*f_2[i][1]+a_5*f_2[i][1-1]+a_6*f_2[i][1-2+100];
for(j=2;j<=96;j++)
{
L_2[i][j]=a_1*f_2[i][j+3]+a_2*f_2[i][j+2]+a_3*f_2[i][j+1]+a_4*f_2[i][j]+a_5*f_2[i][j-1]+a_6*f_2[i][j-2];
}
L_2[i][97]=a_1*f_2[i][97+3-100]+a_2*f_2[i][97+2]+a_3*f_2[i][97+1]+a_4*f_2[i][97]+a_5*f_2[i][97-1]+a_6*f_2[i][97-2];
L_2[i][98]=a_1*f_2[i][98+3-100]+a_2*f_2[i][98+2-100]+a_3*f_2[i][98+1]+a_4*f_2[i][98]+a_5*f_2[i][98-1]+a_6*f_2[i][98-2];
L_2[i][99]=a_1*f_2[i][99+3-100]+a_2*f_2[i][99+2-100]+a_3*f_2[i][99+1-100]+a_4*f_2[i][99]+a_5*f_2[i][99-1]+a_6*f_2[i][99-2];
s[i]=L_2[i];
}
return s;
}
float
(fun(float A[3][100]))
{
float y[3][100];
for(j=0;j<=99;j++)
{
pr[j]=A[0][j];
u[j]=A[1][j]/pr[j];
p[j]=(A[2][j]-A[0][j]*A[1][j]/2)*(r-1);
c[j]=sqrt(r*p[j]/pr[j]);
y[0][j]=u[j];y[1][j]=u[j]-c[j];y[2][j]=u[j]+c[j];
}
for(i=0;i<=2;i++)
{
v[i]=y[i];
}
return v;
}
float
(FUN1(float B[3][100]))
{
float y_1[3][100],f_1[3][100];
for(i=0;i<=2;i++)
{
for(j=0;j<=99;j++)
{
y_1[i][j]=(B[i][j]+sqrt(pow(B[i][j],2)+pow(e,2)))/2;
}

}
for(j=0;j<=99;j++)
{
f_1[0][j]=(2*(r-1)*y_1[0][j]+y_1[1][j]+y_1[2][j])*pr[j]/(2*r);
f_1[1][j]=(2*(r-1)*u[j]*y_1[0][j]+(u[j]-c[j])*y_1[1][j]+(u[j]+c[j])*y_1[2][j])*pr[j]/(2*r);
f_1[2][j]=((r-1)*y_1[0][j]*pow(u[j],2)+y_1[1][j]*pow(u[j]-c[j],2)/2+y_1[2][j]*pow(u[j]+c[j],2)/2+(3-r)*(y_1[1][j]+y_1[2][j])*pow(c[j],2)/(2*(r-1)))*pr[j]/(2*r);
}
for(i=0;i<=2;i++)
{
w1[i]=f_1[i];
}
return w1;
}
float
(FUN2(float C[3][100]))
{
float y_2[3][100],f_2[3][100];
for(i=0;i<=2;i++)
{
for(j=0;j<=99;j++)
{
y_2[i][j]=(C[i][j]-sqrt(pow(C[i][j],2)+pow(e,2)))/2;
}

}
for(j=0;j<=99;j++)
{
f_2[0][j]=(2*(r-1)*y_2[0][j]+y_2[1][j]+y_2[2][j])*pr[j]/(2*r);
f_2[1][j]=(2*(r-1)*u[j]*y_2[0][j]+(u[j]-c[j])*y_2[1][j]+(u[j]+c[j])*y_2[2][j])*pr[j]/(2*r);
f_2[2][j]=((r-1)*y_2[0][j]*pow(u[j],2)+y_2[1][j]*pow(u[j]-c[j],2)/2+y_2[2][j]*pow(u[j]+c[j],2)/2+(3-r)*(y_2[1][j]+y_2[2][j])*pow(c[j],2)/(2*(r-1)))*pr[j]/(2*r);
}
for(i=0;i<=2;i++)
{
w2[i]=f_2[i];
}
return *w2;
}

  • 写回答

3条回答

  • threenewbee 2015-04-23 12:28
    关注

    不关心你的算法,这个说明你的结果无限小。

    评论

报告相同问题?

悬赏问题

  • ¥15 全志H618ROM新增分区
  • ¥20 jupyter保存图像功能的实现
  • ¥15 在grasshopper里DrawViewportWires更改预览后,禁用电池仍然显示
  • ¥15 NAO机器人的录音程序保存问题
  • ¥15 C#读写EXCEL文件,不同编译
  • ¥15 MapReduce结果输出到HBase,一直连接不上MySQL
  • ¥15 扩散模型sd.webui使用时报错“Nonetype”
  • ¥15 stm32流水灯+呼吸灯+外部中断按键
  • ¥15 将二维数组,按照假设的规定,如0/1/0 == "4",把对应列位置写成一个字符并打印输出该字符
  • ¥15 NX MCD仿真与博途通讯不了啥情况