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

一个求解偏微分方程问题的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+=aij*aij*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+)uhi=new doubleN；u0=new double*N；for(i=0；i N；i+)u0i=new doubleN；u=new double*N；for(i=0；i N；i+)ui=new doubleN；e=new double*N；for(i=0；i N；i+)ei=new doubleN；x=new doubleN；y=new doubleN；for(i=0；i N；i+)xi=x1+i*h；yi=xi；for(i=0；i N；i+)for(j=0；j N；j+)if(xi-pi)*(xi-pi)+(yj-pi)*(yj-pi)0.25)uij=1；else uij=-1；fprintf(fp0,x=)；for(i=0；i N；i+)fprintf(fp0,%lft,xi)；fprintf(fp0,；ny=)；for(i=0；i N；i+)fprintf(fp0,%lft,yi)；fprintf(fp0,；nz=)；for(i=0；i N；i+)for(j=0；j N；j+)fprintf(fp0,%et,uij)；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+)uhij=uij；while(1)for(i=0；i N；i+)for(j=0；j N；j+)u0ij=uij；u00=(uh00/tau/gamma+2/h/h*(u01+u10)-u00*u00*u00/vareps/vareps)/(1/tau/gamma+4/h/h-1/vareps/vareps)；uN-10=(uhN-10/tau/gamma+2/h/h*(uN-11+uN-20)-uN-10*uN-10*uN-10/vareps/vareps)/(1/tau/gamma+4/h/h-1/vareps/vareps)；u0N-1=(uh0N-1/tau/gamma+2/h/h*(u0N-2+u1N-1)-u0N-1*u0N-1*u0N-1/vareps/vareps)/(1/tau/gamma+4/h/h-1/vareps/vareps)；uN-1N-1=(uhN-1N-1/tau/gamma+2/h/h*(uN-2N-1+uN-1N-2)-uN-1N-1*uN-1N-1*uN-1N-1/vareps/vareps)/(1/tau/gamma+4/h/h-1/vareps/vareps)；for(i=1；i N-1；i+)ui0=(uhi0/tau/gamma+1/h/h*(2*ui1+ui-10+ui+10)-ui0*ui0*ui0/vareps/vareps)/(1/tau/gamma+4/h/h-1/vareps/vareps)；uiN-1=(uhiN-1/tau/gamma+1/h/h*(2*uiN-2+ui-1N-1+ui+1N-1)-uiN-1*uiN-1*uiN-1/vareps/vareps)/(1/tau/gamma+4/h/h-1/vareps/vareps)；u0i=(uh0i/tau/gamma+1/h/h*(u0i+1+u0i-1+2*u1i)-u0i*u0i*u0i/vareps/vareps)/(1/tau/gamma+4/h/h-1/vareps/vareps)；uN-1i=(uhN-1i/tau/gamma+1/h/h*(uN-1i+1+uN-1i-1+2*uN-2i)-uN-1i*uN-1i*uN-1i/vareps/vareps)/(1/tau/gamma+4/h/h-1/vareps/vareps)；for(i=1；i N-1；i+)for(j=1；j N-1；j+)uij=(uhij/tau/gamma+1/h/h*(ui+1j+ui-1j+uij+1+uij-1)-uij*uij*uij/vareps/vareps)/(1/tau/gamma+4/h/h-1/vareps/vareps)；for(i=0；i N；i+)for(j=0；j N；j+)eij=fabs(uij-u0ij)；if(norm(e,N,h)eps)fprintf(fp,x=)；for(i=0；i N；i+)fprintf(fp,%lft,xi)；fprintf(fp,；ny=)；for(i=0；i N；i+)fprintf(fp,%lft,yi)；fprintf(fp,；nz=)；for(i=0；i N；i+