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

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

第 3卷第4 2期 21 0 0年 4月

红外技术I fa e e h o o y n rdT c n l g r

Vb13 N O. .2 4

Ap . 2 1 r 00

<材料与器件>

利用有限差分和 MA L B矩阵运算直接求解二维泊松方程 T A王忆锋，唐利斌(昆明物理研究所，云南昆明 6 0 2 ) 5 2 3

摘要：根据有限差分法原理，将求解范围用等间距网格划分为一系列离散节点后，二维泊松方程可转化为用一个矩阵方程表示的关于各未知节点的多元线性方程组。利用 MA L B提供的矩阵左除命 TA

令，即可得到各未知节点的函数近似值。该方法概念简单，使用方便，不需要花费较多精力编程即可以求解大型线性方程组。 关键词：半导体；泊松方程；有限差分法；MA L T AB中图分类号：T 0 N3 1文献标识码：A文章编号： 10—8 12 1)40 1—4 0 18 9 (0 00—2 30

Di e tSo uto fTwo di e i na is n Eq to r c l ino - m nso l Po s o ua i n

wih Fi ieDi e e c n ATLAB a rx Co p t to t n t f r n ea d M M t i m u a i nW AN G—e g, TAN G— n Yif n Libi

( u mi stt o P yi, u mig6 0 2, hn ) K n n I tue f hs sK n n 5 2 3 C ia gn i cAbs r c:Ba e n fn t if r n ep n i l s a e o ut e i n i i i e n o a s reso ic ee n de ta t s d o ied fe e c r c p e, f rs l i r g o sd v d d i t e i fd s r t o s i i t onwih a v n y s a e n e v 1 h wo— i e i a is qu ton c n b o et d i t uliee e t t n e e l p c d i t r a.t e t d m nson lPo s on e ai a e c nv re n o m t—l m n

ln a q t s a o n o o swh c a e e p e s d i t x e uai n. e a p o mai n o i e r e uai b utu kn wn n de i h C b x r s e n a ma r q to Th

p r xi to f on n ie c n no n n a u to a e o t i e t a rX l f vii n c a h u k w od lf nci n c n b b a n d wih m ti e tdi so omm a n M ATLAB . ti i p e nd i I ssm l

i o c pt o e i n n o e ai n a d c n b e o s l e l ge ln a q to t o tmo e e ors i n c n e,c nv n e ti p r to n a e us d t o v a i e e uai nswih u r f t n r rp o r m mi g. rga n

Ke r: s mi on u t r Po s o qu ton, fn t i e e c t od, M ATLAB y wo ds e c d co, i s n e ai i ied f r n eme h

摹半导体器件模拟中的静电场问题可归结为在给定电荷分布和边界条件下求解泊松方程…。和薛定谔方程类似 L, 2作为半导体中常用的基本控制方程之一， JV一一Q( e / Or C )

和() 1

式()为电势 1称的泊松方程，式中 Q为电荷量，

C为真空介电常数，为相对介电常数，另外： o V=

泊松方程、尤其是二维以上的泊松方程一般没有解析解，需要根据所用算法编写程序做数值计算。差分法

导熹萼++“

c 2() 3

是常用数值解法之一。本文介绍了一种利用有限差分和 MA L B矩阵运算直接求解二维泊松方程的方法， TA 即将求解区问离散为网格化的节点，通过差分将泊松方程转化为以矩阵方程表示的代数方程组，调用一条矩阵左除命令，即可求出其近似数值解。

为拉普拉斯算子。二维泊松方程可以写为：+a。v

)

考虑如图 1示的边值问题。以 h为步长将求解所区域等间距划分为 NXN个正方形网格，计有 (Ⅳ+1 ) ×( 1个节点，其中待求未知节点( 1×(一 1Ⅳ+ )Ⅳ一 ) J ) V 个，己知边界节点 4 N个。

1二维泊松方程在正方形网格节点上的有限差分在低频时对于介电常数均匀的情况，利用

