偏微分方程数值解法

一、问题用有限元方法求下面方程的数值解inonin二、问题分析第一步利用Green公式，求出方程的变分形式变分形式为：求，使得（*）以及.第二步对空间进行离散，得出半离散格式对区域进行剖分，构造节点基函数，得出有限元子空间：，则（*）的Galerkin逼近为：，求，使得（**）以及，为初始条件在中的逼近，设为在中的插值.则，有，=，代人（**）即可得到一常微分方程组.第三步进一步对时间进行离散，得到全离散的逼近格式对用差分格式.为此把等分为n个小区间，其长度,.这样把求时刻的近似记为，是的近似.这里对（**）采用向后的欧拉格式，即(***)i=0,1,2…,n-1.=由于向后欧拉格式为隐式格式且含有非线性项，故相邻两时间步之间采用牛顿迭代，即：其中用作为牛顿迭代的初值.三、计算流程取作为方程的准确解，算得初值函数：右端函数：+.取.对作三角剖分：分别沿x,y方向对[0,1]区间作N等分并将每个小正方形沿对角线分为两个三角形.对节点进行总体编码（这里按先沿x方向再沿y方向进行顺序编码），并建立局部编码和总体编码之间的对应关系.计算各个节点的实际坐标，找出边界点.构造单元形状函数，计算雅可比矩阵并利用高斯求积公式计算单元刚度矩阵和单元列阵.合成总刚度矩阵和总列阵.处理本质边界条件，求解线性方程组.在第n个时间步利用牛顿迭代法算得第n+1个时间步的数值解，取.计算各个时间步的有限元数值解和L2误差.四、误差分析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阶数为附C程序#include"stdio.h"#include"stdlib.h"#include"math.h"#defineN25//沿xy方向的等分数#defineND3//一个元素的节点个数#defineLEE2*N*N//元素总数#defineNG(N+1)*(N+1)//节点总数#defineTSTP0.01//时间步长#defineTN100//时间迭代步数#defineJ1.0/(N*N)//雅可比行列式的绝对值doubleu0(doublex,doubley)//初值函数u0{doublez;z=100*x*y*(x-1)*(y-1);returnz;}doublef(doublex,doubley,doublet)//右端函数f{doublez;z=100*x*y*(x-1)*(y-1)-200*(t+1)*y*(y-1)-200*(t+1)*x*(x-1)+10000*(t+1)*(t+1)*x*x*y*y*(x-1)*(x-1)*(y-1)*(y-1);returnz;}doubleu(doublex,doubley,doublet)//方程的准确解{doublez;z=100*(t+1)*x*y*(x-1)*(y-1);returnz;}voidII(int**a)//节点的局部编码与总体编码{inti;for(i=1;i<LEE+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; }}}structxy{doublex;doubley;};voidXY(structxyb[])//节点实际坐标{ inti; 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); b[i].y=(i/(N+1))*(1.0/N); }}voidboundnote(intbd[],structxyb[])//找出边界点{ inti,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) { bd[j]=i; j++; } }}voiddealwithbd(double**uk,intbd[])//近似处理本质边界条件{ inti; for(i=bd[1];bd[i]!=0;i++) uk[bd[i]][bd[i]]=uk[bd[i]][bd[i]]*10e20;}voidGAUSSELIMINATE(double**E,doubleRHS[NG+1])//利用高斯消元法解线性方程组{ inti; intk; for(k=1;k<NG;k++) { //选主元 doublebmax=0.0; intik; 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("主元太小"); //主元太小 } //交换第ik行和第k行的元素 if(ik!=k) { doublet; 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(E[i][k]!=0) { doublelk=E[i][k]/E[k][k]; intj; for(j=k+1;j<NG+1;j++) { E[i][j]=E[i][j]-lk*E[k][j]; } RHS[i]=RHS[i]-lk*RHS[k]; } } } if(fabs(E[NG][NG])<1.0e-10) { printf("主元太小"); //主元太小 } //开始回代 RHS[NG]=RHS[NG]/E[NG][NG]; for(i=NG-1;i>0;i--) { doubles=0.0; intj; for(j=i+1;j<NG+1;j++) { s=s+E[i][j]*RHS[j]; } RHS[i]=(RHS[i]-s)/E[i][i]; }}voidAIJ(double**aij,doublearyk[],int**a)//计算单元刚度矩阵{ inti; 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]]);//11 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]]);//22 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]]);//33 aij[i][4]=1.0/(24*N*N)+(-0.5*TSTP)+2*TSTP*1.0/(54*N*N)*(aryk[a[i][1]]+aryk[a[i][2]]+aryk[a[i][3]]);//12aij[i][5]=1.0/(24*N*N)+(-0.5*TSTP)+2*TSTP*1.0/(54*N*N)*(aryk[a[i][1]]+aryk[a[i][2]]+aryk[a[i][3]]);//13aij[i][6]=1.0/(24*N*N)+2*TSTP*1.0/(54*N*N)*(aryk[a[i][1]]+aryk[a[i][2]]+aryk[a[i][3]]);//23 }}voidUK(double**uk,int**a,double**aij)//合成总刚度矩阵{inti; 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];uk[a[i][1]][a[i][3]]+=aij[i][5];uk[a[i][3]][a[i][1]]+=aij[i][5];uk[a[i][2]][a[i][3]]+=aij[i][6];uk[a[i][3]][a[i][2]]+=aij[i][6]; } }voidFE(double**fe,doublearyn[],doublearyk[],intn,int**a,structxyb[])//计算单元列阵{ inti; 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);}voidUFE(doubleufe[],double**fe,int**a)//合成总列阵{ inti; for(i=1;i<LEE+1;i++) { ufe[a[i][1]]+=fe[i][1];ufe[a[i][2]]+=fe[i][1];ufe[a[i][3]]+=fe[i][1]; }}voidNEWTONITERATIVE(doublearyn1[],doublearyn[],intn,int**a,structxyb[],intbd[])//牛顿迭代{ inti; double*aryk; aryk=(double*)calloc(NG+1,sizeof(double)); for(i=1;i<NG+1;i++) aryk[i]=aryn[i];doubleboot;double**uk; double*ufe; double**aij;double**fe; do { uk=(double**)calloc(NG+1,sizeof(double*)); if(uk==NULL)//判断内存分配是否成功 { printf("ukfail\n"); exit(1); }for(i=0;i<NG+1;i++) { uk[i]=(double*)calloc(NG+1,sizeof(double)); if(uk[i]==NULL) { printf("ukifail\n"); exit(1); } } ufe=(double*)calloc(NG+1,sizeof(double)); if(ufe==NULL) { printf("ufefail\n"); exit(1); }aij=(double**)calloc(LEE+1,sizeof(double*)); if(aij==NULL) { printf("aijfail\n"); exit(1); }for(i=0;i<LEE+1;i++) { aij[i]=(double*)calloc(7,sizeof(double));if(aij[i]==NULL) { printf("aijifail\n"); exit(1); } }fe=(double**)calloc(LEE+1,sizeof(double*)); if(fe==NULL) { printf("fefail\n"); exit(1); }for(i=0;i<LEE+1;i++) { fe[i]=(double*)calloc(2,sizeof(double)); if(fe[i]==NULL) { printf("feifail\n"); exit(1); } } AIJ(aij,aryk,a);UK(uk,a,aij);FE(fe,aryn,aryk,n,a,b);UFE(ufe,fe,a);dealwithbd(uk,bd); GAUSSELIMINATE(uk,ufe); boot=0; for(i=1;i<NG+1;i++) if(fabs(ufe[i]-aryk[i])>boot)boot=fabs(ufe[i]-aryk[i]); for(i=1;i<NG+1;i++) aryk[i]=ufe[i]; for(i=0;i<NG+1;i++) free(uk[i]); free(uk); free(ufe);for(i=0;i<LEE+1;i++) free(aij[i]); free(aij); for(i=0;i<LEE+1;i++) free(fe[i]); free(fe); } while(boot>1.0e-8);for(i=1;i<NG+1;i++) aryn1[i]=aryk[i]; free(aryk);}/***********以下为主函数*******************/main(){inti,j;//节点的局部编码与总体编码int**a;a=(int**)calloc(LEE+1,sizeof(int*));for(i=0;i<LEE+1;i++) a[i]=(int*)calloc(ND+1,sizeof(int));II(a);//节点实际坐标structxy*b;b=(structxy*)calloc(NG+1,sizeof(structxy));XY(b);//边界点int*bd;bd=(int*)calloc(NG+1,sizeof(int));boundnote(bd,b);//计算各个时间步的数值解double*fesolution;fesolution=(double*)calloc(NG+1,sizeof(double));intn;double*aryn;aryn=(double*)calloc(NG+1,sizeof(double));for(i=1;i<NG+1;i++) a