c语言偏微分方程函数包,[转载]一个求解偏微分方程问题的C程序

#include

#include

#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

for(j=0;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

uh[i]=new double[N];

u0=new double*[N];

for(i=0;i

u0[i]=new double[N];

u=new double*[N];

for(i=0;i

u[i]=new double[N];

e=new double*[N];

for(i=0;i

e[i]=new double[N];

x=new double[N];

y=new double[N];

for(i=0;i

{

x[i]=x1+i*h;y[i]=x[i];

}

for(i=0;i

for(j=0;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

fprintf(fp0,"%lft",x[i]);

fprintf(fp0,"];ny=[");

for(i=0;i

fprintf(fp0,"%lft",y[i]);

fprintf(fp0,"];nz=[");

for(i=0;i

{

for(j=0;j

fprintf(fp0,"%et",u[i][j]);

fprintf(fp0,"n");

}

fprintf(fp0,"];nfigure,surf(x,y,z)nfigure,contour(x,y,z)n");//输出初始条件

for(k=1;k

{

for(i=0;i

for(j=0;j

uh[i][j]=u[i][j];

while(1)

{

for(i=0;i

for(j=0;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/vareps);

for(i=1;i

{ u[i][0]=(uh[i][0]/tau/gamma+1/h/h*(2*u[i][1]+u[i-1][0]+u[i+1][0])

-u[i][0]*u[i][0]*u[i][0]/vareps/vareps)/(1/tau/gamma+4/h/

h-1/vareps/vareps);

u[i][N-1]=(uh[i][N-1]/tau/gamma+1/h/h*(2*u[i][N-2]+u[i-1][N-1]+

u[i+1][N-1])-u[i][N-1]*u[i][N-1]*u[i][N-1]/vareps/vareps)

/(1/tau/gamma+4/h/h-1/vareps/vareps);

u[0][i]=(uh[0][i]/tau/gamma+1/h/h*(u[0][i+1]+u[0][i-1]+2*u[1][i])-

u[0][i]*u[0][i]*u[0][i]/vareps/vareps)/(1/tau/gamma+4/h/h-

1/vareps/vareps);

u[N-1][i]=(uh[N-1][i]/tau/gamma+1/h/h*(u[N-1][i+1]+u[N-1][i-1]+2*