偏微分方程数值解法.doc

个人收集整理勿做商业用途个人收集整理勿做商业用途个人收集整理勿做商业用途一、问题用有限元方法求下面方程的数值解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（do