利用有限差分和MATLAB矩阵运算直接求解二维泊松方程

第 32卷 第4期 2010年4月 红 外 技 术 Infrared Technology Vb1．32 NO．4 Apr． 2010 < 材料 关于××同志的政审材料调查表环保先进个人材料国家普通话测试材料农民专业合作社注销四查四问剖析材料 与器件 > 利用有限差分和 MATLAB矩阵运算直接求解二维泊松方程 王忆锋 ，唐利斌 (昆明物理研究所，云南 昆明 650223) 摘要：根据有限差分法原理，将求解范围用等间距网格划分为一系列离散节点后，二维泊松方程可 转化为用一个矩阵方程表示的关于各未知节点的多元线性方程组。利用 MATLAB提供的矩阵左除命 令，即可得到各未知节点的函数近似值。该方法概念简单，使用方便，不需要花费较多精力编程即可 以求解大型线性方程组。 关键词：半导体；泊松方程；有限差分法；MATLAB 中图分类号：TN301 文献标识码：A 文章编号：1001—8891(2010)04—0213—04 Direct Solution of Two-dimensional Poisson Equation with Finite Difference and MATLAB M atrix Computation WANG Yi—feng，TANG Li—bin (Kunming Institute ofPhysics,Kunming 650223，China) Abstract：Based on finite difference principles，after solution region is divided into a series of discrete nodes with an evenly spaced interva1．the two—dimensional Poisson equation can be converted into multi—element linear equations about unknown nodes which Can be expressed in a matrix equation．The approximation of each unknown nodal function can be obtained with matriX left division command in MATLAB．It is simple in concept，convenient in operation and can be used to solve large linear equations without more efforts in programming． Key words： semiconductor，Poisson equation， finite difference method，M ATLAB 摹 和 半导体器件模拟中的静电场问题可归结为在给 定电荷分布和边界条件下求解泊松方程⋯。和薛定谔 方程类似L2 J，作为半导体中常用的基本控制方程之一， 泊松方程、尤其是二维以上的泊松方程一般没有解析 解，需要根据所用算法编写程序做数值计算。差分法 是常用数值解法之一。本文介绍了一种利用有限差分 和MATLAB矩阵运算直接求解二维泊松方程的方法， 即将求解区问离散为网格化的节点，通过差分将泊松 方程转化为以矩阵方程表示的代数方程组，调用一条 矩阵左除命令，即可求出其近似数值解。 1 二维泊松方程在正方形网格节点上的有限 差分 在 低频时对 于介 电常数均匀 的情况 ，利用 V 一一Q／(COer) (1) 式(1)称为电势 的泊松方程，式中Q为电荷量， Co为真空介电常数， 为相对介电常数，另外： V =导+熹+萼 c2 为拉普拉斯算子。 二维泊松方程可以写为： + ) (3) 。a v “ 考虑如图 1所示的边值问题。以h为步长将求解 区域等间距划分为 NXN个正方形网格，计有(Ⅳ+1) ×(Ⅳ+1)个节点，其中待求未知节点(Ⅳ一1)×(JV一1) 个，己知边界节点 4N个。 收稿日期：2009．10—21；修订日期：2010-04．10． 作者简介：王忆锋 (1963一)，男，湖南零陵人，高级 工程 路基工程安全技术交底工程项目施工成本控制工程量增项单年度零星工程技术标正投影法基本原理 师，主要研究方向为器件仿真。 213 Lin 铅笔 Lin 铅笔 Lin 铅笔 Lin 铅笔 Lin 铅笔 第32卷 第4期 2010年4月 红 外 技 术 Infxared Technology 、，ol_32 No．4 Apr． 2010 (f-1 1) (f 1) (f+l 1) (1 1) ⋯ ⋯ 一 。__— — — ⋯ ⋯ 一 -9 l 1) (tj+0 — I (1 q-一 l O,i-l1 6一一 (4D 【 广 ⋯ ⋯ -4,(Ⅳ十l +1) l I ⋯ 一 一 (Ⅳ+I ) l I l ⋯ ⋯ (M l 一1) (1，1)6-⋯ ⋯ 一一6———＆——_6⋯ ⋯ 一一 (Ⅳ+1，1) O-I，1) (I，I) 【r+I，1) ▲Y l 图 1 有限差分的正方形网格节点 Fig．1 Square grid nodes of finite difference 对函数 ，)，)做泰勒级数展开，可以写出： (4) 将式(4)中的各式相加，略去 h4以上的高阶项， 可得： G= 1 -1／4 0 ⋯ 0 - 1／4 1 -1／4 ⋯ 0 O ⋯ 0 -1／4 1 -1／4 0 ⋯ 0 -1／4 1 1 1 = —— — — 4 1 O 0 1 ： ． ● ● 0 ⋯ · · · 0 · · · 0 ‘ ． 0 0 1 ． Z= O 0 O O ： ‘． ● ● 0 ⋯ ⋯ 0 · · · 0 ‘ ． 0 0 O (9) (10) 式(7)中的妫 所有未知节点呜构成的(Ⅳ一1) ×I 阶矩阵， 本身为由各行未知节点构成的(Ⅳ一1)×1 阶矩阵： = 缟 ： ● ， = ， ： ● ，J [注：式 (11)中两个矩阵标注做了调整】 另外，B为由 构成的(Ⅳ一1) ×1阶矩阵， 为 由相关节点函数值和边界值构成的(Ⅳ一1)×1阶矩 阵： + ： ￡ (5) ： 。 Oy h 于是式(3)可以写为： -- Ⅵ+~'--l,j"1- +O~,j-1卜 (6) 在如图 1所示的正方形网格节点中，从第二行 ( =2)开始，对每一行中的各未知节点按式(6)写出相 应的差分表达式，经整理后可得下列形式的矩阵方 程 ： =B (7) 式中： 为(Ⅳ一1) ×(Ⅳ一1) 阶矩阵，并有： K ： G I I G ： ： ● ● Z ⋯ Z ⋯ Z ⋯ Z ， ⋯ Z Z I G I Z I G 式中：G、，、z均为(Ⅳ一1)×(Ⅳ一1)阶矩阵， 分别为： 214 ： ● ， 1 ， = 一 4 = 1 = 1 - h q2 ， 2+g2 ， l +gl ，2 - h q3 ，2 +g1I3 ； - h2q．Ⅳ，2+g2 ， Ⅳ+2 + g1 ，Ⅳ+l 一 q2，f+gf，l 一 ^ ．i - h qⅣ ． f+gⅣ (2>K=matrixK(41 结果为 K = 1 2 0 O 2 1 2 0 0 2 1 2 0 O 2 1 将 复制到一个 Word文档，将 l替换为 G、2 替换为 I、0替换为z，再将整个矩阵复制回 MATLAB 的 文档即可。矩阵 G 可直接在程序中写一条 spdiags()命令生成。 在 MATLAB程序中，B矩阵可以表示为： B=【bl；b2；b3；b4；b5；b6；b7；b8；b9】； 该矩阵可通过调用下列程序辅助生成： function B=matrixB(N) symsb B=ones(1，N)芈b；C=zeros(1，N)； fori=l：N C(1，i)=i； end B=B+C： 例如输入： >>B=matrixB(9、 结果为： B=[b+l，b+2，b+3，b+4，b+5，b+6，b+7，b+8，b+9】 将上述结果复制到 Word 文档，并将其相应替换 形成 B=[bl；b2；b3；b4；b5；b6；b7；b8；b9]； 再复制粘贴回MATLAB的M文档中即可。 考虑对于不同 Ⅳ值的计算情况，MATLAB程序 中bN的计算表达式可调用下述程序生成： function bN=bNexpression(N) disp([～b num2s~(N)’=ones(N，1)； 】) disp([～b num2s~(N) ’(1)=1／4 (fq+Vl(N+1，1)+vl(N+2，2))；’]) disp([’b’ num2s~(N) ’(N)=l／4 (flq+V1(N+1，N+2)+ vl(N+2，N+1))；’]) disp([’for i=2：N一1’]) disp([’ b’num2s~(N)’(i)=1／4 (fq+vl(N+2，i+1))；’]) disp([’end ]) 此外，bi(2 表格 关于规范使用各类表格的通知入职表格免费下载关于主播时间做一个表格详细英语字母大小写表格下载简历表格模板下载 四周为边界条件值；这些结果与 文献[6]中用超松弛迭代法计算得到的结果相同。 表 1 求解区域内各节点函数值的计算结果 Table 1 Calculated results offunctional value for nodes in a solution region O 50 50 50 50 50 50 50 50 50 100 0 26．66l 38．17 44．398 48．519 51．77 54．853 58．479 63．93 74．446 100 0 18．475 31．619 40．905 47．906 53．71 59．163 65．134 72．794 83．855 100 0 l5．622 28．924 39．698 48．489 56．001 62．953 70．101 78．256 88．181 100 0 15．089 28．759 40．472 50．352 58．85 66．549 74．06 81．949 90．6l1 10O 0 15．974 30．55 43．O8l 53．596 62．5 70．334 77．639 84．87 92．315 l00 O 18．256 34．387 47．704 58．451 67-22 74．648 81．292 87．578 93．778 100 0 22．663 41．039 54．899 65．283 73．28 79．747 85．302 90．371 95．221 100 0 31．356 52．206 65．571 74．501 80．87 85．758 89．8 93_38l 96．735 1oo 0 50．554 70．859 80．678 86-28 89．941 92．614 94．759 96．62 98．339 l00 0 1oo 100 100 100 100 100 1O0 100 10o 100 4 结束语 在二维程序中，所用网格数一般为 32×32到 256 ×256r 。当节点数大于65536时，MATLAB的Array Editor己无法显示，这时可用类似于 t=vl(：，3)这样的 命令分别显示各行或列的计算结果。若在上述程序中 插入相应的计时函数命令，如 clock或 cputime等， 可以发现程序的运行时间主要耗费在矩阵准备上，而 花在矩阵除法上的时问并不多。 本文介绍的方法充分利用了不同软件的特点，概 216 念直观、结构简单，涉及的编程技巧很少。相应的代 价是程序很长，例如当有65536个网格点时，MATLAB 程序有五百多行，粘贴到 Word中显示有近九十页之 多。但是其中绝大多数是由软件辅助自动生成的矩阵 定义，对于使用者来说只是几次复制、粘贴之类的操 作，工作量增加不多，总体上可使二维泊松方程的数 值求解过程大为简化。 (下转第230页) 第32卷 第4期 2O10年 4月 红 外 技 术 Infrared Technology 、，0l-32 No．4 Apr． 2010 (上接第216页) 参考文献： 【1】 王忆锋，毛京湘．用MATLAB和打靶法实现平面PN结一维泊松方程 的简捷计算fJ1．红外，2010，30(2)：44．46． [2] 王忆锋，唐利斌．利用转移矩阵和 MATLAB直接求解一维薛定谔方 程的一种简捷方法【j]．红外技术，2010，32(3)：177．180． 【3】 马文淦．计算物理学[M】．北京：科学出版社，2005． 【4】 任玉杰．数值分析及其 MATLAB 实现【M】．北京：高等教育出版社， 2f }7 230 [5] E．B．Magrab，S．Azarm，B．Balachandran．MATLAB原理与工程应用 [M]．北京：电子工业出版社，2002． [6】 何红雨．电磁场数值计算法与MATLAB实现[M】．武汉：华中科技大 学出版社，2004． 【7】 邵福球．等离子体粒子模拟[IV1]．北京：科学出版社，2002．