一个求解偏微分方程问题的C程序.docVIP

一个求解偏微分方程问题的C程序 #include stdio.h#include math.h #define vareps 0.02#define gamma 0.1#define pi 3.1415926535898#define eps 1.0e-12#define x1 0#define x2 2*pi//参数取值 double norm(double*a,int N,double h){double b=0；int i,j；for(i=0；i N；i++)for(j=0；j N；j++)b+=a[i][j]*a[i][j]*h*h；b=sqrt(b)；return b；}//误差计算，其中空间网格数N，步长h void main(){double*u0,*u,*uh,*e,*x,*y,r,h,tau,T；int i,j,k,K,N；FILE*fp=fopen(data.m,w)；FILE*fp0=fopen(data1.m,w)； printf(Printf T：\n)；scanf(%lf,T)；//输入时间参数T printf(Printf the length of tau：\n)；scanf(%lf,tau)；//时间步长tau printf(Printf N：\n)；scanf(%d,N)；//空间网格数N K=(int)(T/tau)+1；h=(x2-x1)/N；//空间步长h N++；r=tau/(h*h)；//网比r uh=new double*[N]；for(i=0；i N；i++)uh[i]=new double[N]； u0=new double*[N]；for(i=0；i N；i++)u0[i]=new double[N]； u=new double*[N]；for(i=0；i N；i++)u[i]=new double[N]； e=new double*[N]；for(i=0；i N；i++)e[i]=new double[N]； x=new double[N]；y=new double[N]；for(i=0；i N；i++){x[i]=x1+i*h；y[i]=x[i]；} for(i=0；i N；i++)for(j=0；j N；j++)if((x[i]-pi)*(x[i]-pi)+(y[j]-pi)*(y[j]-pi)0.25)u[i][j]=1；else u[i][j]=-1； fprintf(fp0,x=[)；for(i=0；i N；i++)fprintf(fp0,%lf\t,x[i])； fprintf(fp0,]；\ny=[)；for(i=0；i N；i++)fprintf(fp0,%lf\t,y[i])； fprintf(fp0,]；\nz=[)；for(i=0；i N；i++){for(j=0；j N；j++)fprintf(fp0,%e\t,u[i][j])；fprintf(fp0,\n)；}fprintf(fp0,]；\nfigure,surf(x,y,z)\nfigure,contour(x,y,z)\n)；//输出初始条件 for(k=1；k K；k++){for(i=0；i N；i++)for(j=0；j N；j++)uh[i][j]=u[i][j]； while(1){for(i=0；i N；i++)for(j=0；j N；j++)u0[i][j]=u[i][j]； u[0][0]=(uh[0][0]/tau/gamma+2/h/h*(u[0][1]+u[1][0])-u[0][0]*u[0][0]*u[0][0]/vareps/vareps)/(1/tau/gamma+4/h/h-1/vareps/vareps)；u[N-1][0]=(uh[N-1][0]/tau/gamma+2/h/h*(u[N-1][1]+u[N-2][0])-u[N-1][0]*u[N-1][0]*u[N-1][0]/vareps/vareps)/(1/tau/gamma+4/h/h-1/vareps/vareps)；u[0][N-1]=(uh[0][N-1]/tau/gamma+2/h/h*(u[0][N-2]+u[1][N-1])-u[0][N-1]*u[0][N-1]*u[0][N-1]/vareps/vareps)/(1/tau/gamma+4/h/h-1/vareps/vareps)；u[N-1][N-1]=(uh[N-1][N-1]/tau/gamma+2/h/h*(u[N-2][N-1]+u[N-1][N-2])-u[N-1][N-1]*u[N-1][N-1]*u[N-1][N-1]/vareps/vareps)/(1/tau/gamma+4/h/h-1/vareps/varep