Possion方程的差分法求解

1、Possion方程的差分方法课程名称：微分方程数值解学生姓名：张弘一、问题描述Posskm方程的差分方法设C星如下图所示的十字形区域.由b作机塔的单柱正方形加心1A用力.点差分恬式求解卜由的Pzim右程第边*宜问期:并以表格或阳彤的方式表于你得到的数值解.需望你的和序能适值仟何十氏.提办二出区域和施分方程的对聊性，只黑在如图阴影所小的八分之一区域5上求解期呼（你能证曳这,点吗？L_二、问题分析I.偏微分方程的数值解法主要有有限差分法和Galerkin有限元法。用差分法和有限元法将连续问题离散化的步骤是：1、对求解区域做网格剖分用有限个网格节点代替连续区域。2、微分算子离散化。3、把微分方程的定

2、解问题化为线性代数方程组的求解问题。差分法和有限元方法的主要区别是离散化的第二步，差分法从定解问题的微分或积分形式出发，用数值微商或数值积分公式到处相应的线性代数方程组，有限元法从定解问题的变分形式出发，用Ritz-Galerkin法导出相应的线性代数方程组。差分法的基本问题有：(1) 对求解区域做网格剖分一维情形为把区间分为等距或不等距的区间，二维情形是把区域分割成均匀或不均匀的矩形，其边与坐标轴平行，也可分割成小三角形或凸四边形。(2) 构造逼近微分方程定解问题的差分格式有两种构造差分格式的方法：直接差分化法和有限体积法。(3) 差分解的存在唯一性，收敛性和稳定性研究(4) 差分方程的解法

3、：逐次超松弛法，交替方向法，共轲梯度法。II(1)由题目可知，本题需要考虑矩形网的五点差分格式。题目给出的为第一边值条件，将十字形图形的中心放于坐标原点处，如下图所示：由图形可知，所给区域可以看出是由8个梯形G1通过旋转、翻转拼接而成。因此为了方便计算、减少计算量，只针对G1进行网格剖分，用5点差分格式进行求解。但是由于G1是直角梯形，进行网格剖分时会出现x轴方向网格点个数不同的现象，不利于有差分形成的系数矩阵的生成，所以将三角形S1和梯形G1合在一起形成一个矩形，其区域为0,3/2x02。如果采用等距差分，并且有hx=hy=h,设步长为h,x=0:h:3/2;y=0:h:1/2;nx=len

4、gth(x)-1;%为所求区域中x轴上内点的个数ny=length(y)-1;%为所求区域中y轴上内点的个数对于原来的梯形G1来说，有网格内点nx*ny-(nyA2-my)/2对于矩形区域S1+G1而言，网格内点数为nx*ny,所以在后面的程序中要比在梯形G1中多计算了(nyA2-my)/2个点的函数值，但对程序效率的影响并不是很大，可以接受。以下具体说明只需在G1上求解poisson方程的原因所求方程为：u0anC设直线l是经过原点o的任意一条直线，其方程为：y=kx设p(x,y)是区域内任一点，则其关于l对称的点为p'(s,t)S=(k-1)A2+|2y)/(H2+1)t=(2kx

5、+(kA2-1)y)/(kA2+1),1k22kuxussxuttxus2ut2k1k11k222k(1k2)2k2uxxuss(/,2)2ust2.2utt(),2)1k(1k)1k同理可得：_222k2o2k(k1)k12uyyuss(1k2)2ust(1k2)2utt(1k2将u(s,t)代替u(x,y)得:Uxx+uyy=uss+utt=1其依然满足上述方程。令0=arctan(k)岚p的横坐标x=rcos(a)y=rsin(a)则p关于直线l的对称点为p'(rcos(2-0c),rsin(2-a0)由上述证明可知u(p)=u(p')。由。和p点的任意性知，对于函数u图

6、像上的任意一点p,其关于任意一条经过原点直线的对称点p'都在u的图像上，即u(a+8尸u(即)u关于原点是旋转对称的。当l为x轴时，有u(x,y)=u(x,-y)l为y轴时,u(x,y)=u(-x,y)坐标轴旋转不改变方程的形式，所以函数在直角坐标系中u关于原点是旋转对称的，又所求十字形区域关于x,y轴是轴对称和关于原点中心对称的，因此可通过直求解区域G1,就可以知道函数在整个十字形区域的图像。(2)对S1+G1形成的矩形进行正方形网格剖分，则求解Poisson方程的差分格式和化为如下形式：对于正则内点其差分如下：-Auj=1/hA2*(-u(i,j+1)-u(i-1,j)+4u(i,

7、j)-u(i+1,j)-u(i,j-1)=fij(1)对于矩形区域S1+G1,nx=(xb-xa)/h-1ny=(yb-ya)/h-1按从左到右，从下到上的次序排列网点(ij):j=1,1<i<n)j=2,1wiwnx;,rny,1,标i乎朝量Uh=ui1,U21，Ux-1；u、x-1,u2,nx-1,ny-1,nx-lT利用边界条件可以将(1)式写成如下形式：1 A一2 AuhghBIIBI其中A=为nx*ny阶矩阵，I为nx阶单位矩阵，B为nx阶单位IBIIB矩阵。41141B=.14114右端向量g的元素，依赖于边值a和右端项f,显然A是对称正定矩阵，也是稀疏矩阵因此可以采用

