偏微分方程的数值模拟的程序，输出的解为0，求修改
#include <cstdlib>#include <cmath>
#include <string>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
using namespace std;
const int Q = 4;
const int NX =400;
int e[Q] = {1,-1,2,-2};
double s[Q] ={0.7,19.0/30.0,-11.0/60.0,-3.0/20.0};
double u[NX+3],u1[NX+3],u2[NX+3],Du[NX+3],u0[NX+3],ut[NX+3],utt[NX+3],ut0[NX+3],uxx[NX+3],f[NX+3][Q],F[NX+3][Q];
double g[NX+3],sour_xx[NX+3];
double c,dx,dt,tau,error,maxu;
double x_start,x_end;
double a;
int i,j,ip;
void init();
double feq(int k,double u,double ut);
void evolution(double t);
void output(int m);
void Error();
double exact_solution(double x,double t);
double exact_solution_ut(double x,double t);
inline double fai_one(double u);
inline double sour(double u);
inline double sour_t(double u,double ut);
inline double sour_tt(double u,double ut,double utt);
void u_time(double t);
void prepare();
void max();
int main()
{
int m=1;
int n=1;
double t=0;
x_start=-2.0;
x_end=2.0;
dx=(x_end-x_start)/NX;
c=1.0;
dt=dx/c;
ostringstream name;
name<<"GRE for different tau"<<".dat";
ofstream out(name.str().c_str());
out<<"dx="<<dx<<",dt="<<dt<<",c="<<c<<",a="<<a<<endl;
out<<"s[Q]="<<s[0]<<","<<s[1]<<","<<s[2]<<","<<s[3]<<endl;
tau=0.5;
for(m=1;m<=20;m++)
{
tau=tau+0.001;
t=0;
init();
for(n=1;n<=2000;n++)
{
t=n*dt;
evolution(t);
if(n%200==0)
{
u_time(t);
max();
Error();
out<<tau<<"t"<<n<<"t"
<<resetiosflags(ios::fixed)<<setiosflags(ios::scientific)<<setprecision(6)
<<error<<"t"<<maxu<<resetiosflags(ios::scientific)<<endl;
}
}
}
return 0;
}
void init()
{
a=1.0;
double x=0;
double t0=0;
u_time(t0);
for(i=0;i<=NX+2;i++)
{
x=x_start+double(i-1)*dx;
u[i]=exact_solution(x,t0);
u1[i]=u[i];
u2[i]=u[i];
ut0[i]=4.0/cosh(x);
utt[i]=0.0;
for(j=0;j<Q;j++)
{
f[i][j]=feq(j,u[i],ut[i]);
F[i][j]=0.0;
}
}
output(0);
}
void prepare()
{
for(i=0;i<=NX+2;i++)
{
g[i]=-a*a*u[i]*sin(u[i]);
}
for(i=2;i<=NX;i++)
{
sour_xx[i]=(-g[i+2]+16.0*g[i+1]-30.0*g[i]+16.0*g[i-1]-g[i-2])/(12.0*dx*dx);
}
sour_xx[0]=(-35.0*g[4]+104.0*g[3]-114.0*g[2]+56.0*g[1]-11.0*g[0])/(12.0*dx*dx);
sour_xx[1]=(-35.0*g[5]+104.0*g[4]-114.0*g[3]+56.0*g[2]-11.0*g[1])/(12.0*dx*dx);
sour_xx[NX+1]=(35.0*g[NX+1]-104.0*g[NX]+114.0*g[NX-1]-56.0*g[NX-2]+11.0*g[NX-3])/(12.0*dx*dx);
sour_xx[NX+2]=(35.0*g[NX+2]-104.0*g[NX+1]+114.0*g[NX]-56.0*g[NX-1]+11.0*g[NX-2])/(12.0*dx*dx);
for(i=0;i<=NX+2;i++)
{
ut[i]=(u[i]-u1[i])/dt;
utt[i]=(u[i]+u2[i]-2.0*u1[i])/(dt*dt);
}
}
void evolution(double t)
{
prepare();
int i=0;
int j=0;
int ip=0;
for(i=2;i<NX+1;i++)
{
for(j=0;j<Q;j++)
{
double f_s1,f_s2,f_s3,f_s4;
ip=i-e[j];
f_s1=dt*sour(u[ip])*s[j];
f_s2=dt*dt*sour_t(u[ip],ut[ip])*0.5*s[j];
f_s3=dt*dt*dt*sour_tt(u[ip],ut[ip],utt[ip])*s[j]/6.0;
f_s4=dt*dt*(2.0*tau*tau-2.0*tau+0.25)*sour_xx[ip]*s[j]/(tau-0.5);
F[i][j]=f[ip][j]+(feq(j,u[ip],ut[ip])-f[ip][j])/tau+f_s1+f_s2+f_s3+f_s4;
}
}
for(i=0;i<=NX+2;i++)
{
u2[i]=u1[i];
u1[i]=u[i];
}
for(i=2;i<NX+1;i++)
{
ut[i]=0.0;
for(j=0;j<Q;j++)
{
f[i][j]=F[i][j];
ut[i]+=f[i][j];
}
}
ut[1]=exact_solution_ut(x_start,t);
ut[0]=exact_solution_ut(x_start-dx,t);
ut[NX+1]=exact_solution_ut(x_end,t);
ut[NX+2]=exact_solution_ut(x_end+dx,t);
f[1][0]=feq(0,u[1],ut[1])+(f[2][0]-feq(0,u[2],ut[2]));
f[1][2]=feq(2,u[1],ut[1])+(f[2][2]-feq(2,u[2],ut[2]));
f[0][2]=feq(2,u[0],ut[0])+(f[1][2]-feq(2,u[1],ut[1]));
f[NX+1][1]=feq(1,u[NX+1],u[NX+1])+(f[NX][1]-feq(1,u[NX],ut[NX]));
f[NX+1][3]=feq(3,u[NX+1],u[NX+1])+(f[NX][3]-feq(3,u[NX],ut[NX]));
f[NX+2][3]=feq(3,u[NX+2],u[NX+2])+(f[NX+1][3]-feq(3,u[NX+1],ut[NX+1]));
}
double feq(int k,double u,double ut)
{
double feq=0;
switch(k)
{
case 0:feq=(4.0*ut-fai_one(u))/6.0;
return feq;
case 1:feq=(4.0*ut-fai_one(u))/6.0;
return feq;
case 2:feq=(fai_one(u)-ut)/6.0;
return feq;
case 3:feq=(fai_one(u)-ut)/6.0;
return feq;
default:
return 0;
}
}
void output(int m)
{
double x=0;
int i=0;
ostringstream name;
name<< "D1Q4_"<<m<<".dat";
ofstream out(name.str().c_str());
out << "Title = "heat conduction ep"n"
<< "VARIABLES = "X","U","U0"n"
<< "tau="<<tau<<" "<<NX+1<<"POINTs"<<" "<<"dx="<<dx<<" "<<"dt="<<dt<<" "<<endl;
for(i=0;i<=NX+2;i++)
{
x=x_start+double(i-1)*dx;
out << x << "t" <<u[i]<<"t"<<u0[i]<<endl;
}
}
void Error()
{
double temp1,temp2;
temp1=0;
temp2=0;
for(i=1;i<=NX+1;i++)
{
temp1+=fabs(u[i]-u0[i]);
temp2+=fabs(u0[i]);
}
error=temp1/temp2;
}
double exact_solution(double x,double t)
{
double u;
u=4.0*atan(t/cosh(x));
return u;
}
double exact_solution_ut(double x,double t)
{
double u;
u=4.0/(cosh(x)+t*t/cosh(x));
return u;
}
inline double fai_one(double u)
{
double fai_one;
fai_one=u*pow(a,2.0)/(dt*(tau-0.5)*c*c);
return fai_one;
}
inline double sour(double u)
{
double sour;
sour=-sin(u);
return sour;
}
inline double sour_t(double u,double ut)
{
double sour_t;
sour_t=-ut*cos(u);
return sour_t;
}
inline double sour_tt(double u,double ut,double utt)
{
double sour_tt;
sour_tt=-utt*cos(u)-ut*sin(u)*ut;
return sour_tt;
}
void u_time(double t)
{
double x=0;
for(i=0;i<=NX+2;i++)
{
x=x_start+double(i-1)*dx;
u0[i]=exact_solution(x,t);
}
}
void max()
{
double abu;
maxu=0;
for(i = 0;i <= NX+2 ;i++)
{
abu=fabs(u[i]-u0[i]);
if (abu>maxu)
maxu=abu;
}
}
输出的解为0，请问怎么修改
