偏微分方程数值解法.doc

个人收集整理 勿做商业用途

一、 问题
用有限元方法求下面方程的数值解
in
on
in
二、 问题分析
第一步 利用Green公式，求出方程的变分形式
变分形式为：求,使得
（＊）
以及 .
第二步 对空间进行离散，得出半离散格式
对区域进行剖分，构造节点基函数，得出有限元子空间：,则（＊)的Galerkin逼近为:
，求，使得
（＊*）
以及,为初始条件在中的逼近，设为在中的插值.
则，有,代人（＊＊）即可得到一常微分方程组。
第三步 进一步对时间进行离散，得到全离散的逼近格式
对用差分格式.为此把等分为n个小区间，其长度,。
这样把求时刻的近似记为,是的近似。这里对（*＊）采用向后的欧拉格式，即
(**＊)
i=0，1，2…,n—1。 =
由于向后欧拉格式为隐式格式且含有非线性项，故相邻两时间步之间采用牛顿迭代，即：

其中用作为牛顿迭代的初值.
三、 计算流程
取作为方程的准确解，算得初值函数：

右端函数：
+
。
取。
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.h"
＃include ”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/(N＊N) //雅可比行列式的绝对值

double u0（double x，double y) //初值函数u0
｛
double z；
z=100*x＊y＊(x—1)＊(y—1）；
return z；
｝

double f(double x,double y,double t）//右端函数f
｛
double z；
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);
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｝；

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)；
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）
{
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]）//利用高斯消元法解线性方程组
｛
int i;
int k;
for(k=1;k〈NG；k++)
{
// 选主元
double bmax=0.0;
int ik;
for(i=k；i太小
｝
// 交换第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(E[i］［k］!=0）
｛
double lk=E［i］[k］/E[k］［k］；
int j；
for(j=k+1；j0;i—-）
｛
double s=0。0;
int j;
for(j=i+1;j〈NG+1；j++）
｛
s=s+E[i]［j］＊RHS［j］；
}
RHS[i］=(RHS[i］-s）/E［i][i］；
}
｝

void AIJ(double ＊*aij，double aryk[］，int ＊＊a) //计算单元刚度矩阵
{
int i;
for（i=1；ia[i]［1]］+aryk［a［i］[2]］+aryk［a[i］［3]］)；// 1 1
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
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］］);// 1 2
aij［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］］);// 1 3
aij［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]］）；// 2 3
}
}

void UK（double ＊*uk，int **a,double ＊*aij）//合成总刚度矩阵
｛
int i;
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];
｝

}

void FE(double **fe，double aryn［]，double aryk[］，int n,int ＊*a,struct xy b［］） //计算单元列阵
{
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）;
｝

void UFE(double ufe［］，double **fe，int ＊*a)//合成总列阵
{
int i;
for（i=1；iaryk;
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;

do
{
uk=（double **)calloc(NG+1,sizeof(double *));
if（uk==NULL) //判断内存分配是否成功
{
printf("uk fail\n");
exit(1）;
}

for（i=0；i1)；
｝
aij=（double ＊＊）calloc(LEE+1，sizeof(double ＊））；
if(aij==NULL）
{
printf（"aij fail\n"）;
exit(1）;
｝
for(i=0;i〈LEE+1;i++）
｛
aij[i］=(double ＊)calloc（7，sizeof（double））;
if（aij[i]==NULL)
｛
printf（"aiji fail\n”);
exit(1）；
｝
｝
fe=（double *＊）calloc（LEE+1,sizeof（double ＊));
if（fe==NULL）
｛
printf("fe fail\n"）;
exit（1）；
｝
for(i=0；ij，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＊**＊＊＊*＊*/

main(）
{
int i，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);

//节点实际坐标

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；iTN;n++）
｛

aryn1=(double *）calloc（NG+1,sizeof（double））;

NEWTONITERATIVE(aryn1,aryn,n,a，b，bd);

for(i=0；i 误差 \n\n”）；
for(j=N/2；j〈NG+1;j+=2*（N+1））
printf("%8d ％15f %15f ％15f\n",j，fesolution[j],u（b[j]。x,b［j］。y，TSTP*TN）,fesolution［j]—u（b［j］。x,b[j].y，TSTP＊TN）)；

//计算L2范数意义下的误差

double L2=0；
for(i=1;i〈LEE+1；i++)
L2+=J＊1.0/2*pow(（u((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*TSTP)—1。0/3*(fesolution［a［i][1]]+fesolution[a［i]［2]］+fesolution［a［i］[3］］）),2）;
L2=sqrt(L2）；
printf(”\n 第 %d 步的L2范数意义下的误差为： ”,n）；
printf（"%f\n\n"，L2）；

return 0;
}