Possion方程的差分法求解(共10页)

1、Possion 方程的差分方法课程名称： 微分方程数值解 学生姓名： 张弘 一、问题描述二、问题分析I偏微分方程的数值解法主要有有限差分法和Galerkin有限元法。用差分法和有限元法将连续问题离散化的步骤是：1、对求解区域做网格剖分用有限个网格节点代替连续区域。2、微分算子离散化。3、把微分方程的定解问题化为线性代数方程组的求解问题。差分法和有限元方法的主要区别是离散化的第二步，差分法从定解问题的微分或积分形式出发，用数值微商或数值积分公式到处相应的线性代数方程组，有限元法从定解问题的变分形式出发，用Ritz-Galerkin法导出相应的线性代数方程组。差分法的基本问题有：（1） 对求解区域

2、做网格剖分一维情形为把区间分为等距或不等距的区间，二维情形是把区域分割成均匀或不均匀的矩形，其边与坐标轴平行，也可分割成小三角形或凸四边形。（2） 构造逼近微分方程定解问题的差分格式有两种构造差分格式的方法：直接差分化法和有限体积法。（3） 差分解的存在唯一性，收敛性和稳定性研究（4） 差分方程的解法：逐次超松弛法，交替方向法，共轭梯度法。II（1）由题目可知，本题需要考虑矩形网的五点差分格式。题目给出的为第一边值条件，将十字形图形的中心放于坐标原点处，如下图所示：OS1G1S2由图形可知，所给区域可以看出是由8个梯形G1通过旋转、翻转拼接而成。因此为了方便计算、减少计算量，只针对G1进行网格

3、剖分，用5点差分格式进行求解。但是由于G1是直角梯形，进行网格剖分时会出现x轴方向网格点个数不同的现象，不利于有差分形成的系数矩阵的生成，所以将三角形S1和梯形G1合在一起形成一个矩形，其区域为0,3/20,1/2。如果采用等距差分，并且有hx=hy=h，设步长为h，x=0:h:3/2;y=0:h:1/2;nx=length(x)-1;%为所求区域中x轴上内点的个数ny=length(y)-1; %为所求区域中y轴上内点的个数对于原来的梯形G1来说，有网格内点nx*ny-(ny2-my)/2对于矩形区域S1+G1而言，网格内点数为nx*ny，所以在后面的程序中要比在梯形G1中多计算了(ny2-

4、my)/2 个点的函数值，但对程序效率的影响并不是很大，可以接受。以下具体说明只需在G1上求解poisson方程的原因所求方程为：设直线l是经过原点o的任意一条直线，其方程为：y=kx设p(x,y)是区域内任一点，则其关于l对称的点为p(s,t)S=(k-1)2+|2y)/(k2+1) t=(2kx+(k2-1)y)/(k2+1)则同理可得：将u(s,t)代替u(x,y)得：Uxx+uyy=uss+utt=1其依然满足上述方程。令=arctan(k)点p的横坐标x=rcos() y=rsin()则p关于直线l的对称点为p(rcos(2-),rsin(2-)由上述证明可知u(p)=u(p)。由和

5、p点的任意性知，对于函数u图像上的任意一点p，其关于任意一条经过原点直线l的对称点p都在u的图像上，即u(+)=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方程的差分格式和化为如下形式：对于正则内点其差分如下：-uij=1/h2*(-u(i,j+1)-u(i-1

6、,j)+4u(i,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,1inx；j=2,1inx;,;j=ny,1inx，定义向量Uh=u11,u21,unx-1;u1,nx-1,u2,nx-1,uny-1,nx-1T利用边界条件可以将（1）式写成如下形式:其中A=为nx*ny阶矩阵，I为nx阶单位矩阵，B为nx阶单位矩阵。B=右端向量g的元素，依赖于边值a和右端项f，显然A是对称正定矩阵，也是稀疏矩阵因此可以采用逐次超松弛法，共轭梯度法和交替方向法莱求解，但为

7、了方便程序设计采用了matlab的运算符来求解u。针对本题的Poisson方程，对S1+G1形成的矩形网格的五点差分的具体分析。对S1+G1形成的矩形五点差分的具体分析。ny+1ny3211 2 3 4 i nx-1 nx nx+1G1S2S11/2O3/21/2红色线条围成的区域为G1，L为红色斜边，S1为L左侧的三角形，S2为L右侧的三角形。由对称性知，S1和S2中关于L相互对称的网格点其上U的函数值是相同的。按从左到右，从下到上的次序排列网点(ij)。（1）对于三角形S1中的网点有u(i,j)， nyij1，有u(i,j)-u(j,i)=0 所以S1中网点的差分就为：u(i,j)-u(j

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

9、)=h2*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)=h2*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,j

10、)-u(nx,j-1)-u(nx,j+1)-u(nx-1,j)=h2*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的边界点，但是整个求解区域的内点。三、程序设计及分析function poisson5spot(h)if narginj%因此令A(i,j)=1 A(j,i)=-1%所以在本程序中多计算了(ny2-my)/2 个点的函数值%但对程序的影响并不是很大for i=2:

11、ny A(i-1)*nx+1:(i-1)*nx+i-1,:)=0; for j=1:i-1 A(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);for i=1:ny for j=1:nx u(i,j)=x(i-1)*nx+j); endend%根据对称性，给网格最上面一行的点赋值u(ny+1,1:ny)=u(1:ny,ny+1); %=作Poisson方

12、程在区域上的图形=x,y=meshgrid(0:h:3/2,0:h:1/2);hold onsurf(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.Lhuij0，对任意(xi,yj)Gh，uij不可能在内点取得负的极小，与极值定理相符合。3、翻转后图形如下：4 网格间距h=0.01时得到的图形： 五、实验分析本次实验将求解区域先利用对称性将其缩小为原区域的1/