源代码为:
#include
#include
float C[20];
int main()
{
float dx,dt,t,u[20],u_1[20],u_2[20];
int n,;float L[20],*p;
float*TIME(float A[]);
dt=0.01;dx=2*M_PI/20;
for(n=0;n<=19;n++)
{
u[n]=sin(n*dx);
}
for(t=dt;;t=t+dt)
{
p=TIME(&u[0]);
for(n=0;n<=19;n++)
{
L[n]=*(p+n);
u_1[n]=u[n]+dt*L[n];
}
p=TIME(&u_1[0]);
for(n=0;n<=19;n++)
{
L[n]=*(p+n);
u_2[n]=u[n]*3/4+(u_1[n]+dt*L[n])*1/4;
}
p=TIME(&u_2[0]);
for(n=0;n<=19;n++)
{
L[n]=*(p+n);
u[n]=u[n]*1/3+(u_2[n]+dt*L[n])*2/3;printf("%f,%f,%f\n",t,n*dx,u_2[n]);
}
if (t==20||t==50)
{ for(n=0;n<20;n++)
printf("%f,%f,%f\n",t,n*dx,u[n]);
}
}
return 0;
}
float*TIME(float A[])
{
float a_1,a_2, a_3,a_4,a_5,a_6,dx;
int n;
dx=2*M_PI/20;
a_1=-6/(30*dx);
a_2=-5*a_1+1/(12*dx);
a_3=10*a_1-2/(3*dx);
a_4=-10*a_1+2/(3*dx);
a_5=5*a_1+2/(3*dx);
a_6=-a_1-1/(12*dx);
C[0]=a_1*A[0-3+20]+a_2*A[0-2+20]+a_3*A[0-1+20]+a_4*A[0]+a_5*A[0+1]+a_6*A[0+2];
C[1]=a_1*A[1-3+20]+a_2*A[1-2+20]+a_3*A[1-1]+a_4*A[1]+a_5*A[1+1]+a_6*A[1+2];
C[2]=a_1*A[2-3+20]+a_2*A[2-2]+a_3*A[2-1]+a_4*A[2]+a_5*A[2+1]+a_6*A[2+2];
for(n=3;n<=17;n++)
{
C[n]=a_1*A[n-3]+a_2*A[n-2]+a_3*A[n-1]+a_4*A[n]+a_5*A[n+1]+a_6*A[n+2];
}
C[18]=a_1*A[18-3]+a_2*A[18-2]+a_3*A[18-1]+a_4*A[18]+a_5*A[18+1]+a_6*A[18+2-20];
C[19]=a_1*A[19-3]+a_2*A[19-2]+a_3*A[19-1]+a_4*A[19]+a_5*A[19+1-20]+a_6*A[19+2-20];
return &C[0];
}
这是浮点类型,浮点类型不精确,0.01不一定刚好是20,可以是19.9999之类。建议把0.01放大100倍后在循环使用
if (t==20||t==50)
你把步进数据都有int类型,然后都扩大100倍,这样步进就是1,然后在循环中才处理转换double
#include
#include
float C[20];
int main()
{
float dx,dt,t,u[20],u_1[20],u_2[20];
int n;float L[20],*p;
float*TIME(float A[]);
dt=0.01;dx=2*M_PI/20;
for(n=0;n<=19;n++)
{
u[n]=sin(n*dx);
}
for(t=dt;;t=t+dt)
{
p=TIME(&u[0]);
for(n=0;n<=19;n++)
{
L[n]=*(p+n);
u_1[n]=u[n]+dt*L[n];
}
p=TIME(&u_1[0]);
for(n=0;n<=19;n++)
{
L[n]=*(p+n);
u_2[n]=u[n]*3/4+(u_1[n]+dt*L[n])*1/4;
}
p=TIME(&u_2[0]);
for(n=0;n<=19;n++)
{
L[n]=*(p+n);
u[n]=u[n]*1/3+(u_2[n]+dt*L[n])*2/3;printf("%f,%f,%f\n",t,n*dx,u_2[n]);
}
if (t==20||t==50)
{ for(n=0;n<20;n++)
printf("%f,%f,%f\n",t,n*dx,u[n]);
}
}
return 0;
}
float*TIME(float A[])
{
float a_1,a_2, a_3,a_4,a_5,a_6,dx;
int n;
dx=2*M_PI/20;
a_1=-6/(30*dx);
a_2=-5*a_1+1/(12*dx);
a_3=10*a_1-2/(3*dx);
a_4=-10*a_1+2/(3*dx);
a_5=5*a_1+2/(3*dx);
a_6=-a_1-1/(12*dx);
C[0]=a_1*A[0-3+20]+a_2*A[0-2+20]+a_3*A[0-1+20]+a_4*A[0]+a_5*A[0+1]+a_6*A[0+2];
C[1]=a_1*A[1-3+20]+a_2*A[1-2+20]+a_3*A[1-1]+a_4*A[1]+a_5*A[1+1]+a_6*A[1+2];
C[2]=a_1*A[2-3+20]+a_2*A[2-2]+a_3*A[2-1]+a_4*A[2]+a_5*A[2+1]+a_6*A[2+2];
for(n=3;n<=17;n++)
{
C[n]=a_1*A[n-3]+a_2*A[n-2]+a_3*A[n-1]+a_4*A[n]+a_5*A[n+1]+a_6*A[n+2];
}
C[18]=a_1*A[18-3]+a_2*A[18-2]+a_3*A[18-1]+a_4*A[18]+a_5*A[18+1]+a_6*A[18+2-20];
C[19]=a_1*A[19-3]+a_2*A[19-2]+a_3*A[19-1]+a_4*A[19]+a_5*A[19+1-20]+a_6*A[19+2-20];
return &C[0];
}
我的源代码是这样的,题目要求是:已知上一时刻的u,然后按以下的递推式子:u1=u+dt*L(u); u2=3/4*u+1/4*(u1+dt*L(u1)); u=1/3*u+2/3*(u2+dt*L(u2))求下一时刻的u,最终得出20s和50s时的u,L(u)的式子是已知的,就是调用函数里的式子,可不可以帮我看一下程序有没有错误,运行出来的结果很奇怪、、、