偏微分方程报告

1、成绩2009级数学与应用数学和信息与计算科学专业偏微分方程数值解上机实验实验题目利用有限兀方法和有限差分方法求解偏微分方程完成日期2012年12月17日学生姓名张灵刚所在班级1102090任课教师王晓东西北工业大学理学院应用数学系目录实验目的. （ 2）二实验要求.(2)三实验题目3）四实验二. ( 4)2.1. 实验内容 . ( 4)实验原理算法流程结果分析5 总结讨论6 源程序五. 实验三1. 实验内容2.实验原理. ( 4).5). ( 5). ( 6)6)17). (17) .17. ( 18)算法流程结果分析. .( 18) .5 总结讨论. ( 21)6 源程序21)偏微分方程数值

2、解上机实验报告 实验地点 ：数学系机房实验时间 ：第 1315 周，周一、四下午 5、6 节实验分数 ：占期末考试成绩的 30%一、实验目的及意义掌握有限元方法和有限差分方法的程序实现；学会选择合适 的有限差分格式求解一维非线性对流占优的非定常对流扩散问 题；学会使用三角线性元和四边形线性元的有限元方法求解二维 椭圆方程边值问题， 并对计算结果进行收敛性分析； 尝试采用有 限元方法或有限差分方法实现二维初边值抛物型方程的大规模 数值求解。 通过实验可以提高学生的动手能力， 加深学生对算法 的理解。二、实验要求在下列给出的三个问题中， 最少选择两个问题进行编程实现。 要求给出格式的推导过程、算法

3、流程、实现程序、选取的网格参 数、以表格或图形的方式给出计算结果、对计算结果进行分析、 最后对实验进行总结和讨论。问题2:用三角线性元和四边形线性元的有限元方法求解方程:u 8 2 cos(2 x)sin(2y),(x, y) (0,1) (0,1)u(x,0)u(x,1)0;u(0, y) u(1,y) sin(2 y).取h = 1/4, 1/8,1/16, 1/32, 1/64,比较两种方法的计算精度，并给出数值问题3 :选用合适的数值方法求解方程:(4 4 x2 4 y2) 1 u, (x,y) (0,1) (0,1)u(0,y, t) u(1,y,t)Uy (x,0,t)Uy(x,1

4、,t)0;u(x, y,0) sin( x2)cos( y2).47 49O.1,时点嚼,讣、越9煜鲁、（1281287137、 z 43 67、)、()、128128（詈毘、（环燼、（筹鬆）、（茲龛处的数值解。上机实验(二)一、实验内容用三角线性元和四边形线性元的有限元方法求解方程：2u 8 cos(2 x)sin(2 y), (x, y) (0,1) (0,1) u(x,0) u(x,1) 0; u(0, y) u(1,y) sin(2 y).取h = 1/4, 1/8,1/16, 1/32, 1/64,比较两种方法的计算精度，并给出数值 收敛率二、实验原理由于计算机只能存储有限个数据和做

5、有限次运算，所以必须把连续问题离散化，有限元法通过把求解偏微分问题转化为变分 问题，实质上就是Ritz-Galerkin法，利用有穷维空间近似代替无穷维空间，从而转化为求多元二次函数的极值问题。然后选定单元形状，对求解域进行剖分。构造基函数或单元形状函数，形 成有限元空间了，再形成有限元方程，并提供求解有限元方程的 有效方法。三、算法流程四、计算结果及分析空间步三角元实四边形元实际长h际误差误差140.1298160.027308180.0166760.0036591160.0020810.0004571320.0002600.000057丄640.0000320.000007结果分析：通过给

6、出的上述结果可以发现三种方法相比四边形元 的实际求解误差最小，这是由于四边形型元中的基函数中含有 x y项，使得四边形元的误差远比三角形元小，且当空间步长不 断加细时，向前差分的误差逐渐和四边形元接近。而对于收敛率的求解，如果舍去误差，三种方法的收敛率都约为3.五、总结和讨论有限元计算的有关问题有：把初值问题化为变分形式， 对求解域作网络分割，构造基函数（或单元形状函数），形成有限元方程。 通过本次实验， 我懂得了用有限元方法求解初值问题的基 本数学思想 , 理解了线性有限元法的基本原理和方法。另外，我 也懂得了按 Galerkin 方法推导有限元方程的优点， 它比 Ritz 法 更加方便直接

