偏微分方程报告

成绩成绩2009级数学与应用数学和信息与计算科学专业偏微分方程数值解上机实验实验题目 利用有限元方法和有限差分方法求解偏微分方完成日期 2012年12月17日学生姓名 张灵刚所在班级 1102090任课教师 王晓东西北工业大学理学院应用数学系PAGEPAGE24目录一．实验目的… （2）二．实验要求… （2）三．实验题目… （3）四．实验二… （4）1.实验内容… （4）2.实验原理… （4）.3算法流程… （5）4结果分析… （5）5总结讨论… （6）6源程序… （6）五.实验三1.实验内容… （17）2.实验原理… （17）.3算法流程… （18）4结果分析… 5总结讨论… （21）6源程序… （21）偏微分方程数值解上机实验报告实验地点：数学系机房实验时间13—155、6实验分数30%一、实验目的及意义掌握有限元方法和有限差分方法的程序实现；学会选择合适的有限差分格式求解一维非线性对流占优的非定常对流扩散问限元方法或有限差分方法实现二维初边值抛物型方程的大规模的理解。二、实验要求要求给出格式的推导过程、算法流程、实现程序、选取的网格参数、以表格或图形的方式给出计算结果、对计算结果进行分析、最后对实验进行总结和讨论。问题2：用三角线性元和四边形线性元的有限元方法求解方程：u2cos(2x)sin(2y), (x,y)(0,1)(0,1)u(x,0)u(x,1)0; u(0,y)u(1,y)sin(2y).取h=1/4,1/8,1/16,1/32,1/64,收敛率.问题3：选用合适的数值方法求解方程：(4x2y2)1u, (x,y)(0,1)(0,1)tu(0,y,t)u(1,y,t)u(x,0,t)u(x,1,t)0;y yu(x,y,0)sin(x2)cos(y2).求t0.1,0.5,1.033

31(151747

49)(71

37)(43

67)、64 64 64 64 64 64 128128 128128109 89( ,

、(127,129)

、(63,91

、(31,33)

处的数值解。128128 256 256 256 256 256 256上机实验（二）一、实验内容用三角线性元和四边形线性元的有限元方法求解方程：u2cos(2x)sin(2y), (x,y)(0,1)(0,1)u(x,0)u(x,1)0; u(0,y)u(1,y)sin(2y).取h=1/4,1/8,1/16,1/32,1/64,收敛率二、实验原理Ritz-Galerkin对求解域进行网格剖分三、算法流程对求解域进行网格剖分构造求解问题的基函数构造求解问题的基函数空间步空间步1418116132164形成有限元方程求解有限元方程求解数值收敛率给出实验结果分析三角元实际误差四边形元实际误差0.1298160.0273080.0166760.0036590.0020810.0004570.0002600.0000570.0000320.000007结果分析通过给出的上述结果可以发现三种方法相比四边形元的实际求解误差最小，这是由于四边形型元中的基函数中含有x y项，使得四边形元的误差远比三角形元小，且当空间步长不断加细时向前差分的误差逐渐和四边形元接近而对于收敛率的求解，如果舍去误差，三种方法的收敛率都约为3.五、总结和讨论解域作网络分割，构造基函数（或单元形状函数GalerkinRitzGalerkin附录：程序源代码1.三角形元的型函数function[phi,dxphi,dyphi]=shape(point,element,ei,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);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)=(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-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);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);forei=1:en[gx,w]=gauss(point,element,ei);[phi,dxphi,dyphi]=shape(point,element,ei,gx);fori=1:3ii=element(ei,i);b(ii)=b(ii)+w*8*pi*pi*cos(2*pi*gx(1))*sin(2*pi*gx(2))*phi(i);forj=1:3jj=element(ei,j);A(ii,jj)=A(ii,jj)+w*(dxphi(i)*dxphi(j)+dyphi(i)*dyphi(j));endendendfori=1:xnifpoint(i,2)==0||point(i,2)==1b(i)=0;A(i,:)=zeros(1,xn);A(i,i)=1;endendfori=1:xnifpoint(i,1)==0||point(i,1)==1b(i)=sin(2*pi*point(i,2));A(i,:)=zeros(1,xn);A(i,i)=1;endendx=A\b;%x=inv(A)*bexact=cos(2*pi*point(:,1)).*sin(2*pi*point(:,2));%x=inv(A)*berror=0;k=1;fori=1:xnumforj=1:ynumz(i,j)=x(k);zz(i,j)=exact(k);error=error+(exact(k)-x(k))^2*1/(xnum-1);k=k+1;endendfprintf('%f',error);deltax=1.0/(xnum-1);deltay=1.0/(ynum-1);mesh([0:deltax:1.0],[0:deltay: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);forei=1:en[gx,w]=gauss1(point,element,ei);[phi,dxphi,dyphi]=shape1(point,element,ei,gx);fori=1:4ii=element(ei,i);b(ii)=b(ii)+w*8*pi*pi*cos(2*pi*gx(1))*sin(2*pi*gx(2))*phi(i);forj=1:4jj=element(ei,j);A(ii,jj)=A(ii,jj)+w*(dxphi(i)*dxphi(j)+dyphi(i)*dyphi(j));endendendfori=1:xnifpoint(i,2)==0||point(i,2)==1b(i)=0;A(i,:)=zeros(1,xn);A(i,i)=1;endendfori=1:xnifpoint(i,1)==0||point(i,1)==1b(i)=sin(2*pi*point(i,2));A(i,:)=zeros(1,xn);A(i,i)=1;endendx=A\b;exact=cos(2*pi*point(:,1)).*sin(2*pi*point(:,2));%x=inv(A)*berror=0;k=1;fori=1:xnumforj=1:ynumz(i,j)=x(k);zz(i,j)=exact(k);error=error+(exact(k)-x(k))^2*1/(xnum-1);k=k+1;endendfprintf('%f',error);deltax=1.0/(xnum-1);deltay=1.0/(ynum-1);mesh([0:deltax:1.0],[0:deltay:1],z)figure;上机实验（三）一、实验内容选用合适的数值方法求解方程：(4x2y2)1u, (x,y)(0,1)(0,1)tu(0,y,t)u(1,y,t)u(x,0,t)u(x,1,t)0;y yu(x,y,0)sin(x2)cos(y2).求t0.1,0.5,1.033

31(151747

49)(71

37)(43

67)、64 64 64 64 64 64 128128 128128109 89( ,

、(127,129)

、(63,91

、(31,33)

处的数值解。128128 256 256 256 256 256 256二、实验原理（微分方程的边值分方程的定解问题转化为线性代数方程组问题进行求解。其求解的微分格式如下：un1un

1 un

2*un

un

un 2*un

unj,k

j,k

j

j,k

j1,k

j,k

j,k

j,k1) 44**x*xj j

4**y*y h2 h2k k其中unj,k

表示第n层横坐标为j，纵坐标为k的点处的函数值，进行稳定性分析发现，需将区域按256*256来剖分。对求解区域作网格剖分构造逼近偏微分方程的差分格式三、算法流程对求解区域作网格剖分构造逼近偏微分方程的差分格式稳定性分析稳定性分析确定剖分格数显式差分格式进行求解四、计算结果及分析(33(33,31)64 64(15,17)64 64(47,49)64 64(71,37)128 128(43,67)128 128(10989128 128,)t=0.1t=0.5t=10.45390.21610.08570.22690.18190.0825-0.2317-0.1200-0.04940.68590.35900.14980.21130.12680.05300.0112-0.0189-0.01140.40360.19320.07660.21620.16550.07410.11620.10630.0492(127,(127,129)256 256(63,