8、逐次超松弛法，共轲梯度法和交替方向法莱求解，但为了方便程序设计采用了matlab的V运算符来求解u。针对本题的Poisson方程，对S1+G1形成的矩形网格的五点差分的具体分析。对S1+G1形成的矩形五点差分的具体分析。1/2ny+1ny321O1/232红色线条围成的区域为G1,L为红色斜边，S1为L左侧的三角形，S2为L右侧的三角形。由对称性知，S1和S2中关于L相互对称的网格点其上U的函数值是相同的。按从左到右，从下到上的次序排列网点(ij)。(1)对于三角形S1中的网点有u(i,j),ny>i>j6盾u(i,j)-u(j,i)=0所以S1中网点的差分就为：u(i,j)-u(

9、j,i)=0(2)由于原点。点的特殊性，其邻点中u(1,2)=u(1,-1)=u(-1,1)=u(2,1)所以其差分为：u(1,1)-4u(1,2)=hA2*f(i,j)程序中语句：A(1:nx,nx+1:2*nx)=2*I;和A(1,2)=-2;保障了上面差分方程的系数。(3)对于网格上最下面一行上除了原点。的所有正则网格内点，由对称性得u(1,i)的邻点中u(1,i-1)=u(1,i+1)所以其差分为:4u(i,j)-u(i-1,j)-u(i+1,j)-2*u(i,j+1)=hA2*f(i,j)对于右下角的非正则内点u(1,nx)其差分为：4u(i,j)-u(i-1,j)-2*u(i,j+

10、1)=hA2*f(i,j)程序中的相关语句为：A(1:nx,nx+1:2*nx)=2*I;A(nx+1:2*nx,nx+1:2*nx)=B;(4)对于G1中的所有正则内点，其差分为：4u(i,j)-u(i-1,j)-u(i+1,j)-u(i,j-1)-u(i,j+1)=hA2*f(i,j)程序中相关语句：A(i*nx+1:(i+1)*nx,i*nx+1:(i+1)*nx尸B;A(i-1)*nx+1:i*nx,i*nx+1:(i+1)*nx)=I;A(i*nx+1:(i+1)*nx,(i-1)*nx+1:i*nx)=I;(5)对于网格最后一列的所有非正则内点u(:,nx),其差分为：4u(nx,

11、j)-u(nx,j-1)-u(nx,j+1)-u(nx-1,j)=hA2*f(i,j)(6)在求出了所有的内点后，对于剩下的边界点赋值有：U(ny+1,ny+1:nx)=0%最上面一行上的边界点U(1:ny+1,nx+1)=0%最右侧一列的边界点u(ny+1,1:ny)=u(1:ny,ny+1);%三角形S1和S2边界上的网点，它们是S1的边界点，但是整个求解区域的内点。三、程序设计及分析functionpoisson5spot(h)ifnargin<1%默认步长h=h=;endx=0:h:3/2;y=0:h:1/2;nx=length(x)-1;%取区域的内点ny=length(y)-

12、1;%=勾造矢I阵B=B=eye(nx,nx)*4;fori=1:nx-1B(i,i+1)=-1;endfori=2:nxB(i,i-1)=-1;endI=-eye(nx,nx);%=勾造系数矩阵a=A=zeros(nx*ny,nx*ny);A(1:nx,1:nx)=B;%由区域的对称性，正方形网格最下面一行的差分形式%改为4u(i,j)-u(i-1,j)-u(i+1,j)-2*u(i,j+1)%因为u(i,j+1)=u(i,j-1)A(1:nx,nx+1:2*nx)=2*I;A(nx+1:2*nx,nx+1:2*nx)=B;%矩形网格左下角第一个点的差分形式改为：%4u(i,j)-u(i-1

13、,j)-u(i+1,j)-u(i,j+1)-u(i,j-1)%=4u(i,j)-2u(i,j+1)-2u(i+1,j)A(1,2)=-2;%为了方便，本程序往梯形中增加了一个三角形，以方便编程A(nx+1:2*nx,1:nx)=I;fori=2:ny-1A(i*nx+1:(i+1)*nx,i*nx+1:(i+1)*nx)=B;A(i-1)*nx+1:i*nx,i*nx+1:(i+1)*nx)=I;A(i*nx+1:(i+1)*nx,(i-1)*nx+1:i*nx)=I;end%=%bi=h*h*f;右端向量b=h*h*ones(nx*ny,1);%由于对称性，左侧三角形中有u(i,j)-u(j

14、,i)=0i>j%因此令A(i,j)=1A(j,i)=-1%所以在本程序中多计算了(nyA2-my)Z2个点的函数值%但对程序的影响并不是很大fori=2:nyA(i-1)*nx+1:(i-1)*nx+i-1,:)=0;forj=1:i-1A(i-1)*nx+j,(i-1)*nx+j)=1;A(i-1)*nx+j,(j-1)*nx+i)=-1;b(i-1)*nx+j)=0;endendx=Ab;%为了方便采用左除求解网格点上的函数值%x=gmres(A,b,100);%按顺序将x赋值给uu=zeros(ny+1,nx+1);fori=1:nyforj=1:nxu(i,j)=x(i-1)*

15、nx+j);endend%根据对称性，给网格最上面一行的点赋值u(ny+1,1:ny)=u(1:ny,ny+1);%=Poisson方程在区域上的图形=x,y=meshgrid(0:h:3Z2,0:h:1Z2);holdonsurf(x,y,u)%11第一象限的第一块区域，下面的以此类推surf(y,x,u)%12surf(-y,x,u)%21surf(-x,y,u);%22surf(-x,-y,u)%31surf(-y,-x,u)%32surf(y,-x,u)%41surf(x,-y,u);%42四、实验结果1.在区域G1上求解后的图形显示：图形表示了在S1+G1区域上的函数图像，而不是单纯的G1上的函数图像。2通过拼接后图形显示如下：由上图可知，除了边界点外网格点上的函数值都有u(i,j)>0.LhUij>0,对任意(xi,yj)CGh,uij不可能在内点取得负的极小，与极值定理相符合。3、翻转后图形如下：时。Edri'itwJciokQeioo岌刖6u._：jiftt：与、k