（电磁场与微波技术专业论文）蒙特卡罗法在三维电磁场数值计算中的应用.pdf

声明尸明 本人郑重声明：此处所提交的硕士学位论文蒙特卡罗法在三维电磁场数值计算中 的应用，是本人在华北电力大学攻读硕士学位期间，在导师指导下进行的研究工作和 取得的研究成果。据本人所知，除了文中特别加以标注和致谢之处外，论文中不包含其 他人已经发表或撰写过的研究成果，也不包含为获得华北电力大学或其他教育机构的学 位或证书而使用过的材料。与我一同工作的同志对本研究所做的任何贡献均已在论文中 作了明确的说明并表示了谢意。 学位论文作者签名：荔兰互殓 日期：兰! ! ! ! ! ! ! ! 关于学位论文使用授权的说明 本人完全了解华北电力大学有关保留、使用学位论文的规定，即：学校有权保管、 并向有关部门送交学位论文的原件与复印件；学校可以采用影印、缩印或其它复制手 段复制并保存学位论文；学校可允许学位论文被查阅或借阅：学校可以学术交流为 目的，复制赠送和交换学位论文；同意学校可以用不同方式在不同媒体上发表、传播学 位论文的全部或部分内容。 ( 涉密的学位论文在解密后遵守此规定) 一一 华北电力大学硕十学位论文摘要 摘要 蒙特卡罗法( m o n t ec a r l om e t h o d ，简称m c m ) 是以概率统计理论为基础， 以随机抽样为其主要手段的一种概率统计方法，近些年随着电子计算机的快速发 展，并且由于蒙特卡罗法能够比较逼真地描述事物的特点及物理实验过程，解决一 些数值方法难以解决的问题，因而该种方法的应用领域日趋广泛。本文应用蒙特卡 罗法重点解决了三维静态电磁场边值问题，其中重点解决了第二类边值问题的求 解。在m a t l a b 环境下自行设计编写了所要求解问题的m c m 程序，并且通过典型 算例，并将蒙特卡罗法与有限元法计算结果进行对比，验证了算法及程序的正确性 和可靠性。 关键词：蒙特卡罗法，三维电磁场，数值计算分析，随机游动，概率统计 a b s t r a c t m o n t ec a r l om e t h o di sas t a t i s t i c a lm e t h o dw h i c hi sb a s e do nt h et h e o r yo f p r o b a b i l i t ya n ds t a t i s t i c s ，a n dw h o s ep r i m a r ym e a n si sr a n d o ms a m p l i n g i nr e n c e n ty e a r s w i t ht h er a p i dd e v e l o p m e n to fc o m p u t e r a n db e c a u s em o n t ec a r l om e t h o dc o u l db e m o r er e a l i s t i ct od e s c r i b et h ec h a r a c t e r i s t i c so ft h i n g sa n dp h y s i c a le x p e r i m e n t ，a n d c o u l ds o l v es o m ed i f f i c u l tp r o b l e m sw i c ha r ei n t r a c t a b l et oo t h e rn u m e r i c a lm e t h o d s ，t h e a p p l i c a t i o nf i e l do f t h i st y p eo fm e t h o di si n c r e a s i n g l yw i d e i nt h i sp a p e r , w ef o c u so n s o l u t i o no fo nt h e3 - ds t a t i ce l e c t r o m a g n e t i cf i l e d i np a r t i c u l a rt h es o l u t i o no fn e u m a n n b o u n d a r yv a l u ep r o b l e m i nt h em a t l a be n v i r o n m e n t ，t h ep r o j e c to fm c m s o l v i n gt h e d i r i c h l e tb o u n d a r yc o n d i t i o no fs t a t i cf e l i di sd e s i g n e d ，a n dt h em o n t ec a r l om e t h o di s c o m p a r e dw i t hf i n i t ee l e m e n tm e t h o dt h r o u g ht h et y p i c a le x a m p l e s t h ec o r r e c t n e s sa n d r e l i a b i l i t yo fm o n t ec a r l om e t h o di sv e r i f i e d j i a n gq i a o l i n g ( e l e c t r o m a g n e t i cf i e l da n dm i c r o w a v et e c h o l o g y ) d i r e c t e db yp r o f l i uj i a n g x i n k e yw o r d s ：m o n t ec a r l om e t h o d ，t h r e e d i m e n s i o n a le l e c t r o m a g n e t i cf i e l d ， n u m e r i c a la n a l y s i s ，r a n d o mw a l k , p r o b a b i l i t ya n ds t a t i s t i c s 华北电力大学硕+ 学位论文目录 目录 中文摘要 英文摘要 第一章引言l 1 1 选题背景及意义l 1 2 电磁场数值计算近况1 1 3 蒙特卡罗法发展近况。2 1 4 本论文主要研究内容一4 第二章蒙特卡罗法基础原理5 2 1 用蒙特卡罗法计算的基本思想5 2 1 1 概率统计模型的建立5 2 1 2 随机数抽样6 2 1 3 问题的求解6 2 1 4m c m 计算误差6 2 2 三维问题的蒙特卡罗基本思想7 2 3 随机数的产生一8 第三章三维电磁场蒙特卡罗计算方法理论分析1 0 3 1 蒙特卡罗法解三维拉普拉斯方程1o 3 1 1 区域剖分和方程离散一1 0 3 1 2 概率模型的建立1 1 3 2 蒙特卡罗法解三维泊松方程1 3 3 2 1 区域剖分13 3 2 2 方程离散1 3 3 3 3 概率模型的建立1 4 3 3 本章小结1 5 第四章三维电磁场计算有限元法理论分析o 16 4 1 有限元法解三维电磁场泊松方程第二类边值问题1 6 4 2 有限元法解三维电磁场拉普拉斯方程第二类边值问题计算1 8 4 3 本章小结19 第五章计算程序2 0 5 1m a t l a b 开发环境简介2 0 5 2 程序结构与流程一2 1 5 3 本章小结2 3 第六章典型算例分析2 4 6 1 三维静电场拉普拉斯方程第二类边值问题算例2 4 6 2 三维静电场泊松方程第二类边值问题算例2 6 i 华北电力大学硕士学位论文目录 6 3 蒙特卡洛法影响因素分析2 8 6 3 1 游动步长对m c m 计算的影响2 9 6 3 2 游动次数对m c m 计算的影响3 0 6 3 3 待求点坐标对计算时间的影响3 2 6 4 本章小结3 3 七章结论3 4 7 1 结j 沧3 4 7 2 展望3 5 考文献3 6 j 射3 9 学期间发表的学术论文和参加科研情况4 0 华北电力大学硕十学位论文 1 1 选题背景及意义 第一章引言 电磁场计算的数值方法主要有：有限差分法，有限元法、积分方程法、边界元法和 蒙特卡罗法。其中有限元法自从1 9 6 5 年，w i n s l o w 首先将其应用于电气工程问题，其 后，1 9 6 9 年s i l v e s t e r 将其推广应用于时谐电磁场问题【l 】。发展至今，有限元法凭借容易 求解，收敛性好，占用计算机内存较少，计算程序通用性比较强等优点，已经成为各类 电磁场、电磁波工程问题定量分析与优化设计的主导数值方法。虽然有限元法是现今电 磁场数值计算问题的主导方法，但是仍然存在一定的缺点，这种方法的计算程序复杂冗 长，由于它是区域性解法，分割的元素数和节点数较多，导致需要的初始数据复杂繁多， 最终得到的方程组的元数很大，这使得计算时间长。 蒙特卡罗法也叫随机抽样方法或统计实验方法。半个多世纪以来，由于科学技术的 发展和电子计算机的发明，这种方法作为一种独立的方法被提出来，并首先在核武器的 实验之中得到了应用。蒙特卡罗法是一种计算方法，但与一般数值方法有很大的却别。 它是以概率统计理论为基础的一种方法。由于蒙特卡罗法能够比较逼真地描述事物的特 点及物理实验过程，解决一些数值方法难以解决的问题，因而该方法的应用领域日趋广 泛2 吲。 与有限元法等常用的数值方法相比，蒙特卡罗法有这样几个优点：( 1 ) 收敛速度与 问题维数无关，也就是说，要达到同一精度，用蒙特卡罗选取得的点数与维数无关，计 算时间仅与维数成正比。但是一般数值方法，比如在计算多重积分时，达到同样的误差， 点数与维数的幂次成正比，即计算量要随维数的幂次方而增加，这一特性决定了对多维 问题的适用性。( 2 ) 受问题的条件限制影响小。( 3 ) 程序结构简单，在电子计算机上实 现蒙特卡罗时，程序结构清晰简单，便于编制和调试。( 4 ) 对于模拟像粒子运输等物理 问题具有其他数值计算方法不能替代的作用。 由于蒙特卡罗法的这些优点，现在这种方法已经在核科学、统计物理学、医学、材 料学、经济学等各种领域都得到了广泛的应用。但是在电磁场领域中只得到了初步的研 究，仍然存在许多有待解决的问题。 1 2 电磁场数值计算近况 从1 8 6 4 年m a x w e l l 建立了统一的电磁场理论，建立了著名的m a x w e l l 方程组以来， 所有的宏观电磁问题都可归结为m a x w e l l 方程组在各种边界条件下的求解。在2 0 世纪 6 0 年代前的经典电磁学阶段，电磁场理论和工程中的许多问题多采用解析或渐进法处 理，得到解析解。这类方法能够获得精确解且计算效率比较高，但它的适用范围较窄， l 华北电力火学硕十学位论文 仅能求解有规则边界的简单问题，对于形状复杂边界则无能为力或者需要较高的数学技 巧。 6 0 年代后，随着矩量法、差分法和有限元法等数值方法的广泛应用标志着计算电磁 学时代的到来。它综合了电磁场理论、数值计算方法和计算机软件技术，使得解决各种 类型的复杂问题和大型数值计算成为可能口1 0 | 。 在电磁场数值计算方法中，有限差分法是应用最早的一种方法。有限差分法具有简 单、直观的特点，在电磁场数值分析领域得到广泛的应用。虽然现阶段电磁场数值计算 方法发展很快，有限差分法以其固有的特点仍然是一种不可忽视的数值分析方法。 矩量法是一种将连续方程离散化成代数方程组的方法，它既适用于微分方程，又适 用于积分方程。矩量法只有在计算机可供使用以后才能得到广泛的应用，同时也得到更 快的发展。 有限元法以变分原理为基础，是所以它广泛应用于拉普拉斯方程和泊松方程所描述 的各种物理场中。用有限元法解非线性场及分层介质中的电磁场都有其特点，处理比较 简单。并且这种方法不受场域边界形状的限制，且对第二类、第三类及不同媒质分界面 的边界条件不必作单独处理，因而用有限元法就电磁场数值解有很多优点。 边界元法是近十年来发展形成的数值计算方法，它计算精度高，数据处理工作量小， 适合于解无限域问题，但在处理多介质问题和复杂线性问题等方面不如有限元法方便。 1 3 蒙特卡罗法发展近况 蒙特卡罗方法思想的提出可以追溯到十八世纪著名的蒲丰问题，1 9 7 7 年，法国 科学家蒲丰( b u f f o n ) 提出用投针试验计算圆周率万的问题，但是投针试验求得万的 近似值的方法，是进行真正的试验，并统计实验结果，要使获得的频率值和概率值 偏差小，就要进行大量的试验，在实际中难以做到。可以设想，若要进行l o 万次 投针试验，每次试验用时5 秒，则需要5 0 万秒，大约1 3 9 小时，对于简单的概率 模型需要耗费这么长时间，那么对于实际的例如核裂变等问题，用这种方法近似求 出是不可能的，所以在现代计算机出现之前，这种雏形时代的蒙特卡罗法并没有得 到实际上的应用。 表1 1 用蒙特卡罗法计算的万的近似值 投针次数1 0 0 0 02 0 0 0 01 0 0 0 0 02 0 0 0 0 0 万的近似值3 1 6 2 2 3 33 1 3 7 9 9 33 1 4 1 1 7 93 1 4 1 3 5 4 由表1 1 表明，随着模拟次数的增加，所得到的值越来越接近真值，而要进行 这样的模拟需要很大的计算量，只有通过计算机才能实现，所以直n - 十世纪四十 华北电力人学硕士学位论文 年代以后，随着电子计算机的迅速发展，利用计算机模拟概率过程，实现多次模拟 试验并统计计算结果，进而可获得所求问题的近似结果。计算机的大存储量，高运 算速度使得在短时间内获得精度高且内容丰富的模拟结果，使得蒙特卡罗法在实际 中得到迅速的发展和应用。在第二次世界大战中蒙特卡罗法首先被美国的科学家应 用于原子弹的研制中。蒙特卡罗法的发展大体分为三个阶段。 第一阶段主要发展直接仿真的方法，也就是对具有随机性质的问题进行直接仿 真，没有形成一种普遍和有效的解决实际问题的方法。这段时间中主要的成果的发 现和论证随机变量的统计性质和对随机变量进行定量分析的理论和方法。 蒙特卡罗方法的第二阶段主要发展了数学仿真方法，是蒙特卡罗方法的真正确 立的阶段。数学仿真方法就是针对具体的随机性问题，把对随机现象进行直接仿真 改为用数学方法进行仿真的解答，这种方法显然比直接仿真法有很多的优越性。 蒙特卡罗方法的第三阶段主要发展了模拟仿真的方法，该方法使蒙特卡罗法迅 速发展和广泛的应用。人们为了更加有效的解决实际问题，可以认为地将所要解决 的问题变换为与原问题答案有一定关系的人为“随机现象 ，然后利用“人为随机 现象与原来问题的关系，再利用计算机仿真方法得到原问题的解答。 目前，这种方法被广泛应用于核科学、统计物理学、医学、材料学、经济学等， 学科中，在我国已经成功举办了九届全国蒙特卡罗方法及其应用学术交流会，并且。 成立了学者无涯等蒙特卡罗学术交流论坛，但是在电磁场计算领域的应用相对滞 后，特别是在多媒质电磁场、三维电磁场领域发展非常缓慢。 蒙特卡罗法在电磁场领域的应用是在1 9 4 4 年，k a k u t a n i 提出电位理论和随机游 动之间的关系。在1 9 5 6 年，m e m u l l e r 在求解l a p l a c e 方程边值问题随时正式提 出了不用差分近似的蒙特卡罗法。在1 9 9 0 年，m n o s a d i k u 首次将蒙特卡罗法引 入电磁场数值计算课本教学中，开创了电磁场数值计算和蒙特卡罗法的结合的起 点。1 9 9 5 年，翟景春等建立了l a p l a c e 方程边值问题的概率模型，并且分析了蒙特 卡罗法的误差和收敛速度，并且给出了计算机实现的方法。同年，唐象能首次应用 蒙特卡罗法计算了长方体区域内的静电场边值问题。1 9 9 6 年，m a t t h e wn o s a d i k u 和k e m i n gg u 提出了三角形网格的随机游动模型，解决了椭圆形方程的混合边值 问题，通过二维算例验证了该方法的效率。2 0 0 1 年，马汶金在计算物理书采用离散 型蒙特卡罗法以二维单一媒质静电场为例，说明了蒙特卡罗法的应用，指出了蒙特 卡罗法是求解部分节点数的一种有效的方法。2 0 0 6 年吕英华在计算电磁学一书中在 以二维l a p l a c e 方程边值问题为例，说明了蒙特卡罗法在电磁场数值计算中的应用。 2 0 0 7 年，武卫革在其硕士论文中，以l a p l a c e 第一类边界为例应用蒙特卡罗方法解 决了二维静电场边值问题的应用，和多媒质边值问题的蒙特卡罗法 蒙特卡罗方法的精度问题一直是蒙特卡罗方法受关注的问题。对于这一问题 3 华北电力大学硕士学位论文 1 9 9 2 年，m n o s a d i k u 和d h u n t 提出了一种新的概率游动模型用于计算 d i r i c h l e t 边值问题，并且该模型能在较少的时间达到更高的精度要求。1 9 9 5 年， 唐象能等首次应用蒙特卡罗方法计算了长方体区域内静电场的边值问题，并且提供 了精度估算公式。2 0 0 0 年，裴鹿成研究了蒙特卡法与问题维数的课题，指出蒙特卡 罗方法与一般的数值方法不同，其误差与问题维数无关。同年，赵涛在其硕士论文 中用蒙特卡罗法和有限元法对流速数值求解，用蒙特卡罗方法得出的结果和物理模 型试验结果比较吻合，其误差比有限元法得到结果的误差小很多，体现出了蒙特卡 罗方法的可靠性和优越性。 现在，蒙特卡罗法仍然在进一步发展中，特别是随着非线性物理学的发展，蒙 特卡罗法已经成为了必备的定量分析研究的工具和仿真模拟的理论基础。 1 4 本论文主要研究内容 ( 1 ) 本论文回顾介绍了电磁场数值计算的发展近况，列举并介绍了电磁场数值计算 各种方法。介绍了蒙特卡罗法的基本原理和用蒙特卡罗法解决问题的基本步骤。 ( 2 ) 在三维静电场中，推导了三维静电场中拉普拉斯方程第二类边值问题的蒙特卡 罗解法，并在此基础上推导出三维静电场中泊松方程第二类边值问题的蒙特卡罗解法。 【3 ) 介绍了电磁场数值计算中典型方法有限元法解三维静电场拉普拉斯方程、 泊松方程第二类边值问题解法。 ( 4 ) 在m a t l a b 环境下设计编写了典型算例的蒙特卡罗法计算程序，在a n s y s 环 境下的有限元法对典型算例的计算。 ( 5 ) 将两种方法的计算结果进行比对分析，得出结论。 4 华北电力大学硕士学位论文 第二章蒙特卡罗法基础原理 蒙特卡罗法( m o n t ec a r l om e t h o d ) 是一种根据待求随机问题的变化规律，或者物 理现象本身的统计规律，人为地构造一个合适的概率模型，依照该模型进行大量的统计 实验，使它的某些统计参量正好是待求问题的解。在计算机上通过某种数字进行的假想 “试验 的得到事件出现的频率，或用这个随机变量具体取值的算术平均值作为问题的 近似解。 2 1 用蒙特卡罗法计算的基本思想 蒙特卡罗方法以随机模拟和统计试验为手段，是一种从随机变量的概率分布中，通 过选择随机数的方法产生一个符合该随机变量概率分布特征的随机数数列，将其作为输 入变量序列进行求解的方法。蒙特卡罗方法的基本原理基于统计学中的大数定理和中心 极限定理。大数定理反映了大量随机数之和的性质。如果函数厂在a , b 1 区间内，以均 匀的概率分布密度随机地取刀个数，对每个u i 计算出函数值厂化) ，大数定理告诉我 们这些函数值之和除以一所得的值将收敛于函数厂的期望值，即： 1 l1 t i m 一芝：f ( u ，) = i f ( u ) d u n - - k 。o 甩一i = l b - a 。 ( 2 一1 ) 公式( 2 1 ) 的左边正是公式右边积分的蒙特卡罗估计值。根据这个法则我们在抽取足 够多的随机样本之后，计算得到的积分的蒙特卡罗估计值将收敛于该积分的正确结果 7 1 。 用蒙特卡罗法解决问题时，需要几个过程【1 1 1 5 】： ( 1 ) 描述或构造概率过程，对于本身就具有随机特性的问题，主要是正确的描述和模 拟这个概率过程，对于确定性问题，就必需事先构造一个人为的概率过程，使它的某些 统计参量是所求问题的解，即要将不具有随机性质的问题转化为随机性质的问题，这是 蒙特卡罗法研究与应用的重要组成部分。 ( 2 ) 实现从已知概率分布抽样，构造了概率模型之后，由于各种概率模型都可以看作 是由各种各样的概率分布构成的，因此产生已知概率分布的随机变量就成为实现蒙特卡 罗法模拟试验的基本手段，这也是蒙特卡罗法被称为随机抽样的原因。 ( 3 ) 建立各种统计量，使其期望值是所要求解问题的解。在计算机上通过随即变量具 体取值的算术平均值作为问题的近似解。最后根据所构造的概率模型编制计算程序进行 计算，获得计算结果。 2 1 1 概率统计模型的建立 概率统计模型的建立是蒙特卡罗法的基础，根据实际问题建立一个简单且便于实现 5 华北电力大学硕士学位论文 的概率统计模型，确定随即变量f ，使得随机变量的数学期望值e ( 孝) 等于所要求解问题 的解缈，即 缈2e ( 孝) ( 2 2 ) 本文所要研究问题为不具有随机性的确定性问题，所以事先必须构造一个认为的概 率过程。 2 1 2 随机数抽样 为了实现概率过程的数字模拟，对该模型的随机变量建立抽样方法，抽取足够的随 机数缶乞磊靠，进行模拟试验。 2 1 3 问题的求解 用磊磊磊磊的算术平均值痧作为所求量缈的估计值 缈庐= 专善缶 根据大数定理，得 因此，当n 充分大时，有 p ( 伊2 牌痧) = 1 缈叫伊痧= 专善缶 这说明可以用算术平均值痧作为所求量矿的估计值。 2 1 4m c m 计算误差 ( 2 - 3 ) ( 2 4 ) ( 2 5 ) 痧为一个算数平均值，缶乞磊磊是互相独立，则弓焉( 仃为均方差) 服从数 学期望为0 ，方差为i 的标准正态分布，根据中心极限定理，当足够大时，对于任意 五 0 ，有： 兄) _ 击驴捌 陆6 ， 其中口为显著水平，1 一口为置信水平，这说明蒙特卡罗法与实际值之间的误差范围 6 华北电力大学硕士学位论文 力盯 缈一例 丽 ( 2 7 ) 的概率为1 一口。 口与力的关系可以通过式( 2 6 ) 求得，也可以通过数学用表求得。例如：取置信水平 1 1 一口= 9 9 ，可以查得t = 3 ，可以解释为不等式i 痧一伊i 熹成立的概率是9 9 。 v v 由上可知，为了减小蒙特卡罗法计算误差，就应该选择最优的随机变量，使得其方 差最小。在方差固定时增加模拟次数，可以有效地减小误差。如果试验次数增加1 0 0 倍， 精度提高l o 倍，但是这样做就会增加计算机的计算时间，提高费用。所以在考虑蒙特 卡罗法精确度时，不能只是简单的减少方差和增加模拟次数，还要同时顾及耗费时间【1 6 】。 2 2 三维问题的蒙特卡罗基本思想 ( 工，j ，z + ) ( j ，j ，+ ， 户 h , y ，z ) ( 工，y ，z )( 工+ h ， ， j x , y i t , z ) 图2 一l 三维场剞分相邻结点不葸图 场内一点m 与其相邻的六个节点的关系如图2 一l 与图2 2 所示。为求解内部节点m 的点位值缈( m ) ，构造概率模型如下【1 7 1 9 1 ：考虑场域内一节点在区域内作随机游动，设 质点从m 出发，第一步以的概率向与节点m 相邻的六个网格节点游动，到达其中一 个节点后如果该节点仍是内部节点，与第一步相似，质点继续以的概率向与之相邻 的六个节点游动，直到质点到某一个边界点q 时才停止游动，如果质点本来就在边界节 点上，则质点不游动。按照下面的方式定义一个随机变量孝，设随机游动任意一条游动 路线为 ：m j m m 2 专专心一l 专q 1 1 ( 2 7 ) 7 华北电力人学硕士学位论文 定义随机变量孝的取值为： 图2 - 2 三维场剖分结点示意图 孝_ g ( 匕) = 厂( q ) ( 2 - 8 ) 若质点出发点是边界点q ，则停留在边界点q 处，定义随机变量f 的取值为： 孝2g ( v o ) = 厂( 9 ( 2 7 ) 可以证明这样定义的随机变量孝的数学期望值 e ( f ) 2f p ( m ) ( 2 8 ) 一次随机游动可以得到一个变量孝的值，经过n 次随机游动可以得到n 个随机变 量磊，彘，磊氕，用平均值 醋1 善n 磊( i = 1 , 2 , 3 - - n ) ( 2 - 9 ， 作为伊( m ) 的近似值。 2 3 随机数的产生 蒙特卡罗法是一种概率统计方法，对于不具有随机性的确定性问题的解决，随机数 的产生是利用蒙特卡罗方法求解的基础之一。随机数可以分为3 种不同的类型：真随机 数、准随机数、伪随机数。需要指出的是在实际应用中的随机数通常都是通过某些数学 公式计算而产生的伪随机数，这种伪随机数从数学意义上讲已经不是随机的了，但是只 8 华北电力大学硕十学位论文 要能够通过随机数的一系列的统计检验，我们就能将它当作真随机数而放心的使用。这 些统计性检验包括啪1 ：均匀性检验、独立性检验、组合规律检验、无连贯性检验、参数 检验等等。其中最基本的是均匀性检验和独立性检验。所谓均匀性是指在 0 ，1 区域中 等长度区间的随机数分布的个数应该相等。独立性是按先后顺序出现的若干个随机数 中，每一个随机数的出现都和它前后的各个随机数没有关系。 随着蒙特卡罗方法的不断广泛的应用，随机数的产生方法不断的更新。最古老的方 法就是手工方法，可想而知这种方法耗时巨大，在进行大量的随机模拟试验时显然不能 满足需求。随着电子计算机的广泛利用，在计算机上安装一台物理随机数发生器，把具 有随机性质的物理过程变换为随机数，这样得到的随即数随机性和均匀性都很好，完全 符合随机数的特性。但是这种方法得到的实验结果不能重复，很难推广到不具有随机 性的确定性问题中。 还有另夕卜一种利用计算机来产生随机数的数学方法，是目前使用最广泛，发展最快 的方法，这种方法产生的随机数，即为前面所讲的伪随机数，如果能够通过一系列的统 计检验，即如果能具有真正随即数的一些统计性质，也就可以把它们当作真正的随机数 放心使用了。 本文研究的三维电磁场属于不具有随机性的确定性问题，需要将不具有随机性质的 问题转化为随机性质的问题，在本文中即为对所求区域剖分，构造质点的随即游动，转 化为随机性问题，随要一个随机数列来决定质点的游动路线。本文以m a t l a b 为卡发环 境，在m a t l a b 产生随机数的语句很简单，提供了2 3 种随机变量分布类型的随机数发生 器，如正态分布、对数正态分布、泊松分布、威布尔分布等( 这些分布基本包括了工程实 际中出现的变量分布情况) ，可直接产生变量x 以代入功能函数，省去了可能会带来很 大麻烦的求分布函数反函数这一步，极大地提高了效率。常用指令举例如下昭妇： r = r a n d ( m ，n ) 产生m 行n 列的( o 1 ) 间均匀分布随机数组r 。 r = n o r m m d ( m u ，s i g m a , m ，n ) 产生服从n ( p ，0 2 ) 分布的m 行n 列随机变量数组r 。 r = l gn r n d ( m u ，s i g m a , m ，n ) 产生l nr 服从n ( 0 2 ) 分布的m 行n 列随机变量数组r 。 本文应用r a n d 语句产生( o 1 ) 间均匀分布的随机数组，来决定质点游动路线，将要 解决的确定性问题转化为随机性问题。 9 华北电力大学硕+ 学位论文 第三章三维电磁场蒙特卡罗计算方法理论分析 蒙特卡罗法已经在核科学、统计物理学、医学、材料学、经济学等各种领域都得到 了广泛的应用，但是在计算电磁学领域发展缓慢，特别是三维电磁场的计算，本章重点 介绍了蒙特卡罗方法解三维电磁场中拉普拉斯方程和泊松方程的第二类边值问题的基 本思想，为程序的设计编写奠定了理论基础瞳羽。 3 1 蒙特卡罗法解三维拉普拉斯方程 静电场三维拉普拉斯方程第二类边值问题为： v 2 = 窘+ 窘+ 窘= 。c 旷c v 。c 幺。 卦他) 烈q ) = 甜( 奶) p q q e f i qe f 2 ( 3 - l a ) ( 3 一坳 ( 3 1 0 其中q 为求解区域，f , f ：为q 的两个边界，0 l 是边界i ，上的点，0 2 是r ：上的点， p 是q 中的内点。 3 1 1 区域剖分和方程离散 将q 等步长立方体剖分，步长为h ，在接近边界处，取与边界最近的结点，剖分后 的边界为i 。剖分后区域内任意一结点与相邻六个网格结点的关系如图( 3 1 ) ，则公式的 差分形式为【2 3 之9 】： 伊( x ，y ，z ) = e + 伊( x + ，y ，z ) + 一伊( x 一 ，少，z ) + + 驴( 毛y + h ，z ) + e f o ( x , y h ，z ) + p + 缈( x ，y ，z + 矗) + 一地y ，z 一| i i ) ( 3 2 ) 其中：巳、只，厶、墨、已、一为随机游动的粒子在点o ，y ，z ) 向相邻结点o + j i i ，y ，z ) 、 o h ，y ，z ) 、( z ，y + j l i ，z ) 、o ，y - h ，z ) 、o ，y ，z + 五) 、o ，y ，z j 1 1 ) 移动的概率，并且 + = 巴= o + = 弓一= + = 一= 丢( 3 - 3 ) 1 0 华北电力大学硕士学位论文 ( 五j ，z + | 1 1 ) ( t y + h ， ， h , y ，z ) ( 工，y ，z )( x + h j 1 ， x , y 一 ，z ) 图3 - 1 三维场剖分相邻结点示意图 图3 - 2 三维静态场区域剖分节点示意图 3 1 2 概率模型的建立 如图3 2 所示，为求解内部结点p 的电位值，构造概率模型如下：可以构造一个随 机游动：质点从区域内结点p 出发，每一步以1 6 的概率向相邻结点游动，到达其中一 个结点后，如果该节点仍是内部节点，与第一部相似，质点仍以1 6 的概率向与之相邻 的六个节点游动【3 0 3 2 1 。 l 、当它达到第一类边界r ：上的点q 2 时停止游动。定义一个随即变量f ，并且设随 即游动的任意一条路线为【3 2 】 y 。：p 专m l _ 鸠寸一m k i 专q 2 r 2 定义随即变量f 的取值为 1 1 华北电力大学硕士学位论文 f = g ( v v ) = “( q 2 ) 这样定义的随即变量的数学期望值即为内点p 的电位值，即： e 皓) = 伊( 一 ( 3 q ( 3 5 ) 一次随机游动得到一个f 的取值，则经过n 次随即游动就得到n 个随即变量 磊，磊，氏，用平均值 芋= 2 ：i l i i ；n ；j ；：( 1 = = 1 ，2 ，3 。】n ) ( 3 6 ) 作为矽( p ) 的近似值。 2 、当它达到第二类边界r 。上的点q 时，或者游动的最初点就是q i 时，它不一定立 即终止游动，而是以概率p 回到与q 相邻的结点上，继续进行游动，并以概率l p 终 止在q 点。p 是满足从0 到1 的任意实数【3 3 】。 按照下面的方式定义一个随机变量g ，设随机游动任意一条游动路线为： 咋：p 专m i 专m 2 专专m 纠- - ) q f j 对于上述游动定义一个随机变量孝，在每一条自m 点出发的游动路线中，当质 点由区域内某个结点游动到q 时，定义善的值为 善= g ( ) = 南他) ( 3 - 7 ) 当到了边界点q 又回到内结点p e 继续进行游动，则定义善的值为 孝2 9 ( ) 2 万1g ( ) ( 3 8 ) 这就是说，游动又从以处开始，最后将终止于边界点，但孝对于自点p 处出发的路线 的值，应等于孝对于原来就从点以处出发的路线k 的值与1 p 的乘积1 2 彤1 。 这样定义的随机变量的数学期望值即为p 点的电位值缈( p ) ，即： e ( 孝) = 缈( p )( 3 9 ) 1 2 华北电力大学硕+ 学位论文 3 2 蒙特卡罗法解三维泊松方程 对于任一三维电场问题，如电位函数缈= 认x , y ，z ) ，来描述场的分布，则电位缈通 常应满足如下的边值问题 3 2 1 区域剖分 p ( 五y ，z ) = 一一 g = f o ( p ) 2 = 0 ( 3 - 1 0 a ) ( 3 - 1 0 b ) ( 3 - 1 0 c ) 对q 的剖分采用正六面体剖分，步长为a x = a y = a z = h ，任葸区域内一点m 【x ，y ，z ) 与其相邻的六个点的关系与图3 - 1 一致。 3 2 2 方程离散 方程( 3 10 a ) l 拘差分形式为： 妒( z ，y ，_ z ) 一篓p ( 五y ，z ) ：只+ 认x + j i i ，j ，z ) + 只一烈石一五，少，z ) + e + 认x ，j ，+ j i l ，z ) o 占 + e 一烈五y 一 ，z ) + + 烈x ，y ，z + 而) + e 一烈工，y , z - h ) ( 3 1 1 ) 其中 只+ = 只一2e + 2 e 一= 十= 一2 否1 ( 3 1 2 ) 设函数 ( z ，y ，z ) ：缈( 石，y ，z ) 一竺p ( x ，y ，z ) o 占 ( 3 一1 3 ) 则式( 3 1 1 ) 表示为 ( x ，y ，z ) = 否1 缈( x + 办，y ，z ) + 否1 缈( x 一办，y ，z ) + i 1 伊( x ，y + j l l ，z ) + 否1 缈( x , y - h , z ) + 否1 缈( x , y , z + h ) + l f a ( x , y , z - h ) ( 3 1 4 ) 边界条件( 3 10 b ) 变化为 吨= 咀一万h 2 九= 五一石h 2 以 ( 3 1 5 ) 边值条件( 3 - 1 0 c ) 变化为 一 华北电力人学硕十学位论文 其中 当p 为常数时 得到 筹妒警i 岛一等等妒。一z = 。3 舶， 石= 等瓢 石= 0 抛 ，、 丽k 。u ( 3 - 1 7 ) ( 3 - 1 8 ) ( 3 - 1 9 ) 这就将泊松方程形式转化为和前面公式( 3 1 ) 相似的形式，在程序中的体现为根据边 界条件( 3 - 1 5 ) n 1 ( 3 1 9 ) ，先用蒙特卡罗方法计算出，根据公式( 3 1 3 ) 得到缈。 3 3 3 概率模型的建立 随机游动的路线与前面拉普拉斯方程蒙特卡罗解法相同。定义随机变量孝设随机游 动任意一条游动路线为【3 4 1 咋：m ijm 2 专专m 川一q f 1 、在每一条自m 点出发的游动路线中，当质点由区域内某个结点游动到第一 类边界条件的边界面j 。上的边界点q 时停止游动，定义随即变量孝的取值为： 孝= g ( 1 ，p ) = e a ( q ) 多= g 1 ，。) =踢 这样定义的随机变量的数学期望值即为内点户的电位值，即： e 偕) = ( p ) ( 3 - 2 1 ) ( 3 2 1 ) 一次随机游动得到一个善的取值，则经过n 次随即游动就得到n 个随即变量 螽，受，磊，用平均值 作为o ( p ) 的近似值。 手= 2 j ；1 i ：萎n j ；：( 1 = ：1 ，2 ，3 。】q ) ( 3 2 2 ) 1 4 华北电力大学硕士学位论文 2 、在每一条自m 点出发的游动路线中，当质点由区域内某个结点游动到q 2 时， 定义善的值为 孝_ g ( ) = 击f ( q 2 ) ( 3 - 2 3 ) 当到了边界点q 2 又回到内结点见继续进行游动，则定义f 的值为 孝2g ( ) 2 1g ( 嘲 ( 3 2 4 ) 这样证明随机变量f 的数学期望值 由公式( 3 1 3 ) 可以得到 3 3 本章小结 e ( 孝) = e ( g ( 巧) ) = ( p ) 缈( 尸) = ( 一+ 万h 2p ( 尸) ( 3 - 2 5 ) ( 3 - 2 6 ) 本章主要介绍了三维电磁场蒙特卡罗法，首先在蒙特卡罗法理论分析的基础上，计 算三维拉普拉斯方程第二类边界条件和三维泊松方程第二类边界条件的蒙特卡罗法，包 括区域的剖分、方程的离散、概率模型的建立、随机游动位置确定等。在其求解过程中 可以看出：蒙特卡罗法可以单独求解计算区域内我们感兴趣点的点位值，其电位值不需 要通过联解其他几点电位值求解，只需要输入边界条件即可，这是蒙特卡罗法的优势之 一。此外，蒙特卡罗法求解过程中只涉及到几个节点的计算，不需要计算刚度矩阵，所 用计算容量较小，这是蒙特卡罗法又一优势。 1 5 华北电力大学硕士学位论文 第四章三维电磁场计算有限元法理论分析 有限元法是电磁场数值计算方法之一，凭借容易求解，收敛性好，占用计算机内存 较少，计算程序通用性比较强等优点，已经成为各类电磁场、电磁波工程问题定量分析 与优化设计的主导数值方法。本章内容主要介绍有限元法解决三维电磁场拉普拉斯方程 和泊松方程第二类边值问题的推导。 三维场有限元法，就是将空间剖分成有限个互不重叠且无裂缝的空间单元，在三维 场离散化时，剖分的单元经常采用四面体、五面体、六面体等各种形式，本文采用了正 六面体剖分。构造空间单元的插值函数和形状函数，然后利用等价变分问题，得到空间 单元的系数矩阵，并累加成总系数矩阵，最后计算有限元方程，即可得到场的数值解口和矧。 4 1 有限元法解三维电磁场泊松方程第二类边值问题 对于任一三维电场问题，如电位函数伊= 钗x ，y z ) ，来描述场的分布，则电位缈通 常应满足如下的边值问题 所对应的等价变分问题为 卜州2 饼饼卜批一 【q s 2 = f o ( p ) 采用和蒙特卡罗相同的正立方体剖分，则其差值函数为 缈2 + 口2 x + a 3 y + a 4 z 4 a 5 砂+ a 6 y z + 口7 z x + x y z 则可以得到单元内部的点位函数值 o ( x ，y ，z ) = 孵仍 其中形状函数m 分别为 1 6 ( 4 一l a ) ( 4 1 6 ) ( 4 - l c ) ( 4 _ 2 ) ( 4 - 3 ) ( 4 4 ) 丝 ， 、i “ 垆k 叽儿翱 t f，；，牲，j铲j科们带o j ，m 华北电力人学硕士学位论文 吖：( h - x ) ( h i - y ) ( h 一- z ) h 孵：x ( h - y _ ) - ( h 一- z ) 2 h 孵= 丁x y ( h - z ) 弘羞芝 卜5 ， 肟= 学 一 孵= 学 n ；= i x y z 孵= 丁( h - x ) y z 用矩阵形式表不电位函数的表达式： 驴( x ，y ，z ) = 【l 【缈l ( 4 - 6 ) 将整个区域q 离散成个元素n o 个结点其能量泛函，( 缈) 可表示为各个元素e 的能量 泛函c ( 妒) 的总和 ，( 纠圭f ( 动= e ( 刃 卢1 ( 4 一7 ) 变分问题( 3 i ) 被转化为一多元二次函数的极值问题 r o f ：羔要：o ( 2 ，氓) a 仍智a 仍 、 一。 ( 4 8 ) 将痧( x ，y ，z ) = 【l 【伊l 代入泛函 t ( 矽) 圭c ( 乒) 州( 州荆针叫蚴 ， 并对仍求导，可得 ( 4 一l o ) 2 = o = 7 p一 缈 k = 1j幔鬲 。l 泊松方程( 4 1 a ) 式中p ( x ，y ，z ) 为0 ，方程即为拉普拉斯方程，即 v 2 舯，= 等+ 守+ 窘= 。 孰川q ) 缈( q 2 ) = “( q 2 ) 与此等价的变分问题为 1 8 p q q f t q r 2 ( 4 1 4 a ) ( 4 -