7、。 我也对虚功原理有了初步的认识。 因为 Galerkin 方法基于虚功原理， 所以不但可用于保守场问题， 也可使用于非 保守场即非驻定问题。附录：程序源代码1. 三角形元的型函数function phi, dxphi, dyphi=shape(point, element, ei, gx)x1=point(element(ei,1),1);x2=point(element(ei,2),1);x3=point(element(ei,3),1);y1=point(element(ei,1),2);y2=point(element(ei,2),2);y3=point(element(ei,3),2

8、);S=max(y2-y1,x2-x1)*max(y3-y1,x3-x1);ksi=(x3*y1-y3*x1)+(y3-y1)*gx(1)+(x1-x3)*gx(2)/S;eta=(x1*y2-y1*x2)+(y1-y2)*gx(1)+(x2-x1)*gx(2)/S;phi(1)=1-ksi-eta;phi(2)=ksi;phi(3)=eta;dxphi(1)=-(y3-y1)/S-(y1-y2)/S;dxphi(2)=(y3-y1)/S;dxphi(3)=(y1-y2)/S;dyphi(1)=-(x1-x3)/S-(x2-x1)/S;dyphi(2)=(x1-x3)/S;dyphi(3)=(

9、x2-x1)/S;2. 四边形元的型函数function phi, dxphi, dyphi=shape1(point, element, ei, gx)x1=point(element(ei,1),1);x2=point(element(ei,2),1);x3=point(element(ei,3),1);x4=point(element(ei,4),1);y1=point(element(ei,1),2);y2=point(element(ei,2),2);y3=point(element(ei,3),2);y4=point(element(ei,4),2);S=max(y2-y1,x2-

10、x1)*max(y4-y1,x4-x1);ksi=(gx(1)-x1)/(x2-x1);eta=(gx(2)-y1)/(y3-y1);phi(1)=(1-ksi)*(1-eta);phi(2)=(1-ksi)*eta;phi(3)=ksi*eta;phi(4)=ksi*(1-eta);dxphi(1)=-1/(x2-x1)*(1-eta);dxphi(2)=-1/(x2-x1)*eta;dxphi(3)=1/(x2-x1)*eta;dxphi(4)=1/(x2-x1)*(1-eta);dyphi(1)=-1/(y3-y1)*(1-ksi);dyphi(2)=1/(y3-y1)*(1-ksi);

11、dyphi(3)=1/(y3-y1)*ksi;dyphi(4)=-1/(y3-y1)*ksi;3. 三角形元的主函数clearxnum=65;ynum=65;pi=3.141592653;point, xn,element,en=mesh2(0,0,1,1,xnum,ynum);A=zeros(xn,xn);b=zeros(xn,1);for ei=1:engx, w=gauss(point , element, ei);phi, dxphi, dyphi=shape(point, element, ei, gx); for i=1:3 ii=element(ei,i);b(ii)=b(ii)

