Poisson方程九点差分格式_米瑞琪.doc

数值实验报告I实验名称Poisson方程九点差分格式实验时间2016年 4 月 15 日姓名米瑞琪班级信息1303学号1309010304成绩一、实验目的，内容 1、理解Poisson方程九点差分格式的构造原理； 2、理解因为网格点的不同排序方式造成的系数矩阵格式的差异； 3、学会利用matlab的spdiags(),kron()函数生成系数矩阵；二、算法描述针对一个Poisson方程问题：u=-f(x,y) 在Poisson方程五点差分格式的基础上，采用Taylor展开分析五点差分算子的截断误差，可以得到：huxi,yj=uxi,yj+112h124uxi,yjx4+h224uxi,yjy4+Oh4=uxi,yj+112h122x2+h222y22uxi,yjx2+2uxi,yjy2-h12+h22124uxi,yjx2y2+Oh4 (1)为了提高算子截断误差的精度，在(1)式中配凑出了差分算子的形式，将原Poisson方程代入（1）式有：huxi,yj=uxi,yj+112h122x2+h222y22u(xi,yj)x2+2u(xi,yj)y2-h12+h22124uxi,yjx2y2+Oh4 (2)考虑4uxi,yjx2y2，有：4uxi,yjx2y2=1h22uxxxi,yj+1-2uxxxi,yj+uxxxi,yj-1+Oh22=1h22uxi+1,yj+1-2uxi,yj+1+uxi-1,yj+1h12-h12124uxi,yj+1x4-2uxi+1,yj-2uxi,yj+uxi-1,yjh12+2h12124uxi,yjx4+uxi+1,yj-1-2uxi,yj-1+uxi-1,yj-1h12-h12124uxi,yj-1x4=1h22Lhuxi,yj+1-2Lhuxi,yj+Lhuxi,yj-1-h12h2125uxi,yj+1x4y+h12h2125uxi,yj-1x4y=1h22Lhuxi,yj+1-2Lhuxi,yj+Lhuxi,yj-1-ch12h22126uxi,yjx4y2(3)将（3）代回（2）可得huxi,yj+h12+h2212h12h224uxi,yj-2uxi,yj+1+uxi,yj-1+uxi+1,yj+uxi-1,yj+uxi+1,yj+1+uxi+1,yj-1+uxi-1,yj+1+uxi-1,yj-1=-fxi,yj-112h122fxi,yjx2+h222fxi,yjy2 得到Poisson方程的九点差分格式：huij+h12+h2212h12h224uij-2ui,j+1+ui,j-1+ui+1,j+ui-1,j+ui+1,j+1+ui+1,j-1+ui-1,j+1+ui-1,j-1=-fij-112h122fijx2+h222fijy2 (4)在计算机上实现（4）式，需要在五点差分格式huij=-fij的基础上在等式两端分别增加一部分，将等式左侧新增的部分写成紧凑格式，有：4-2-21-24-21-21-241-2-214-2-211-21-24-21-211-2-241-2-21-21-214-21-24-21-2-2-2-24u11u12u1,N2-1u21u22u2,N2-1uN1-1,1uN1-1,2uN1-1,N2-1对于该矩阵，可以看成是两个矩阵的组合：C=4-2-24-2-244-2-24-2-244-2-24-2-2-2-24以及D=-211-211-2-21-211-211-211-21-2-21-21-2111-2则生成这两个矩阵可以采用Kroncker生成，方法类似于五点差分格式。对于右端添加的关于f(x,y)的二阶导数，可以采用中心差分格式进行近似代替，即：112h122fijx2+h222fijy2=112fi+1,j-2fij+fi-1,j+fi,j+1-2fij+fi,j-1=112-4fij+fi+1,j+fi-1,j+fi,j+1+fi,j-1写成相应的紧凑格式有：112h122fijx2+h222fijy2=112-4111-4111-411-41111-41111-41111-411-4111111-4f11f12f1,N2-1f21f22f2,N2-1fN1-1,1fN1-1,2fN1-1,N2-1+112f01f02f0,N2-1000fN1,1fN1,2fN1,N2-1+112f100f1,N2f200f2,N2fN1-1,00fN1-1,N2该式中的矩阵又可以分解为两个矩阵的和：E=-411-411-4-411-411-4-411-41111-4F=11111111111111三程序代码 根据上述分析，可以写出程序代码为：clcclear%网格剖分xa=0;xb=1;N1=64;h1=(xb-xa)/N1;x=xa+1:(N1-1)*h1;ya=0;yb=1;N2=64;h2=(yb-ya)/N2;y=ya+1:(N2-1)*h2;%生成网格点与边界点处f的函数值，其中f_b1是中间全部为0的向量，f_b2是每一个分块中间为0的向量f=(x,y)-(2*pi2)*exp(pi*(x+y)*(sin(pi*x)*cos(pi*y)+cos(pi*x)*sin(pi*y);f_b1=zeros(N1-1)*(N2-1),1);f_b2=zeros(N1-1)*(N2-1),1);for i=1:N1-1 f_b2(i-1)*(N2-1)+1)=f(x(i),ya); f_b2(i*(N2-1)=f(x(i),yb); for j=1:N2-1 f_in(i-1)*(N2-1)+j)=f(x(i),y(j); if i=1 f_b1(j)=f(xa,y(j); end if i=N1-1 f_b1(i-1)*(N2-1)+j)=f(xb,y(j); end endendf_in=f_in;%生成迭代矩阵Ae2=ones(N1-1,1);K1=spdiags(e2,-2*e2,e2,-1,0,1,N1-1,N1-1);e3=ones(N2-1,1);I1=spdiags(e3,0,N2-1,N2-1);A=kron(K1,I1);A=A/h12;%生成迭代矩阵B和CI2=spdiags(e2,0,N1-1,N1-1);K2=spdiags(e3,-2*e3,e3,-1,0,1,N2-1,N2-1);B=kron(I2,K2);C=-2*(h12+h22)/(12*h12*h22)*B;B=B/h22;%生成迭代矩阵DS1=spdiags(e2 e2,-1 1,N1-1,N1-1);S2=spdiags(e3,0,N2-1,N2-1);D=(h12+h22)/(12*h12*h22)*kron(S1,K2);%生成f的中心差分格式对应矩阵E和FE=kron(S1,S2);K3=spdiags(e3,-4*e3,e3,-1 0 1,N2-1,N2-1);F=kron(I2,K3);f_vec=-f_in-(1/12)*(E+F)*f_in+f_b1+f_b2);%生成网格点X,Y=meshgrid(x,y);%计算网格函数值u_d=(A+B+C+D)f_vec;u0=zeros(length(f_vec),1);%计算误差u_real=(x,y)exp(pi*(x+y)*sin(pi*x).*sin(pi*y);for i=1:N1-1 u_m(i-1)*(N2-1)+1:i*(N2-1)=u_real(x(i),y);endu_v=u_m;err_d=max(abs(u_d-u_v);sol=reshape(u_d,N2-1,N1-1);mesh(X,Y,sol)四. 数值结果针对课本P93给出的问题，分别采用步长164,1128，将计算出的误差列表如下：步长五点差分格式误差九点差分格式误差1640.02861330711346713.59173530739554e-0611280.007156344071610482.24435417806035e-07可见采用九点差分格式可以进一步缩小误差，达到更高阶的精度。五. 计