PN结泊松方程的一种改进算法及其Matlab验证

邓宏贵，郭  俊，刘利强，周克省

(中南大学 物理科学与技术学院，湖南 长沙，410083)

摘  要：针对在PN结泊松方程求解过程中几种常用方法存在的不足，提出一种改进算法。该算法结合求解非线性方程组的Newton迭代法与SOR(逐次超松弛迭代)法，即用松弛因子对Newton迭代过程的前、后2项进行加权平均，组成新的迭代公式。为进一步完善算法，在迭代公式中修改松弛因子，采用最佳松弛因子形式。根据改进算法的计算思路，运用Matlab7.0编程，对算法进行仿真与模拟。结果表明：算法真实可行，既保持计算的高精度，也明显地减少计算的迭代次数，提高求解过程的收敛速度，且仿真图像与文献图像较吻合。

关键词：泊松方程；Newton迭代法；逐次超松弛迭代法

中图分类号：O475       文献标识码：A         文章编号：1672-7207(2008)05-0913-05

An ameliorative arithmetic of Poisson equation for PN junction and its verification of Matlab

DENG Hong-gui, GUO Jun, LIU Li-qiang, ZHOU Ke-sheng

(School of Physics Science and Technology, Central South University, Changsha 410083, China)

Abstract: There exist some disadvantages in several common methods for solving the Poisson equation of PN junction, and so a new ameliorative arithmetic was presented. It combined Newton iterative method for solving nonlinear equations and SOR (Successive over-relaxation) method, which means weighted average between the two closed items in Newton iterative method. Finally the relaxation factor was revised to the best relaxation factor so as to improve the arithmetic better. According to the thinking of computing, the ameliorative arithmetic was simulated and modeled by Matlab 7.0 to validate its feasibility. The results show that this arithmetic is authentic and feasible. It keeps the high precision, advances the convergence speed markedly and decreases the iteration times of computing significantly as well. What’s more, the image of simulation is consistent with that of literature well.

Key words: Poisson equation; Newton iterative method; successive over-relaxation method

随着集成电路的飞速发展，对半导体器件的计算机模拟变得越来越重要^[1]，尤其是数值模拟方法在半导体器件的研究和生产中得到了越来越广泛的应   用^[2]。半导体器件的数值模拟一般需要求解3个耦合的差分方程：泊松方程、电子连续性方程和空穴连续性方程^[3]。其中，泊松方程是一类椭圆型偏微分方程，在流体力学、渗流理论、电磁学、传热学、结构力学等学科的基础理论研究中占有重要的地位。

针对半导体PN结在强电磁脉冲作用下的响应，国内外已经有不少相关研究^[4-7]。泊松方程反映了PN结各点的电势分布，在半导体PN结的物理特性研究中也有着举足轻重的作用。 但是，由于只有少量的泊松方程有解析解，因此，寻求该类方程高精度的数值解已成为理论研究的重要内容^[8-10]。

泊松方程是非线性耦合方程，这类方程可以通过选择合适的数值减幅迭代求解^[11-14]。PN结泊松方程的常用解法有3种：第1种是有限元法。对于不同的初始条件，有限元法容易发散，其数据处理很复杂；第2种方法先用Newton法将方程线性化，再求解3对角线矩阵方程，最后用迭代法求解线性方程的准确解，称Newton法Ⅰ^[15]。此方法计算思路简单，但编程复杂、计算量大、占用贮存空间大，具有局部收敛性。第3种方法结合有限差分和Newton迭代公式进行求解，称Newton法Ⅱ^[15]。这种计算方法除编程简单、节省贮存单元外，还有全收敛的优点，但其收敛速度低，导致消耗CPU较多，编程计算时间过长。基于这些不足，本文作者提出将Newton法Ⅱ与SOR法相结合的计算方法，既保留上述几种方法的优点，也克服它们的缺陷。在保证计算精度的前提下，减少计算过程中的迭代次数，提高方程求解的收敛速度，降低CUP消耗，节省编程耗时。

1  PN结的物理模型

本文采用的PN结为线性缓变PN结^[15]，其杂质分布可表示为：

其结构图如图1所示^[16]。

图1  线性缓变PN结杂质分布

Fig.1  Doping distribution of linear ease variable PN junction

由半导体物理理论，线性缓变PN结的Poisson方程为：

2  Poisson方程的改进算法

由于用有限差离散法解Poisson方程的过程，不受V取值的影响(为简单起见，令V=0)，方程(2)可改写为：

在线性缓变PN结取L个等间距网格点，则Poisson方程的有限差分式为：

两端网格点的值由边界条件决定^[17]：

对式(4)进行线性化处理，Poisson差分方程变成关于PN结电势的线性方程组：

令

，       (8)

则

。     (9)

