《Possion方程的差分法求解.docx》由会员分享，可在线阅读，更多相关《Possion方程的差分法求解.docx（10页珍藏版）》请在金锄头文库上搜索。
1、Possion方程的差分法求解Possion方程的差分法求解 编辑整理：尊敬的读者朋友们：这里是精品文档编辑中心，本文档内容是由我和我的同事精心编辑整理后发布的，发布之前我们对文中内容进行仔细校对，但是难免会有疏漏的地方，但是任然希望（Possion方程的差分法求解）的内容能够给您的工作和学习带来便利。同时也真诚的希望收到您的建议和反馈，这将是我们进步的源泉，前进的动力。本文可编辑可修改，如果觉得对您有帮助请收藏以便随时查阅，最后祝您生活愉快 业绩进步，以下为Possion方程的差分法求解的全部内容。Possion 方程的差分方法课程名称: 微分方程数值解 学生姓名： 张弘 一、问题描述二、问
2、题分析I偏微分方程的数值解法主要有有限差分法和Galerkin有限元法。用差分法和有限元法将连续问题离散化的步骤是:1、对求解区域做网格剖分用有限个网格节点代替连续区域.2、微分算子离散化。3、把微分方程的定解问题化为线性代数方程组的求解问题。差分法和有限元方法的主要区别是离散化的第二步,差分法从定解问题的微分或积分形式出发,用数值微商或数值积分公式到处相应的线性代数方程组，有限元法从定解问题的变分形式出发，用Ritz-Galerkin法导出相应的线性代数方程组。差分法的基本问题有：（1） 对求解区域做网格剖分一维情形为把区间分为等距或不等距的区间，二维情形是把区域分割成均匀或不均匀的矩形，其
3、边与坐标轴平行，也可分割成小三角形或凸四边形。（2） 构造逼近微分方程定解问题的差分格式有两种构造差分格式的方法：直接差分化法和有限体积法。（3） 差分解的存在唯一性,收敛性和稳定性研究（4） 差分方程的解法：逐次超松弛法，交替方向法，共轭梯度法.II（1）由题目可知，本题需要考虑矩形网的五点差分格式。题目给出的为第一边值条件，将十字形图形的中心放于坐标原点处,如下图所示：由图形可知，所给区域可以看出是由8个梯形G1通过旋转、翻转拼接而成。因此为了方便计算、减少计算量，只针对G1进行网格剖分,用5点差分格式进行求解。但是由于G1是直角梯形,进行网格剖分时会出现x轴方向网格点个数不同的现象，不利
4、于有差分形成的系数矩阵的生成，所以将三角形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（ny2my)/2对于矩形区域S1+G1而言,网格内点数为nxny，所以在后面的程序中要比在梯形G1中多计算了（ny2my）/2 个点的函数值,但对程序效率的影响并不是很大,可以接受。以下具体说明只需在G1上求解poisson方程的原因所求
5、方程为：设直线l是经过原点o的任意一条直线,其方程为：y=kx设p(x,y)是区域内任一点,则其关于l对称的点为p(s，t)S=（k1)2+|2y）/(k2+1） t=(2kx+(k21）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)。由和p点的任意性知，对于函数u图像上的任意一点p，其关于任意一条经过原点直线l的对称点p都在u的图像上，即u（+）=u（)，即u关于
6、原点是旋转对称的。当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(i1,j）+4u(i,j)-u(i+1,j）u（i，j1）)=fij(1）对于矩形区域S1+G1,nx=(xb-xa）/h-1ny=(yb
7、ya）/h1按从左到右，从下到上的次序排列网点（ij）：j=1,1inx;j=2,1inx;,;j=ny,1inx，定义向量Uh=u11，u21,unx1；;u1,nx-1,u2，nx-1，,uny1,nx-1T利用边界条件可以将(1）式写成如下形式：其中A=为nx*ny阶矩阵，I为nx阶单位矩阵,B为nx阶单位矩阵。B=右端向量g的元素，依赖于边值a和右端项f，显然A是对称正定矩阵，也是稀疏矩阵因此可以采用逐次超松弛法,共轭梯度法和交替方向法莱求解,但为了方便程序设计采用了matlab的运算符来求解u。针对本题的Poisson方程，对S1+G1形成的矩形网格的五点差分的具体分析.对S1+G1
8、形成的矩形五点差分的具体分析。红色线条围成的区域为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，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：2nx）=2*I； 和A(1,2)
9、=2；保障了上面差分方程的系数.（3）对于网格上最下面一行上除了原点o的所有正则网格内点，由对称性得u（1,i)的邻点中u（1，i1)=u（1，i+1)所以其差分为：4u（i，j)u（i-1，j）-u（i+1，j）2*u(i，j+1)=h2f(i,j）对于右下角的非正则内点u(1,nx)其差分为：4u(i,j）u（i-1，j)2u(i，j+1）=h2*f（i，j)程序中的相关语句为:A（1：nx，nx+1：2nx）=2*I; A(nx+1:2*nx,nx+1：2nx）=B；(4）对于G1中的所有正则内点，其差分为：4u（i，j)u（i1，j）-u（i+1，j)-u（i，j-1)u(i,j+1）
10、=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(inx+1：（i+1)*nx,（i1）*nx+1：i*nx)=I;(5）对于网格最后一列的所有非正则内点u（：,nx），其差分为:4u(nx，j)u（nx，j-1）-u（nx，j+1）-u(nx1,j）=h2f（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，n
11、y+1)；%三角形S1和S2边界上的网点，它们是S1的边界点，但是整个求解区域的内点。三、程序设计及分析function poisson5spot(h)if narginj%因此令A(i,j)=1 A(j，i)=1所以在本程序中多计算了（ny2my)/2 个点的函数值%但对程序的影响并不是很大for i=2:ny A(i1)nx+1:(i-1)nx+i-1,：）=0; for j=1：i1 A（(i-1）*nx+j,(i1）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方程在区域上的图形=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
