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

第 32卷 第4期 2010年4月 红 外 技 术 Infrared Technology Vb1．32 NO．4 Apr． 2010 <材料与器件 > 利用有限差分和 MATLAB矩阵运算直接求解二维泊松方程 王忆锋 ，唐利斌 (昆明物理研究所，云南 昆明 650223) 摘要：根据有限差分法原理，将求解范围用等间距网格划分为一系列离散节点后，二维泊松方程可 转化为用一个矩阵方程 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 示的关于各未知节点的多元线性方程组。利用 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