收稿日期:2005212210基金项目:辽宁省教育厅科研基金资助项目(05L415)・作者简介:刘大卫(1964-),男,贵州贵阳人,贵州工业大学副教授・第24卷　第2期2006年4月沈阳师范大学学报(自然科学版)Journal of S henyang Norm al U niversity (N atural Science )V ol 124,N o.2Apr.2006文章编号:1673-5862(2006)02-0166-04正方形环域Laplace 方程的简明数值解法刘大卫1,高　明2,3(1.贵州工业大学基础部,贵州贵阳　550003; 2.沈阳师范大学物理科学与技术学院,辽宁沈阳　110034;3.沈阳师范大学实验中心,辽宁沈阳　110034)摘 要:通过正方形环域的Laplace 方程的数值求解过程,详细介绍了使用MA TLAB 求解微分方程的方法・用MA TLAB 的M 文件,生成正方形环域,用函数numgrid 作网格划分,用函数delsq 建立五点差分格式建立并求解拉普拉斯方程第一边值问题・关　键　词:Laplace 方程;差分法;MA TLAB 中图分类号:O 175 文献标识码:A0　引 言Laplace 方程是解决电磁场问题中最常见的方程,在一些具有较复杂边界形状的区域中求出方程的解析解是非常困难的[122]・因此寻求一种有效的、简明的数值解法对于解决实际问题中复杂边界区域中的电磁场分布问题具有非常重要的实际价值・通过一个特殊的方形区域的电场分布问题介绍一种应用MA TLAB 数值求解Laplace 方程的方法・考虑图1所示正方形环域,设区域内满足Laplace 方程Δu =0,内边界处电势u =100,外边界处电势u =0,求区域内的电势分布,易见,这是一个Laplace 方程的第一边值问题・现用差分法求解这个问题,首先把研究区域划分为图2所示的网格,在这个划分中,除去边界点,区域被分为240个网格节点・图1　正方形环域图2　网格的划分差分法求解的基本思想是,在网格节点上用差商代替微商,结合边界条件,把定解问题转化为以未知函数u (x ,y )在节点上的数值为未知量的线性方程组:Ax =b其中,x 为解向量,代表函数u (x ,y )在节点上的数值・A 为系数矩阵,与网格节点的划分和编号方式有关,通常是一个大型的稀疏矩阵・b 为常数向量,由边界条件确定・对上述问题,A 为240×240阶稀疏矩阵,b 为240×1阶稀疏常数向量・下面用MA TLAB 提供的网格划分函数numgrid 和差分格式建立函数delsq 来构造系数矩阵A ・1　编写M 文件,生成正方形环域MA TLAB 提供了一些常见区域如正方形、圆、圆环、L 型域等区域上的网格划分函数numgrid ,但函数的可行域中没有正方形环域・为此,打开网格划分函数numgrid.m ,在其命令行中加入域名为‘N1’的新可行域,根据实际物理问题,编写命令行,生成正方形环域[1]・2　网格的划分和差分格式的建立在MA TLAB 命令窗口键入:µR =‘N1’;定义正方形环域N1・µn =19;图3　函数numgrid 划分的正方形环域定义划分网格数・除去边界点之外,此划分定义了(19-2)×(19-2)-7×7=240个内节点・µg =numgrid(R ,n );在R 上按n 划分区域・函数numgrid 把区域划分为正方形环域网格,生成数字矩阵g ・在g 中,边界节点编号为0・内部节点的编号为1～240,240为内部节点数・节点在域中的位置与节点对应元素在矩阵g 中的位置相同・而内部节点的编号,则表示在方程Ax =b 中,该节点在向量x 中所在的行・现用命令spy 查看矩阵g 的结构图:µspy (g );图4　五点差分格式矩阵A 的结构图结果如图3所示,它直观地显示了所定义的新可行域‘N1’的确将区域划分为正方形环域网格,240个非零元素表明区域被划分为240个格点・然后,用函数delsq 建立五点差分格式:µA =delsq (g );函数delsq 在g 上对240个内节点建立五点差分格式,生成系数矩阵A ,A 中第i 行表示第i 个节点,节点i 在域中的位置可通过数字矩阵g 查询・现用命令spy 查看矩阵A 的结构图:µspy (A );结果如图4所示,240×240阶系数矩阵A 中有1104个非零元素,它表明系数矩阵A 确是一个大型稀疏矩阵・3　常数矩阵b 的建立、稀疏矩阵技术在图2所示的网格中,对内部不邻近边界的节点,其邻近的4个节点都为内部节点,故与之对应方程的常数项为0・对内部邻近边界的节点,其邻近的4个节点中有一个边界点,与之对应方程的常数项应为边界条件值・这里讨论的问题,外侧边界条件均为0,只有内侧邻近边界的28个节点的对应方程的常数项不为0,因而常数矩阵b 是一个240行且其中只有28个非零元素的稀疏列向量,要得到常数矩阵b ,需要MA TLAB 提供的稀疏矩阵技术・首先考查内侧邻近边界的28个节点的位置编号,因它们分别在数字网格划分矩阵g 中的第6行到第14行、第6列到第14列,键入:µg (6:14,6:14)可得到这28个节点的位置编号(用斜体表示)・761第2期 刘大卫等:正方形环域Laplace 方程的简明数值解法7390100110120130140150160740000000161750000000162760000000163770000000164780000000165790000000166801678191101111121131141151168如上侧邻近内边界的7个节点的位置编号分别为:90、100、110、120、130、140、150等等・然后,把这28个节点在常数矩阵b 中对应的行、列(列号为1)、边界条件值(u =100)写成数据文件mydat.dat ,其内容如下:90 1 10011011001201100 (240)10.00该文件指定了即将生成的常数矩阵b 中指定行、列的元素值,最后一行确定常数矩阵b 的阶数・保存mydat.dat 文件,并用load 命令加载到工作间后键入:µb =spconvert (mydat );函数spconvert 把外部数据文件mydat.dat 转化为稀疏矩阵・注意到在mydat.dat 中,最大行号为240,最大列号为1,函数spconvert 将数据文件mydat.dat 转化为240×1的、只在指定的28个位置上有非零元素的稀疏列向量,即把边界条件传递到常数矩阵b 中・4　解的查询对以上建立的大型差分方程Ax =b ,采用G auss 2Seidel 法[325]进行求解,其中的迭代精度为110×10-6,初值设为0,为了与坐标变量x 相区别,记解向量为u ・现查看解的结果・由于区域和边界条件的对称性,我们只需查看如图5所示的35个网格节点上的结果・键入:图5　在35个节点上的解µg (1:7,1:10)查询这35个网格节点的位置编号,结果如下:11835526986961061161936537087971071173754718898108118557289991091197390100110120例如,图5中第5列的5个节点的位置编号是:69、70、71、72、73,通过键入:µu (69:73)可查询解向量u 中这5个节点处的解为:10.6　21.5　33.1　45.5　59.1同此查出其他节点处的解,并标示于图5中・5　解的图示解向量u 为240行的列向量,利用MA TLAB 的数值图形化技术[2],可以将其在坐标网格上直观地显示出来・861沈阳师范大学学报(自然科学版) 第24卷图6　矩阵U 的三维网格图µU =g ;用矩阵g 的结构定义矩阵U ・µU (g >0)=full (u (g (g >0)));在g 的非零元素位置处用u 的值置换g 中元素并显示为满阵型・µmesh (U );colormap (gray/2);以灰度方式显示U 的三维网格图,结果如图6所示・ (注:关于在网格划分函数numgrid.m 中生成新可行域的程序,有兴趣的读者,欢迎与作者联系)参考文献:[1]珀塞尔E M.电磁学[M ].北京:科学出版社,1979.[2]王纪林,向光辉.特殊函数与数学物理方程[M ].上海:上海交通大学出版社,1989.[3]王沫然.MA TLAB 610与科学计算[M ].北京:电子工业出版社,2001.[4]陆金甫,关　治.偏微分方程数值解法[M ].北京:清华大学出版社,2004.[5]周长发.科学与工程数值算法[M ].北京:清华大学出版社,2004.Simple Numerical Computation of Laplace ’s Equationon Square with Smaller Square H oleL IU Da 2wei 1,GA O M i ng2,3(1.Guizhou University of Technology ,Guiyang 550003,China ;2.College of Physics Science and Technology ,Shenyang Normal University ,Shenyang 110034,China ;3.Experimental Center ,Shenyang Normal University ,Shenyang 110034,China )Abstract :A numerical methood of solving Laplace equation by using MA TLAB was discussed ,the main ste ps can be followedas :generate a new available region by editing MA TLAB ’s M 2file ,using numgrid and delsq ,solve the finite difference of Laplace ’s equation with Dirichlet boundary value problem.K ey w ords :Laplace ’s equation ;finite difference ;MA TLAB科技期刊中表格的规范化使用表格简称表,是记录数据或事物分类等的一种有效表达方式,其作用与插图类似,都是代替或补充文字叙述・对表格的设计要求应该是科学、明确、简洁,具有自明性,具体地说就是突出重点、内容简洁、设计科学、表达规范・一个表通常应包括表序、表题、项目栏、表身,必要时可以加表注・即使文章中只有一个表,也应该加上表序,不可省略或用“附表”等・表格有多种形式,在科技期刊中推荐使用三线表・三线表是卡线表的一种,是一般卡线表经简化和改造而成的,通常一个表只有3条线,即顶线、底线和栏目线,其中顶线和底线为粗线,栏目线为细线・三线表不一定只有3条线,必要时可加辅助线・三线表保留了传统卡线表的几乎全部功能,却克服了传统卡线表的缺点,增强了表格的简洁性,所以科技书刊中普遍使用三线表・961第2期 刘大卫等:正方形环域Laplace 方程的简明数值解法。