关于sod激波管的数值解问题

各位大神，我是计算流体力学的初学者，现在有一道关于sod激波管的数值解问题，我用C言编程写的，可就是结果不对，各位大神能不能给看一下、、题目要求是，初始条件x<0.5时，速度密度压强分别为0,、1、1，x>0.5时，速度密度压强分别为0,、0.125、1，求0.14时刻的速度密度和压强，空间差分采用5阶迎风，时间用R_K3阶格式，采用Sterger_Warming分裂，源程序如下：
#include<stdio.h>
#include<math.h>
float dx=1.0/100,e=0.1,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 a_1=0,a_2=0,a_3=-1,a_4=1,a_5=0,a_6=0;
float*q[3];float*s[3];float*v[3];float*w1[3];float*w2[3];
float pr[100],u[100],p[100],c[100],A[3][100],B[3][100],E[3][100],F[3][100],G[3][100];
float f_1[3][100],f_2[3][100],L_1[3][100],L_2[3][100],y[3][100],y_1[3][100],y_2[3][100];
int i,j;
int main()
{
float dt;
//float 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];
float L[3][100],U[3][100],U_1[3][100],U_2[3][100];
float*(TIME(float E[3][100]));
float*(time(float F[3][100]));
float*(fun(float A[3][100]));
float*(FUN1(float B[3][100]));
float*(FUN2(float G[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.0001;
for(j=0;j<=99;j++)
{
if(j<50)
{
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<=1400;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++)
{
f_1[i][j]=*(*(w1+i)+j);f_2[i][j]=*(*(w2+i)+j);
}
}
*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_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++)
{
f_1[i][j]=*(*(w1+i)+j);f_2[i][j]=*(*(w2+i)+j);
}
}
*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[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]-pow(u[j],2)*pr[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 E[3][100]))
{
float L_1[3][100];
for(i=0;i<=2;i++)
{
L_1[i][0]=-(a_1*E[i][0]+a_2*E[i][0]+a_3*E[i][0]+a_4*E[i][0]+a_5*E[i][0+1]+a_6*E[i][0+2]);
L_1[i][1]=-(a_1*E[i][0]+a_2*E[i][0]+a_3*E[i][1-1]+a_4*E[i][1]+a_5*E[i][1+1]+a_6*E[i][1+2]);
L_1[i][2]=-(a_1*E[i][0]+a_2*E[i][2-2]+a_3*E[i][2-1]+a_4*E[i][2]+a_5*E[i][2+1]+a_6*E[i][2+2]);
for(j=3;j<=97;j++)
{
L_1[i][j]=-(a_1*E[i][j-3]+a_2*E[i][j-2]+a_3*E[i][j-1]+a_4*E[i][j]+a_5*E[i][j+1]+a_6*E[i][j+2]);
}
L_1[i][98]=-(a_1*E[i][98-3]+a_2*E[i][98-2]+a_3*E[i][98-1]+a_4*E[i][98]+a_5*E[i][98+1]+a_6*E[i][99]);
L_1[i][99]=-(a_1*E[i][99-3]+a_2*E[i][99-2]+a_3*E[i][99-1]+a_4*E[i][99]+a_5*E[i][99]+a_6*E[i][99]);
q[i]=L_1[i];
}
return *q;
}
float*(time(float F[3][100]))
{
float L_2[3][100];
for(i=0;i<=2;i++)
{
L_2[i][0]=a_1*F[i][0+3]+a_2*F[i][0+2]+a_3*F[i][0+1]+a_4*F[i][0]+a_5*F[i][0]+a_6*F[i][0];
L_2[i][1]=a_1*F[i][1+3]+a_2*F[i][1+2]+a_3*F[i][1+1]+a_4*F[i][1]+a_5*F[i][1-1]+a_6*F[i][0];
for(j=2;j<=96;j++)
{
L_2[i][j]=a_1*F[i][j+3]+a_2*F[i][j+2]+a_3*F[i][j+1]+a_4*F[i][j]+a_5*F[i][j-1]+a_6*F[i][j-2];
}
L_2[i][97]=a_1*F[i][99]+a_2*F[i][97+2]+a_3*F[i][97+1]+a_4*F[i][97]+a_5*F[i][97-1]+a_6*F[i][97-2];
L_2[i][98]=a_1*F[i][99]+a_2*F[i][99]+a_3*F[i][98+1]+a_4*F[i][98]+a_5*F[i][98-1]+a_6*F[i][98-2];
L_2[i][99]=a_1*F[i][99]+a_2*F[i][99]+a_3*F[i][99]+a_4*F[i][99]+a_5*F[i][99-1]+a_6*F[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]/A[0][j];
p[j]=(A[2][j]-pow(u[j],2)*pr[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 G[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]=(G[i][j]-sqrt(pow(G[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;
}