收稿日期：2 0 .02;修订日期：2 1 - .0 0 9 1—1

0 00 1 . 4 作者简介：王忆锋 (9 3),男，湖南零陵人，高级工程师，主要研究方向为器件仿真。 16一

21 3

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

第 3卷第 4期 2 21 00年 4月

红外技术I fa e c ol y nx r dTe hn og.√ . “

、o _2 No4,l 3 . Ap . 2 0 r 01

: ( 1… 1 )(1 1 ( f ) f 1 ( l 1 - ) f )+ …一。 _——……一 - _— 9 l 1 )

=

=

=

九十

九一-

九 -1九 0 1 /4十一

…

0 0

14/+

1+

-1 4… /

(j 0 t+

—

…

…

-

I

4(十+), l1Ⅳ I l(+Ⅳ I

G=++

( q一 1 - l

(D 4…一一

【

广…

)…

lIl

O 0

一

……

一

() 90 -1 4/ 0 1 -1 4/ -1 4/ 1

▲Y

O, l一 i 16一 -

( l一) M 1 l

矿矿矿妒一十一+一 一 一 妒1 O 1=

( 16……一 6—&— 6…一 (+ 1 1) -,一——— _…一Ⅳ 1),O I1 -, (I ) I),【 I1 r,)+

0

O 0

…

0

图 1有限差分的正方形网格节点 +

1 0 1 ————

0.

O O Z=:●

00

: + 4+●

.●

‘

+.

01

‘ .●

‘

(0 1)

.

F g 1 S u r rd n d so n t i e e c i . q a eg i o e f i df r n e i f e

一对函数，, )做泰勒级数展开，可以写出： )~

一一

0

一… 0一

0

… 0 O

式()的妫所有未知节点呜构成的( 1×I 7中Ⅳ一 ) 本身为由各行未知节点构成的(Ⅳ一1×1 )

阶矩阵， 阶矩阵：

一一 一 .

,

缟= , =

:●

:●

,

J

[:式 ( 1中两个矩阵标注做了调整】注 1) 另外，B为由构成的 ( 1×1阶矩阵，为 Ⅳ一 ) () 4

由相关节点函数值和边界值构成的(

1×1阶矩Ⅳ一 )阵：-

将式 () 4中的各式相加，略去 h 4以上的高阶项，可得：+。

h q 2+ g2l 2,,

+ gl 2,

:

￡h

() 5

::●

,

1=一

-

hq2 3,

+g I 1 3 ;

,

O y

4-

( 2 1)

h2 ., g2 2 qⅣ2+Ⅳ+,

于是式()以写为： 3可-

+ g1 Ⅳ+l,

Ⅵ+￣ l" ', 1 -j-

+Oj1￣卜， -

() 61=

一

q,+gf 2f l,一

^.

在如图 1示的正方形网格节点中，从第二行 (所 =2 )开始，对每一行中的各未知节点按式 ()出相 6写应的差分表达式，经整理后可得下列形式的矩阵方程：=B () 7

i

(< f ) 2<Ⅳ,

-

h q f gⅣⅣ+.

2一

,

Ⅳ+ gⅣ 1+ gⅣ 22+1+。

,

式中：为( 1×( 1阶矩阵，并有： Ⅳ一 )Ⅳ一 ) G I Z… Z Z

1=

-

hqⅣ 3,

+ gⅣ+, 23

i-

(3 1)

h q J gⅣ 1+ 2Ⅳ,Ⅳ+ .+. 2+ gN 2 .Ⅳ+,

N+ l

I G,…

K:

:●

:●

( 2大型矩阵的多软件辅助生成 8 )Z I G I 一 i

Z

…

维泊松方程边值问题的有限差分解法，是利用

Z…

Z I G

式中：G、,、z均为(Ⅳ一1×( 1阶矩阵，其定义 )Ⅳ一 )分别为：

网格状模型上离散节点的数值解来逼近原函数的真实解。在计算机软硬件环境允许的情况下，网格划分越细密，或节点选取越多，离散化模型越能较为精确地逼近真实问题，从而获得具有足够精度的数值解。

21 4

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

第3 2卷第 4期 21 0 0年 4月

、 -2 N O. bl3 4

王忆锋等：用有限差分和 MA L B矩阵运算直接求解二维泊松方程 利 TA

A p . 201 r 0

但相伴而来的是大型矩阵。融合使用 Wod和 r MA L B,可使大型矩阵的构建易于实现。 TA 例如，矩阵的可调用下述程序辅助生成： fn t nK= t x N) u ci ma K( o i fp o e ( 1;= n sN,)

’ ) 1 f+ lN

