偏微分方程数值解法.doc

个人收集整理 勿做商业用途一、 问题用有限元方法求下面方程的数值解 in on in 二、 问题分析第一步 利用Green公式，求出方程的变分形式变分形式为：求,使得 （）以及 .第二步 对空间进行离散，得出半离散格式对区域进行剖分，构造节点基函数，得出有限元子空间：,则（)的Galerkin逼近为:，求，使得 （*）以及,为初始条件在中的逼近，设为在中的插值.则，有,=,代人（）即可得到一常微分方程组。第三步 进一步对时间进行离散，得到全离散的逼近格式对用差分格式.为此把等分为n个小区间，其长度,。这样把求时刻的近似记为,是的近似。这里对（*）采用向后的欧拉格式，即 (*)i=0，1，2,n1。 =由于向后欧拉格式为隐式格式且含有非线性项，故相邻两时间步之间采用牛顿迭代，即：其中用作为牛顿迭代的初值.三、 计算流程取作为方程的准确解，算得初值函数：右端函数：+。取。1 对作三角剖分：分别沿x , y方向对0，1区间作N等分并将每个小正方形沿对角线分为两个三角形。2 对节点进行总体编码（这里按先沿x方向再沿y方向进行顺序编码），并建立局部编码和总体编码之间的对应关系。3 计算各个节点的实际坐标，找出边界点。4 构造单元形状函数，计算雅可比矩阵并利用高斯求积公式计算单元刚度矩阵和单元列阵.5 合成总刚度矩阵和总列阵。6 处理本质边界条件，求解线性方程组.7 在第n个时间步利用牛顿迭代法算得第n+1个时间步的数值解，取.8 计算各个时间步的有限元数值解和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.hinclude ”stdlib.h”#include math.h”#define N 25 /沿x y方向的等分数define ND 3 /一个元素的节点个数#define LEE 2*N*N /元素总数define NG (N+1）*（N+1) /节点总数define TSTP 0。01 /时间步长define TN 100 /时间迭代步数define J 1.0/(NN) /雅可比行列式的绝对值 double u0（double x，double y) /初值函数u0 double z； z=100*xy(x1)(y1）； return z； double f(double x,double y,double t）/右端函数f double z； z=100xy（x-1）(y1)-200(t+1)y*(y-1)-200(t+1)*x(x-1)+10000（t+1)（t+1）xxyy*(x-1）*（x-1）(y-1）(y1); return z; double u（double x,double y，double t）/方程的准确解 double z； z=100*（t+1)*x*y*(x1）（y-1）; return z；void II（int *a) /节点的局部编码与总体编码 int i； for（i=1;iLEE+1；i+） if（i2=1） ai1=（i+1）/2+i/（2N); ai2=ai1+1; ai3=ai1+N+1； if(i2=0） ai3=i/2+1+(i-1）/(2*N）; ai1=ai3+N+1; ai2=ai3+N; struct xy double x;double y;；void XY(struct xy b) /节点实际坐标 int i； for(i=1;iNG+1;i+） if(i%(N+1）=0） bi.x=1； bi.y=(i/(N+1)-1)*（1.0/N); else bi.x=（i%(N+1）-1)*（1.0/N)； bi.y=（i/(N+1)（1。0/N）； void boundnote(int bd，struct xy b)/找出边界点int i，j=1；for（i=1；iNG+1;i+）if(bi.x=0|bi.x=1.0|bi.y=0|bi。y=1.0）bdj=i;j+；void dealwithbd(double *uk，int bd) /近似处理本质边界条件int i;for（i=bd1；bdi！=0；i+）ukbdibdi=ukbdibdi*10e20;void GAUSSELIMINATE（double *E,double RHSNG+1）/利用高斯消元法解线性方程组int i;int k;for(k=1;kNG；k+)/ 选主元double bmax=0.0;int ik;for(i=k；iNG+1；i+)if(bmaxfabs（Eik）)bmax=fabs（Eik）;ik=i;if（bmax1。0e-10) printf(”主元太小”); / 主元太小/ 交换第ik行和第k行的元素if(ik！=k)double t;for（i=k;iNG+1；i+)t=Eiki；Eiki=Eki；Eki=t；t=RHSik；RHSik=RHSk；RHSk=t;/ 消元for(i=k+1;iNG+1;i+)if(Eik!=0）double lk=Eik/Ekk；int j；for(j=k+1；j0;i-）double s=0。0;int j;for(j=i+1;jNG+1；j+）s=s+EijRHSj；RHSi=(RHSi-s）/Eii；void AIJ(double *aij，double aryk，int a) /计算单元刚度矩阵int i;for（i=1；iLEE+1；i+) aiji1=1.0/（12*N*N）+TSTP+2*TSTP1。0/（54*NN）*（arykai1+arykai2+arykai3)；/ 1 1aiji2=1.0/（12N*N)+0。5TSTP+2TSTP*1.0/(54N*N)*（arykai1+arykai2+arykai3）；/ 2 2aiji3=1.0/（12*NN)+0。5TSTP+2TSTP*1。0/（54*NN)（arykai1+arykai2+arykai3);/ 3 3aiji4=1。0/（24*NN)+(0.5TSTP)+2TSTP1。0/(54*NN）*（arykai1+arykai2+arykai3);/ 1 2 aiji5=1.0/（24*N*N）+(-0.5TSTP)+2*TSTP1.0/（54*N*N）*（arykai1+arykai2+arykai3);/ 1 3 aiji6=1.0/(24*N*N)+2TSTP1.0/（54*NN）（arykai1+arykai2+arykai3）；/ 2 3 void UK（double *uk，int *a,double *aij）/合成总刚度矩阵 int i;for(i=1；iLEE+1;i+)ukai1ai1+=aiji1; ukai2ai2+=aiji2; ukai3ai3+=aiji3； ukai1ai2+=aiji4;ukai2ai1+=aiji4； ukai1ai3+=aiji5； ukai3ai1+=aiji5； ukai2ai3+=aiji6； ukai3ai2+=aiji6;void FE(double *fe，double aryn，double aryk，int n,int *a,struct xy b） /计算单元列阵int i；for(i=1;iLEE+1;i+）fei1=1。0/(18N*N）（arynai1+arynai2+arynai3）+TSTP*1.0/（54*NN）*(arykai1+arykai2+arykai3)（arykai1+arykai2+arykai3)+TSTP1.0/(6*NN)f（(bai1.x+bai2.x+bai3.x)/3.0，(bai1.y+bai2.y+bai3。y）/3。0,(n+1)*TSTP）;void UFE(double ufe，double *fe，int *a)/合成总列阵int i;for（i=1；iLEE+1；i+）ufeai1+=fei1; ufeai2+=fei1; ufeai3+=fei1;void NEWTONITERATIVE（double aryn1，double aryn,int n，int *a，struct xy b，int bd）/牛顿迭代int i;double *aryk;aryk=(double )calloc(NG+1，sizeof（double);for(i=1；iNG+1；i+） aryki=aryni; double boot； double *uk;double *ufe；double *aij; double fe;douk=（double *)calloc(NG+1,sizeof(double *);if（uk=NULL) /判断内存分配是否成功printf(uk failn);exit(1）; for（i=0；iNG+1；i+) uki=（double *)calloc(NG+1,sizeof(double）); if(uki=NULL) printf(uki failn)； exit（1)；ufe=(double *）calloc（NG+1，sizeof(double)；if（ufe=NULL）printf(ufe failn”);exit（1)； aij=（double ）calloc(LEE+1，sizeof(double ）；if(aij=NULL）printf（aij failn）;exit(1）; for(i=0;iLEE+1;i+） aiji=(double )calloc（7，sizeof（double）; if（aiji=NULL)printf（aiji failn”);exit(1）； fe=（double *）calloc（LEE+1,sizeof（double );if（fe=NULL）printf(fe failn）;exit（1）； for(i=0；iboot）boot=fabs(ufei-aryki)；for（i=1;iNG+1;i+） aryki=ufei;for(i=0；iNG+1；i+）free(uki)；free(uk）;free（ufe）； for（i=0;iLEE+1；i+）free(aiji)；free(aij);for（i=0;iLEE+1;i+）free（fei）；free(fe）；while（boot1.0e-8); for(i=1;iNG+1;i+） aryn1i=aryki；free（aryk);/*以下为主函数*/main(） int i，j;/节点的局部编码与总体编码 int a; a=(int *)calloc（LEE+1，sizeof(int ); for（i=0；iLEE+1；i+） ai=（int )calloc（ND+1，sizeof（int)）; II(a);/节点实际坐标 struct xy *b； b=（struct xy )calloc（NG+1，sizeof(struct xy）； XY(b);/边界点 int bd； bd=(int *)calloc(NG+1，sizeof(int）)； boundnote（bd,b);/计算各个时间步的数值解 double fesolution； fesolution=(double *)calloc(NG+1，sizeof(double）)； int n； double aryn； aryn=（double *）calloc（NG+1,sizeof（double）)； for（i=1；iNG+1;i+) aryni=u0（bi。x，bi.y); double *aryn1； for(n=0;nTN;n+） aryn1=(double *）calloc（NG+1,sizeof（double）; NEWTONITERATIVE(aryn1,aryn,n,a，b，bd); for(i=0；iNG+1；i+)aryni=aryn1i;free（aryn1）; for（i=1；iNG+1；i+) fesolutioni=aryni; printf（”n 第 d 步的结果为：nn,n); printf（节点总体编码 有限元值 精确值 误差 nn”）； for(j=N/2；jNG+1;j+=2*（N+1） printf(%8d 15f %15f 15fn,j，fesolutionj,u（bj。x,bj。y，TSTP*TN）,fesolutionju（bj。x,bj.y，TSTPTN）)；/计算L2范数意义下的误差 double L2=0； for(i=1;iLEE+1；i+) L2+=J1.0/2*pow(（u(bai1。x+bai2.x+bai3。x)/3。0,(bai1.y+bai2。y+bai3。y)/3。0，n*TSTP)1。0/3*(fesolutionai1+fesolutionai2+fesolutionai3）),2）; L2=sqrt(L2）； printf(”n 第 %d 步的L2范数意义下的误差为： ”,n）； printf（%fnn，L2）； return 0;