根据Newton迭代公式，可得Poisson方程Newton法的计算公式：

由上述求解过程可以看出，用Newton迭代法求解泊松方程思路清晰，计算简单。而且根据Newton迭代法的收敛定理^[18]，它是全局收敛的。但是，这种方法存在一个严重缺陷：计算过程中迭代次数多，收敛速度慢。

为了克服Newton迭代法的这一缺陷，本文对Newton法进行改进，将其与SOR法相结合，修改松弛因子，形成新的算法。

根据SOR法的求解思路，在计算出的第k+1个近似解基础上，用参数ω对与进行加权平均，即

式(13)即本文改进算法的计算公式。在改进算法的整个计算过程中，添加和修改松弛因子的目的在于减少计算的迭代次数，加快算法的收敛速度，节省计算机程序运行时间。

3  改进算法的编程验证

根据改进算法的计算思路与公式，应用Matlab 7.0编程求解PN结Poisson方程，对改进算法的可行性进行验证。由于本文中PN结两侧是完全对称的，所以，程序只考虑x≥0一侧的情况。

表1给出了Newton法与改进算法在误差限为10^-5的计算精度下，求解Poisson方程所需的迭代次数。可见，改进算法显著地减少了迭代次数，实现了计算收敛速度的提高。由此证实了改进算法的可行性。

表1  Newton法与改进算法在求解Poisson方程过程中的迭代次数

Table 1  Iteration times of Newton Method and ameliorative arithmetic in process of solving Poisson equation

表1中改进算法的松弛因子ω₀=1.25。可见，改进算法大幅度地缩短了程序的计算机运行时间，提高了方程求解的运算速度。

对于不同的计算精度，Newton法Ⅰ、Newton法Ⅱ和改进算法所用的迭代次数各有不同。计算精度越高，改进算法减少迭代次数的优势越明显，如表2所示。

表2 不同计算精度下几种方法的迭代次数

Table 2  Iteration times of several methods for different computing precisions

表3给出了在相同计算精度和不同松弛因子下改进算法的计算机耗时和迭代次数。由表3可见，收敛速度与松弛因子并没有明显的比例关系。通过多次取值得出，在某一合适的松弛因子附近，收敛速度达到最高，此时的松弛因子称为最佳松弛因子。本文得到的最佳松弛因子ω₀=1.35。

表3  不同收敛因子下改进算法的计算机耗时

Table 3  Computing time of ameliorative arithmetic for different convergence factors

图2所示为在改进算法下运用Matlab7.0编程得出的线性缓变PN结的电势分布图和文献[15]中线性缓变PN结的电势分布图。根据半导体理论及文献曲线，PN结中的电势随着位置从P型到N型逐渐上升，呈非线性增加。图2中曲线2体现了这一变化趁势，表明本程序的仿真结果真实可靠。由于结点的选取有差别，导致改进算法中的初始值比文献[15]中初始值要低，所以，图中曲线2的位置比曲线1的位置略低。

图2  改进算法与文献所得PN结电势分布比较

Fig.2  Comparison between potential distribution in PN junction from program of ameliorative arithmetic and potential distribution in PN junction from Ref.[15]

4  结  论

a. 采用等间距网格，将问题简化，易于求解和编程。

b. 采用Newton法Ⅱ的迭代思路，保留了编程简单、占用贮存量小以及全收敛的优点。

c. 所提出的改进算法结合SOR法，引入和修改松弛因子，弥补了其他方法迭代次数多、收敛速度低、CPU消耗较大、计算耗时长的缺陷。在保持高精度的情况下，显著地提高了收敛速度，缩短了程序的计算机耗时。

d. 此改进算法仅在一维空间得到实现，把它应用到二维、三维空间将成为后续研究工作的发展方向。

参考文献：

[1] 刘 珺. 半导体器件计算机模拟新方法研究[D]. 南宁: 广西大学电气工程学院, 2005.
LIU Jun. A study of a novel method for semiconductor device simulation[D]. Nanning: School of Electrical Engineering, Guangxi University, 2005.

[2] 韩峰岩. 半导体器件模拟中求解载流子方程的混合算法[J]. 空军工程大学学报, 2001, 2(2): 41-44.
HAN Feng-yan. A hybrid algorithm for solving carrier transport equations in semiconductor device[J]. Journal of Air Force Engineering University, 2001, 2(2): 41-44.

[3] Mayergoyz I D. Solution of the nonlinear Poisson equation of semiconductor device theory[J]. J Appl Phys, 1986, 59(1): 195-199.

[4] 余 稳, 蔡新华, 黄文华, 等. 电磁脉冲对半导体器件的电流模式破坏[J]. 强激光与粒子束, 1999, 11(3): 355-358.
YU Wen, CAI Xin-hua, HUANG Wen-hua, et al. The current-mode destroy of semiconductor devices by electromagnetic pulse[J]. High Power Laser and Particle Beams, 1999, 13(3): 355-358.