12、+w*8*pi*pi*cos(2*pi*gx(1)*sin(2*pi*gx(2)*p hi(i);for j=1:3jj=element(ei,j);A(ii,jj)=A(ii,jj)+w*(dxphi(i)*dxphi(j)+dyphi(i)*dyphi(j);endendendfor i=1:xnif point(i,2)=0|point(i,2)=1b(i)=0;A(i,:)=zeros(1,xn);A(i, i)=1;endendfor i=1:xnif point(i,1)=0|point(i,1)=1 b(i)=sin(2*pi*point(i,2); A(i,:)=zeros(1

13、,xn);A(i, i)=1;endend x=Ab;%x=inv(A)*b exact=cos(2*pi*point(:,1).*sin(2*pi*point(:,2);%x=inv(A)*berror=0;k=1;for i=1:xnumfor j=1:ynumz(i,j)=x(k);zz(i,j)=exact(k);error=erro 叶(exact(k)-x(k)八2*1/(x nu m-1); k=k+1;end endfprintf(%f,error);deltax=1.0/(xnum-1);deltay=1.0/(ynum-1);mesh(0:deltax:1.0,0:delt

14、ay:1,z)figure;mesh(0:deltax:1.0,0:deltay:1,zz)4. 四边形元的主函数clearxnum=65;ynum=65;pi=3.141592653;point, xn,element,en=mesh1(0,0,1,1,xnum,ynum);pi=3.141592653;A=zeros(xn,xn);b=zeros(xn,1);for ei=1:engx, w=gauss1(point , element, ei);phi,dxphi, dyphi=shape1(point, element, ei,gx);for i=1:4ii=element(ei,i)

15、;b(ii)=b(ii)+w*8*pi*pi*cos(2*pi*gx(1)*sin(2*pi*gx(2)*p hi(i);for j=1:4jj=element(ei,j);A(ii,jj)=A(ii,jj)+w*(dxphi(i)*dxphi(j)+dyphi(i)*dyphi(j);endendendfor i=1:xnif point(i,2)=0|point(i,2)=1b(i)=0;A(i,:)=zeros(1,xn);A(i, i)=1;endendfor i=1:xnif point(i,1)=0|point(i,1)=1b(i)=sin(2*pi*point(i,2);A(i,

16、:)=zeros(1,xn);A(i, i)=1;endendx=Ab;exact=cos(2*pi*point(:,1).*sin(2*pi*point(:,2);%x=inv(A)*b error=0;k=1;for i=1:xnumfor j=1:ynumz(i,j)=x(k);zz(i,j)=exact(k);error=erro 叶(exact(k)-x(k)八2*1/(x nu m-1); k=k+1;endendfprintf(%f,error);deltax=1.0/(xnum-1);deltay=1.0/(ynum-1);mesh(0:deltax:1.0,0:deltay:

17、1,z) figure;上机实验（三）26一、实验内容选用合适的数值方法求解方程:U(44 x2 4 y2) 1 u, (x,y) (0,1) (0,1)u(0, y, t)U(1,y,t)Uy(X,0,t) Uy(X,1,t)0;u(x,y,0)2 2sin( x )cos( y ).求t o.igi.o时，点理备浮环厲半八已色八仝旦）、64 6464 6464 64128 128128 128嵋爲、（227曇、（算256）、喘炭处的数值解。二、实验原理由于计算机只能存储有限个数据和作有限次运算，所以任何 一种用计算机解题的方法，都必须把连续问题（微分方程的边值 问题、初值问题）离散化，最终

18、化成有限形式的线性代数方程组。其求解步骤如下：首先对求解区域进行网格剖分，用有限个 网格节点代替连续区域；其次利用微商代替微分的思想，将微分 算子离散化，构造逼近微分方程定解问题的差分格式，从而把微 分方程的定解问题转化为线性代数方程组问题进行求解。其求解的微分格式如下：n 1 nUj,kUj,k4 4* * Xj * XjnnnUj 1,k25* uj 1,k(厂4* *yykhnnnUj,k 12 uj ,k uj ,k 1)h2)其中u；k表示第n层横坐标为j ,纵坐标为k的点处的函数值，进 行稳定性分析发现，需将区域按 256*256来剖分。二、算法流程对求解区域作网格剖分弋A构造逼近

19、偏微分 方程的差分格式四、计算结果及分析t取不同值题目中要求的各点的值如下表:t=0.1t=0.5t=133 31（阳64）0.45390.21610.0857（色卩）0.22690.18190.0825（47, 49）64 64-0.2317-0.1200-0.04947137、（128 128）0.68590.35900.1498（43 67）128 1280.21130.12680.0530严89）1281280.0112-0.0189-0.0114(127 129)J/256 2560.40360.19320.0766(6391 )256 2560.21620.16550.07413133(,) 256 2560.11620.10630.0492t=0.1s时的图像:o ot=0.5时的图像:0.6 s-f,0 0t=1时的图像:0.2分析结果可以发现： 随着时间 t 的增加，所得结果的绝对值减小。五、总结和讨论本题中所采用的网格剖分达到了稳定性的要求， 算法的时间 和空间复杂度都较大， 以致计算过程比较复杂。 我们需要找到一 种方法，既能达到稳定性的要求，又能提高计算效率，这是我们 今后努力的方向。附录：程序源代码 function wentisanchafen(