+,+ l 22)] (=/ ( V ( 11 v (,)’ 1 4 q ) N+;) ds (b i[’ n m2￣ N)’=/ f+ 1N+,+ ) p’ u s ( ( l ( V ( 1 2+ N) 4 l q Nv ( 2N+1)’ lN+, );)]

ds (fr=: 1 ) i[o 2N一’ p’ i]

ds ( i[ p’

bn m2￣ N)(=/ f+ lN+, 1)]’u s (’) 1 ( v ( 2+ )’ i 4 q i;)

KI s da s[*, 2 p, 1,】 N)= p i ( pP ]一， 1, g 2,[ 0 N,;%生成稀疏矩阵

ds (e d ) i[n p’]

此外，b 2<J的计算表达式可调用下列程序 i<i V ( )生成：fn t nb= ix r si nN) u ci i be pe s ( o of ri N o _2:

%[,, p为三条对角线元素；一，,] 2 pP2】 [ 01为三条对角 1线元素位置

K=ul )将稀疏矩阵转换为满矩阵 flK1% (输入下列语句：>>K= t x 4 mar K( 1 i

ds (b n m2 t i’o e ( 1 f/; b n r2￣ i i[ u s (= n s p r) N,) q4 ’ un s ( )’ ) 1 f+ l21) (=/ ( v (,; 1 4 q )e nd

b’

n m2￣ i u s ( )

结果为K=

’=/ f+ 1N+, u s (’;) ( l ( V ( 2’ m2V i )’ N) 4 q n ) )】21 2 O

12 0 0

02 1 2

O0 2 1

3直接计算求解二维泊松方程的程序框架式() 7是一个以矩阵方程表示的线性方程组。线性方程组的具体解法可参见各种文献，例如文献[,】 3 4等，本文不再涉及。这些解法的一个共同点是要根据相应算法，编写和调试程序，因而有一定的工作量。利用 MA L T AB强大的矩阵计算功能，分别定义好矩阵和后，在相当大的范围内，采用如下表达式：=

将复制到一个 Wod文档，将 l替换为 G、2 r替换为 I 0替换为 z再将整个矩阵复制回 MA L 、, T AB

的

文档即可。矩阵 G可直接在程序中写一条

sda s )令生成。 p ig(命

在 MA L T AB程序中，B矩阵可以表示为

:B= b;2 b;4 b;6 b;8 b】【 lb;3 b;5 b;7 b;9;

即可直接计算式() 7。其中反斜线代表矩阵除法，

该矩阵可通过调用下列程序辅助生成：fn t nB ma iB N u ci = tx ( ) o rs m sb y

在 MA L T AB中称之为矩阵的左除，它比通过计算逆矩阵来求解式() 7要具有更好的数值稳定性『。 5】直接求解二维泊松方程的 MA L T AB程序框架如下所示，照不同的 J值补齐相应的部分即可运行该按 7 v程序。cos l; la l; l eal c e ra l

B= n s1芈；= eo (,; o e (, b C z r s1N) N)f ri: o=lN

C(,=; 1) i ie d n

N 9%待求未知节点数为 Nx=; N Z zrs N;%生成 ( 1xN—)= eo( ) N, N一) ( 1阶零矩阵 I一y (/;=e eN) 4p o e( 1;= n sN,)

B=B+C:

例如输入：>= tx (、>B ma iB 9 r

%生成( 1xN一) N.)( 1阶单位矩阵

结果为：B=[+l b 2 b,+,+,+,+,+,+9 b,+,+3 b 4 b 5 b 6 b 7 b 8 b】

G= sda s[/,一/]一，,, N)%生成 N N p ig ( p4P p4, 1 1 N,;一，[ 0] x

阶稀疏矩阵 GK= .

将上述结果复制到 Wo文档，并将其相应替换 d r形成B= bl b; 3 b; 5 b; 7 b; 9;[;2 b;4 b;6 b;8 b]

%调用子程序 mar K N)将其运行结果复制到 tx (, i

Wod中替换、调整后粘贴到此处； r;】;

再复制粘贴回 MA L B的 M文档中即可。 TA

考虑对于不同Ⅳ值的计算情况，MA L T AB程序

vl o e ( 2N+ )= n sN+, 2;f ri:+2 o _1N

v (,= n s 1 2 O l 1: o e (, ) 5; ) N+

中b N的计算表达式可调用下述程序生成：fn t nb b x rsinN) u c o N= Ne pes ( i o ds(bn m2￣N)= n sN,;) i[ u s (’ e( 1 p~ o )】 ds(b i[ p~ n m2￣ N) u s(

v ( 2: o e(, 2*

0;%边界条件 l N+,= n s1 ) N+ ) 10 v ( 1=; v (N+ ) 10%边界条件 l, O i ) li 2= 0;,ed n21 5

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

第3卷第4 2期 21 0 0年 4月

红外技术I fa e e h oo y nrrdT c n lg

、 .2 No 4 b1 3 . Ap . 2 1 r 00

h:=l

q y 0 f= h 2 q y x=; q一 x;

v (k= h(;%将计算结果连同边界条件 l,)p i) j i分解拼成 (+ ) N 1阶矩阵 N 1×(+ )ii: -+1 ed n e nd

b= n s 1; l o e( ) N,f ri N— o=2: l

b () l f+ l1 )v (, ) l1/ ( v (,+ 1 1;= 4 q 2 2 )

b (=/ f+ l1 1+ l2 2) lN) 1 ( v (, ) v (, ); 4 q N+ N+ b (=/ f+ l1i1) 1i 1 ( v (, ) ) 4 q+;ed n

作为比较，考虑文献[】 6中给出的一个算例：一十一+’=。 I I

%调用子程序 b x rsinN) Ne pes (,将其运行结果复制 o粘贴到此处； %调用子程序 be pes nN.)将其运行结果复制 ix rsi ( 1, o粘贴到此处；边界条件为：

a0 v

0 .

X=0 =1 0 Y=0 Y=1 0

%调用子程序 maiBN) tx (,将其运行结果复制到 rWod中替换、调整后粘贴到此处； r p iK B; h= ki;=l

Ox Y=f x y= (,) (, )

1 o. 0 5, 0 1 0 0,

%矩阵左除出各未知节点的近似值求

取 N=9即划分为 1×1=1 0个网格，共有， 0 0 0 1×1=1 1 1 1 2个节点，其中待求节点 9=8个，边 X9 1

界节点 1 1 8=4 2 - 1 0个。执行下述程序得到的结果如表 1所示，其中表格四周为边界条件值；些结果与这文献[] 6中用超松弛迭代法计算得到的结果相同。T be1 C lua drsl fu cinl au r o e sl in ein al aclt uto fnt a v e o d sna ou o go e e s o l f n i t r5 0 4 .1 85 9 4 .0 79 6 4 .8 84 9 5 .5 03 2 5 .9 35

6 5 .5 84 16 28 5. 3 7 5 4.01

fr=: I o j2 N+o f rk=2: N+1

表 1求解区域内各节点函数值的计算结果O 0 0 0 0 0 O0 0

5 0 2 .6 66 l 1 .7 84 5 l .2 56 2 1 .8 50 9 1 .7 59 4 1 .5 82 62 66 2. 3 3135 . 6

5 0 3 .7 8 1 316 9 .1 2 .2 89 4 2 .5 87 9 3 .5 05 3 .8 43 7410 .39 52.06 2

5 0 4 .9 43 8 4 .0 09 5 3 .9 96 8 4 .7 04 2 4 .8 3O l 4 .0 77 45 89 4. 9 65571 .

5 0 51 7 . 7 5 . 371 5 .0 60 1 5 .5 88 6 . 25 6 -2 7273. 28 8 8 0.7

5 0 5 .5 48 3 5 .6 913 6 .5 29 3 6 .4 65 9 7 .3 03 4 7 .4 46 87 7 9.47 8 7 5.58

5 0 5 .7 84 9 6 .3 51 4 7 .0 01 1 7 .6 40 7 .3 76 9 8 .9 12 285. 2 30 8 8 9.

5 0 6 .3 39 7 .9 27 4 7 .5 82 6 8 .4 19 9 8 .7 48 8 .7 75 89 3 0.71 9 38l 3_

5 0 7 .4 44 6 8 .5 38 5 8 .8 81 1 9 .l 06 1 9 .1 23 5 9 .7 37 895.21 2 96. 5 73

10 0 10 0 10 0 10 0 1O 0 l0 0 10 01 00 1o o

0 0

5 .5 05 4 1o o

7 .5 08 9 10 0

8 .7 06 8 10 0

8 2 6-8 10 0

8 .4 99 1 10 0

9 .1 26 4 10 0

9 .5 47 9 10 O

9 .2 66 10 0

9 .3 83 9 1o 0

l0 0 10 0

4结束语在二维程序中，所用网格数一般为 3×3 2 2到 26 5×2 6 5 r。当节点数大于 6 5 6时， T A的 A ry 53 MA L B r a

念直观、结构简单，涉及的编程技巧很少。相应的代

价是程序很长，例如当有 656网格点时， T AB 53个 MA L程序有五百多行，粘贴到 Wod中显示有近九

十页之 r多。但是其中绝大多数是由软件辅助自动生成的矩阵

E i r己无法显示，这时可用类似于 tv (3这样的 dt o= l: ),命令分别显示各行或列的计算结果。若在上述程序中

定义，对于使用者来说只是几次复制、粘贴之类的操作，工作量增加不多，总体上可使二维泊松方程的数

插入相应的计时函数命令，如 c c l k或 cui o pt me等， 可以发现程序的运行时间主要耗费在矩阵准备上，而花在矩阵除法上的时问并不多。 本文介绍的方法充分利用了不同软件的特点，概21 6

值求解过程大为简化。

(下转第 20页) 3

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

第 3卷第 4期 221 O 0年 4月

红外技术I r rd c no o nfae Te h l gy

、 l 2 No 4,- 0 3 . Ap . 2 1 r 00

(上接第 26页 ) 1 参考文献：【】王忆锋，毛京湘.用 MAT A 1 L B和打靶法实现平面 P N结一维泊松方程的简捷计算 f .外，0 03 ()4 .6 J红 1 2 1,02:44 . []王忆锋，唐利斌 .用转移矩阵和 MA L 2利 T AB直接求解一维薛定谔方程的一种简捷方法【 .外技术，0 03 ()1710 j红] 2 1,23:7 .8 .【】马文淦.计算物理学[ .北京：学出版社， 0 5 3 M】科 2 0.

[] E B Marb S Azr B. aahnrn M A L B原理与工程应用 5 . . ga, . am, B l ada. T A c[ .北京：电子工业出版社，0 2 M] 20 . [】何红雨 .电磁场数值计算法与 MA AB实现[ .汉：中科技大 6 TL M】武华学出版社， 0 4 20.

【】邵福球 .离子体粒子模拟[1. 7等 I]北京：科学出版社，02 V 2 0.

【】任玉杰 .值分析及其 MA L B实现【 .北京：高等教育出版社， 4数 T A M】2 7 f}

20 3

本文来源：https://www.bwwdw.com/article/anii.html