[5] 余 稳, 蔡新华, 黄文华, 等. 二极管失效和烧毁阈值与电磁波参数关系[J]. 强激光与粒子束, 2000, 12(2): 215-218.
YU Wen, CAI Xin-hua, HUANG Wen-hua, et al. Relationships of diode failure and burnout energy thresholds with RF parameters[J]. High Power Laser and Particle Beams, 2000, 12(2): 215-218.

[6] Orvis W J, McConaghy C F, Yee J H, et al. Modeling and testing for second breakdown phenomena[C]//Electrical Overstress/ Electrostatic Discharge Symposium. Las Vegas, NV, 1983: 108-117.

[7] Yee J H, Orvis W J, Martin L C, et al. Modeling of current- and thermal- mode second-breakdown phenomena[C]//4th Annual Electrical Overstress/Electrostatic Discharge Symposium. Orlando, Florida, 1982: 76-81.

[8] 田振夫. 求解泊松方程的紧致高阶差分方法[J]. 西北大学学报: 自然科学版, 1996, 26(2): 109-114.
TIAN Zhen-fu. Compact high-order difference methods for solving the Poisson equation[J]. Journal of Northwest University: Natural Science Edition, 1996, 26(2): 109-114.

[9] 田振夫. 数值求解Poisson方程的四阶紧致差分格式[J]. 宁夏大学学报: 自然科学版, 1995, 16(3): 22-24.
TIAN Zhen-fu. A kind of H⁴ compact difference schemes for solving Poisson equation[J]. Journal of Ningxia University: Natural Science Edition, 1995, 16(3): 22-24.

[10] 曹卫东, 陆昌根. 泊松方程非等距有限差分的数值求解方法[J]. 河海大学学报: 自然科学版, 2006, 34(2): 123-126.
CAO Wei-dong, LU Chang-gen. Numerical solution of Poisson equation with non-uniform mesh finite difference scheme[J]. Journal of Hehai University: Natural Science, 2006, 34(2): 123-126.

[11] Pourfath M, Kosina H, Selberherr S. A fast and stable Poisson-Schrodinger solver for the analysis of carbon nanotube transistors[J]. J Comput Electron, 2006, 5(2): 155-159.

[12] 周一峰, 唐进元, 何旭辉. 轴向受力梁强非线性超谐波与次谐波共振的能量迭代法[J]. 中南大学学报: 自然科学版, 2005, 36(4): 698-703.
ZHOU Yi-feng, TANG Jin-yuan, HE Xu-hui. Energy-iterative method for solution of super-harmonic and sub-harmonic resonance of strongly non-linear vibrations with an axially forcing beam[J]. Journal of Central South University: Science and Technology, 2005, 36(4): 698-703.

[13] 曹 平, 邓志斌, 陈 枫. 边坡稳定性任意条分法分析中安全系数的计算[J]. 中南大学学报: 自然科学版, 2005, 36(2): 302-306.
CAO Ping, DENG Zhi-bin, CHEN Feng. Calculation of safety factor for slope stability analysis with polygon elements[J]. Journal of Central South University: Science and Technology, 2005, 36(2): 302-306.

[14] 周一峰. 求强非线性动力系统周期解的MH法和EMH法[J]. 中南大学学报: 自然科学版, 2004, 35(1): 151-156.
ZHOU Yi-feng. The MH and EMH method for strong nonlinear dynamic systems[J]. Journal of Central South University: Science and Technology, 2004, 35(1): 151-156.

[15] 赵鸿麟. 半导体器件计算机模拟[M]. 天津: 天津大学出版社, 1989.
ZHAO Hong-lin. Compute simulation of semiconductor devices[M]. Tianjin: Tianjin University Press, 1989.

[16] 施 敏. 半导体器件物理[M]. 黄振岗, 译. 北京: 电子工业出版社, 1987.
Simon M S. Physics of semiconductor devices [M]. HUANG Zhen-gang, translates. Beijing: Electronic Industry Press, 1987.

[17] Orvis W J, Yee J H. A one-dimensional code for modeling solid state devices[R]. San Francisco: Lawrence Livermore National Laboratory, 1985.

[18] 韩旭里, 万 中. 数值分析与实验[M]. 长沙: 中南大学出版社, 2006.
HAN Xu-li, WAN Zhong. Numerical analysis and experiment[M]. Changsha: Central South University Press, 2006.

通信作者：邓宏贵(1965-)，男，湖南邵阳人，博士后，副教授，从事高功率微波与半导体器件相互作用机理研究；电话：13975159948；E-mail: denghonggui@163.com