Poisson方程九点差分格式-米瑞琪

数值实验报告实验名称Poisson方程九点差分格式实验时间实验名称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*piA2)*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);fori=1:N1-1f_b2((i-1)*(N2-1)+1)=f(x(i),ya);f_b2(i*(N2-1))=f(x(i),yb);forj=1:N2-1f_in((i-1)*(N2-1)+j)=f(x(i),y(j));ifi==1f_b1(j)=f(xa,y(j));endifi==N1-1

f_b1((i-1)*(N2-1)+j尸f(xb,y(j));endendendf_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/h1A2;%生成迭代矩阵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*(h1A2+h2A2)/(12*h1A2*h2A2)*B;B=B/h2A2;%生成迭代矩阵DS1=spdiags([e2e2],[-11],N1-1,N1-1);S2=spdiags(e3,0,N2-1,N2-1);D=(h1A2+h2A2)/(12*h1A2*h2A2)*kron(S1,K2);%生成f的中心差分格式对应矩阵E和FE=kron(S1,S2);K3=spdiags([e3,-4*e3,e3],[-101],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);fori=1:N1-1u_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给出的问题，分别采用步长，将计算出的误差列表如下：Jk五点差分格式误差九点差分格式误差—0.02861330711346713.59173530739554e-06—0.007156344071610482.24435417806035e-07可见米用九点差分格式可以进一步缩小误差，达到更高阶的精度。

五.计算中出现的问题，解决方法及体会在生成九点差分格式的时候，等号右端涉及到了对f的二阶偏导，我最初利用符号函数定义了f,随后求出具二