b[i].y=(i/(N+1))*(1.0/N);
}
}
void boundnote(int bd[],struct xy b[])//找出边界点
{
int i,j=1;
for(i=1;i<NG+1;i++)
{
if(b[i].x==0||b[i].x==1.0||b[i].y==0||b[i].y==1.0)
#define TSTP 0.01 //时间步长
#define TN 100 //时间迭代步数
#define J 1.0/(N*N) //雅可比行列式的绝对值
double u0(double x,double y) //初值函数u0
{
double z;
z=100*x*y*(x-1)*(y-1);
return z;
}
}
void AIJ(double **aij,double aryk[],int **a) //计算单元刚度矩阵
{
int i;
for(i=1;i<LEE+1;i++)
{
aij[i][1]=1.0/(12*N*N)+TSTP+2*TSTP*1.0/(54*N*N)*(aryk[a[i][1]]+aryk[a[i][2]]+aryk[a[i][3]]);// 1 1
//主元太小
}
//交换第ik行和第k行的元素
if(ik!=k)
{
double t;
for(i=k;i<NG+1;i++)
{
t=E[ik][i];
E[ik][i]=E[k][i];
E[k][i]=t;
}
t=RHS[ik];
RHS[ik]=RHS[k];
RHS[k]=t;
}
//消元
for(i=k+1;i<NG+1;i++)
{if(i%2==1)
{
a[i][1]=(i+1)/2+i/(2*N);
a[i][2]=a[i][1]+1;
a[i][3]=a[i][1]+N+1;
}
if(i%2==0)
{
a[i][3]=i/2+1+(i-1)/(2*N);
a[i][1]=a[i][3]+N+1;
a[i][2]=a[i][3]+N;
}
void UFE(double ufe[],double **fe,int **a)//合成总列阵
{
int i;
for(i=1;Biblioteka <LEE+1;i++)
{
ufe[a[i][1]]+=fe[i][1];
ufe[a[i][2]]+=fe[i][1];
ufe[a[i][3]]+=fe[i][1];
}
}
void NEWTONITERATIVE(double aryn1[],double aryn[],int n,int **a,struct xy b[],int bd[])//牛顿迭代
其中用 作为牛顿迭代的初值.
三、计算流程
取 作为方程的准确解，算得初值函数：
右端函数：
+
.
取 .
1对 作三角剖分：分别沿x , y方向对[0,1]区间作N等分并将每个小正方形沿对角线分为两个三角形.
2对节点进行总体编码（这里按先沿x方向再沿y方向进行顺序编码），并建立局部编码和总体编码之间的对应关系.
L2误差 的阶数为 其中 为时间步长， 为空间步长
这里取N=25,则元素总数LEE = 2*N*N =1250,节点总数NG = (N+1)(N+1) = 676,
取时间步长TSTP = 0.01,时间步数TN = 100,即在t=1时，算得结果为：
即当 = 0.01, =0.04时，
L2误差为0.052853阶数为
for(i=1;i<LEE+1;i++)
{
uk[a[i][1]][a[i][1]]+=aij[i][1];
uk[a[i][2]][a[i][2]]+=aij[i][2];
uk[a[i][3]][a[i][3]]+=aij[i][3];
uk[a[i][1]][a[i][2]]+=aij[i][4];
uk[a[i][2]][a[i][1]]+=aij[i][4];
{
bd[j]=i;
j++;
}
}
}
void dealwithbd(double **uk,int bd[]) //近似处理本质边界条件
{
int i;
for(i=bd[1];bd[i]!=0;i++)
uk[bd[i]][bd[i]]=uk[bd[i]][bd[i]]*10e20;
}
void GAUSSELIMINATE(double **E,double RHS[NG+1])//利用高斯消元法解线性方程组
return z;
}
double u(double x,double y,double t)//方程的准确解
{
double z;
z=100*(t+1)*x*y*(x-1)*(y-1);
return z;
}
void II(int **a) //节点的局部编码与总体编码
{
int i;
for(i=1;i<LEE+1;i++)
{
int i;
for(i=1;i<LEE+1;i++)
fe[i][1]=1.0/(18*N*N)*(aryn[a[i][1]]+aryn[a[i][2]]+aryn[a[i][3]])+TSTP*1.0/(54*N*N)*(aryk[a[i][1]]+aryk[a[i][2]]+aryk[a[i][3]])*(aryk[a[i][1]]+aryk[a[i][2]]+aryk[a[i][3]])+TSTP*1.0/(6*N*N)*f((b[a[i][1]].x+b[a[i][2]].x+b[a[i][3]].x)/3.0,(b[a[i][1]].y+b[a[i][2]].y+b[a[i][3]].y)/3.0,(n+1)*TSTP);
{
int i;
int k;
for(k=1;k<NG;k++)
{
//选主元
double bmax=0.0;
int ik;
for(i=k;i<NG+1;i++)
{
if(bmax<fabs(E[i][k]))
{
bmax=fabs(E[i][k]);
ik=i;
}
}
if(bmax<1.0e-10)
{
printf("主元太小");
{
int i;
double *aryk;
aryk=(double *)calloc(NG+1,sizeof(double));
for(i=1;i<NG+1;i++)
aryk[i]=aryn[i];
double boot;
double **uk;
double *ufe;
double **aij;
double **fe;
}
}
}
struct xy {double x;double y;};
void XY(struct xy b[]) //节点实际坐标
{
int i;
for(i=1;i<NG+1;i++)
if(i%(N+1)==0)
{
b[i].x=1;
b[i].y=(i/(N+1)-1)*(1.0/N);
}
else
{
b[i].x=(i%(N+1)-1)*(1.0/N);
if(aij==NULL)
{
printf("aij fail\n");
aij[i][2]=1.0/(12*N*N)+0.5*TSTP+2*TSTP*1.0/(54*N*N)*(aryk[a[i][1]]+aryk[a[i][2]]+aryk[a[i][3]]);// 2 2
aij[i][3]=1.0/(12*N*N)+0.5*TSTP+2*TSTP*1.0/(54*N*N)*(aryk[a[i][1]]+aryk[a[i][2]]+aryk[a[i][3]]);// 3 3
if(uk[i]==NULL)
{
printf("uki fail\n");
exit(1);
}
}
ufe=(double *)calloc(NG+1,sizeof(double));
if(ufe==NULL)
{
printf("ufe fail\n");
exit(1);
}
aij=(double **)calloc(LEE+1,sizeof(double *));
}
double f(double x,double y,double t)//右端函数f
{
double z;