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

数值实验报告I实验名称Poisson方程九点差分格式实验时间2016年 4 月 15 日姓名米瑞琪班级信息1303学号1309010304成绩一、实验目的，内容 1、理解Poisson方程九点差分格式的构造原理； 2、理解因为网格点的不同排序方式造成的系数矩阵格式的差异； 3、学会利用matlab的spdiags(),kron()函数生成系数矩阵；二、算法描述针对一个Poisson方程问题： 在Poisson方程五点差分格式的基础上，采用Taylor展开分析五点差分算子的截断误差，可以得到：为了提高算子截断误差的精度，在(1)式中配凑出了差分算子的形式，将原Poisson方程代入（1）式有：考虑，有：将（3）代回（2）可得 得到Poisson方程的九点差分格式：在计算机上实现（4）式，需要在五点差分格式的基础上在等式两端分别增加一部分，将等式左侧新增的部分写成紧凑格式，有：对于该矩阵，可以看成是两个矩阵的组合：以及则生成这两个矩阵可以采用Kroncker生成，方法类似于五点差分格式。对于右端添加的关于f(x,y)的二阶导数，可以采用中心差分格式进行近似代替，即：写成相应的紧凑格式有：该式中的矩阵又可以分解为两个矩阵的和：三．程序代码 根据上述分析，可以写出程序代码为：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*pi^2)*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/h1^2;%生成迭代矩阵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*(h1^2+h2^2)/(12*h1^2*h2^2)*B;B=B/h2^2;%生成迭代矩阵DS1=spdiags([e2 e2],[-1 1],N1-1,N1-1);S2=spdiags(e3,0,N2-1,N2-1);D=(h1^2+h2^2)/(12*h1^2*h2^2)*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给出的问题，分别采用步长，将计算出的误差列表如下：步长五点差分格式误差九点差分格式误差0.02861330711346713.59173530739554e-060.007156344071610482.24435417806035e-07可见采用九点差分格式可以进一步缩小误差，达到更高阶的精度。