2020年南邮计算物理实践报告

　　南 南 京 邮 电 大 学 实

　　验 验 报 告 课程名称:

　　计算物理实践

　　专 专

　　业:

　　应用物理学

　　学 学

　　号:

　　姓 姓

　　名:

　　完成日期:

　　2014 年 7 月

　　南邮计算物理实践报告

　　目

　　录 第一章

　　简单物理实验的模拟及实验数据处理 1

　　1、1 问题描述1 1、2 原理分析1

　　1、2、1 特殊情况1

　　1、2、2 一般情况3 1、3Matlab 程序仿真4 1、4Matlab 仿真结果4 第二章

　　方程组的解 5

　　2、1 问题描述5 2、2 原理分析5

　　2、2、1 迭代公式的建立及其几何意义5

　　2、2、2 解题过程5 2、3 流程图6 2、4Matlab 程序仿真6 2、5Matlab 仿真结果6 第三章

　　静电场问题的计算 7

　　3、1 问题描述7 3、2 原理分析7 3、3Matlab 程序仿真9 3、4Matlab 仿真结果9 第四章

　　热传导方程与波动方程的差分解法 10

　　4、1 问题描述10 4、2 原理分析10 4、3 解题步骤13 4、4Matlab 程序仿真13 4、5Matlab 仿真结果13 第五章

　　矩量法在静电场边值问题计算中的应用 16

　　5、1 问题描述16 5、2 原理分析16 5、3Matlab 程序仿真18 5、4Matlab 仿真结果18

　　结束语 19

　　参考文献 20

　　附录一 21

　　附录二 22

　　附录三 23

　　附录四 25

　　附录五 26

　　南邮计算物理实践报告 第一章

　　简单物理实验的模拟及实验数据处理 1、1 问题描述 模拟电偶极子的场与等位线。

　　设在 ) , ( b a 处有电荷 q ,在 ) , ( b a 处有电荷 q 。那么在电荷所在平面上任何 一 点 的 电 势 与 场 强 分 别 为 )1 1(4) , (0 r rqy x V, V E 。

　　其 中2 2 2 2) ( ) ( , ) ( ) ( b y a x r b y a x r ,9019 104 。

　　又 设 电 荷62 10 q , 5 . 1 a , 5 . 1 b 。

　　1、2 原理分析 电偶极子就是指一对等值异号的点电荷相距一微小距离所构成的电荷系统,它就是一种常见的场源存在形式。

　　1、2、1特殊情况 图(1)表示中心位于坐标系原点上的一个电偶极子,它的轴线与Z轴重合,两个点电荷q 与-q 间的距离为L。此电偶极子在场点 P 处产生的电位等于两个点电荷在该点的电位之与,即

　　)1 1(4) (0 r rqr

　　(1) 其中 r 与 r 分别就是q 与-q 到 P 点的距离。

　　图(1)

　　电偶极子

　　一般情况下,我们关心的就是电偶极子产生的远区场,即负偶极子到场点的距离r 远远大于偶极子长度L的情形,此时可以的到电偶极子的远区表达式

　　204c o s) (rqlr

　　(2) 可见电偶极子的远区电位与 ql 成正比,与r的平方成反比,并且与场点位置矢量r与z轴的夹角β有关。

　　为了便于描述电偶极子,引入一个矢量 p ,模为qL ,方向由-q 指向q ,称之为此电偶极子的电矩矢量,简称为偶极矩,记作

　　p ql

　　(3) 此时(2)式又可以写成

　　20204 4cos) (rperqlrr

　　(4) 电偶极子的远区电场强度可由(4)式求梯度得到。因电位

　　只就是坐标r 与β的函数,于就是有

　　2 30 0cos sin2 4rp pE e er r

　　(5) 从(4)式与(5)式可以瞧到,电偶极子的远区电位与电场分别与r的平方与r的三次方成反比。因此,其电位与场强随距离 的下降比单个点电荷更为迅速,这就是由于两个点电荷q与-q的作用在远区相互抵消的缘故。

　　根据(4)式,电偶极子的等电位面方程可由 204cos) (rqlr

　　为定值得到。

　　将电力线微分方程写成球坐标形式,并注意此时电场只有r与 两个分量,则有:

　　c crErdEdr

　　(6) 把电场表达式(5)带入上式,得:

　　 22sin) (sinsincos2d drdr

　　(7) 解上式得:

　　南邮计算物理实践报告 Cr 2sin1

　　(8) 式(8)即就是电偶极子远区场的电力线方程。

　　图(2)绘出了电偶极子 为常数的平面内(8)式取不同的常数所对应的等电位线与电场线。

　　图(2)

　　电偶极子的场与等电位线 说明:图中准确的只就是电力线的形状,电力线的疏密并不严格与场强成正比,只就是疏的地方场强小些,密的地方场强大些而已。

　　1、2、2一般情况 前面讨论了电偶极子的中点位于坐标系原点且偶极矩方向为Z的情况。对于中点不在原点与偶极矩非Z的方向的一般情况,通过与前面类似的推导,可以得到远区的电位: 204cos) (rqlr

　　(9) 其中,r就是电偶极子中心指向场点P的相对单位位置矢量,偶极矩P=qL,L的方向依然规定为从-q到q 。

　　经推导还可得到远区场的电场强度表达式:

　　03 00304sin4cos 2rprrpV E

　　(10) 由上式可以瞧出,电偶极子的电场线均分布于由r、θ构成的平面上,并且任意一个平面上的电场线分布都相同。

　　从以上几种不同情况下电偶极子在空间激发的电场结果来瞧,电场强度与p 成正比,与源点到场点的距离3r 成反比,电偶极子在远处的性质就是由其电偶极矩

　　来表征的,电偶极矩就是电偶极子的重要特征。

　　设电荷所在平面上任意一点的电势为

　　)1 1(4) , (0 r rqy x V

　　(11) 其中

　　2 2 2 2) ( ) ( , ) ( ) ( b y a x r b y a x r

　　(12) 因此,只要给定空间任意一点的位置坐标P(x,y),就可以算出这一点的电位。

　　1、3Matlab 程序设计仿真 源程序见附录一 1、4Matlab 仿真结果

　　第二章

　　方程组的解法 2、1问题描述 用牛顿法解方程 1 0xxe ,精度自设。

　　2、2原理分析

　　2、2、1 迭代公式的建立及其几何意义 (1)建立公式 将 (x) f 在nx 点 Taylor 展开 """ 2(x )(x) (x ) (x )(x x ) (x x ) ...2!nn n n nff f f

　　"(x) (x ) (x )(x x )n n nf f f ——Taylor 展开线性化 (x) 0 f 近似于"(x ) (x )(x x ) 0n n nf f

　　解出 x 记为1 nx,则1"(x )(x )nn nnfx xf

　　(n=0,1,...) （2）几何意义

　　过 ( , ( ))n nx f x 切线 ( ) "( )( )n n ny f x f x x x 与 0 y 求交点,解出1 nx x ,则1"(x )(x )nn nnfx xf

　　2、2、2 解题过程 令 1 ) ( xxe x f ,有x xxe e x f ) ( " ,那么根据 Newton 迭代法建立迭代公式 1"(x ) 1(x )xnn n nx xnf xex x xf e xe 2、3流程图

　　N Y 开始 x0=0、5 e=0、0001 00 00001xx xx ex xe x e x-x0>e

　　2、4Matlab程序设计仿真 源程序见附录二 2、5Matlab仿真结果 x=0、5671

　　第三章

　　静电场问题的计算 3、1问题描述 长直接地金属槽,如图 3-2 所示,其侧壁与底面电位为零,顶盖电位为x sin 100 ,求槽内电位,并绘出电位分布图。

　　3、2原理分析 (1)原理分析: 二维拉普拉斯方程

　　) , ( ) , (2y x f y xyy xx

　　(1) 有限差分法的网格划分,通常采用完全有规律的分布方式,这样可使每个离散点上得到相同形式的差分方程,有效的提高解题速度,经常采用的就是正方形网格划分。

　　设 网 格 节点 (i,j) 的电 位 为j i,, 其 上 下 左右 四 个 节点 的 电 位分 别 为。

　　， ， ，j i j i j i j i , 1 , 1 1 , 1 , 在 h 充分小的情况下,可以j i,为基点进行泰勒级数展开: 333222, 1 ,6121hyhyhyj i j i

　　 333222, 1 ,6121hyhyhyj i j i 333222, , 16121hxhxhxj i j i 333222, , 16121hxhxhxj i j i 把以上四式相加,在相加的过程中,h 的所有奇次方项都抵消了。得到的结果的精度为 h 的二次项。

　　2 22, 1 , 1 1, 1, ,2 24 ( )i j i j i j i j i jhx y

　　(2) 由于场中任意点 ( , ) i j 都满足泊松方程: 2 222 2=( , ) F x yx y 式中 ( , ) F x y 为场源,则式(2)可变为:

　　2, , 1 , 1 1, 1,1( ) ( , )4 4i j i j i j i j i jhF x y

　　(3) 对于无源场, ( , ) 0 F x y ,则二维拉普拉斯方程的有限差分形式为:

　　) (41, 1 , 1 1 . 1 , , j i j i j i j i j i

　　(4) 上式表示任一点的电位等于围绕它的四个等间距点的电位的平均值,距离 h越小则结果越精确,用式(4)可以近似的求解二维拉普拉斯方程。

　　边界条件: ( , ) 0(0, ) (a, ) ( ,0) 0( ,b) 100sinxx yyx yy y x Vx xV

　　(2)解题过程: 在直角坐标系中,金属槽中的电位函数 满足拉普拉斯方程: 2 22 20x y 其边界条件满足混合型边值问题的边界条件: ( , ) 0(0, ) (a, ) ( ,0) 0( ,b) 100sinxx yyx yy y x Vx xV 取步长 1 h , , x y 方向上的网格数为 16, 10 m n ,共有160个网孔与17 11 187 个节点,其中槽内的节点(电位待求点)有 15 9 135 个,边界节点52个,设迭代精度为610 ,利用MATLAB编程求解。

　　3、3Matlab程序设计仿真 源程序见附录三 3、4Matlab仿真结果

　　第四章

　　热传导方程与波动方程的差分解法 4、1问题描述 求有限空间内的热传导问题:2 22 2( , ,0)u u ut x yu x y xy 的数值解,边界条件如教材中图9、2所示,其她参数可以自取,将计算结果图形化。

　　4、2原理分析 二维热传导方程的初、边值混合问题与一维的类似,在确定差分格式并给出定解条件后,按时间序号分层计算,只就是每一层就是由二维点阵组成,通常称为网格。

　　内部无热源均匀介质中二维热传导方程为: 2 22 2( )u u ut x y

　　( 0 x l

　　0 y s

　　0 t T )

　　(1) 其初始条件为: (x,y,0) (x,y) u

　　(2) 现在设时间步长为 ,空间步长为 h ,如图9、3所示,将 xOy 平面均分为M N 的网格,并使 Nh l

　　Mh s 则有: t k

　　0,1,... k

　　x ih

　　0,1,2,...,N i

　　y jh

　　0,1,2,..., j M

　　对节点 ( , ) i j ,在 k 时刻(即 k 时刻)有: , , , , 1 , ,2, , 1, , , , 1, ,k2 22, , , 1, , , , 1,2 222i j k i j k i j ki j k i j k i j k i ji j k i j k i j k i j ku u utu u u ux hu u u uy h

　　(3) 将差分格式(3)代入偏微分方程(1)中,可得: , , 1 ,j,k 1, , 1, , , 1, , 1,(1 4 ) ( )i j k i i j k i j k i j k i j ku u u u u u

　　(4) 式中2h

　　式(4)为二维热传导方程的显式差分格式,运用式(4)与边界条件就可以由初始条件逐次计算出任意时刻温度的分布。

　　下面讨论边界条件: 如图9、3所示阴影部分,即在 0 x 边界的10 y M h 与2M h y Mh 区域以及整个 x Nh , 0 y Mh 边界均为绝热壁;而在 0 x 边界的1 2M h y M h 区域为与恒温热源相连的口。

　　0 y 与 y Mh 两边界温度始终为0,实际上也就是与恒温源相连的。也就就是说,对于绝热壁应满足: 0, ,0j kux (1 21,2,..., 1, 1,..., 1 j M M M

　　1,2,3,... k ) , ,0N j kux

　　( 1,2,..., 1 j M

　　1,2,3,... k ) 上述边界条件的差分近似式为: 1, , 0, ,0j k j ku uh

　　, , 1,j,k0N j k Nu uh

　　即:0, , 1, , j k j ku u

　　(1 21,2,..., 1, 1,..., 1 j M M M

　　1,2,3,... k )

　　, , 1, , N j k N j ku u

　　( 1,2,..., 1 j M

　　1,2,3,... k )

　　(5) 对于与恒温源相连的边界,在热传导过程中始终有恒定的热流,常可取归一化值,例如高温热源可取“1”,而低温热源可取“0”。按图9、3的情况,边界条件还有: 0, ,1j ku

　　1 2,...,M j M

　　,0,k , ,0i i M ku u

　　0,1,2,..., i N

　　综合上述初值、边值混合问题,并设初始时刻各点温度均为零,则上述差分格式可归纳为: , , 1 , , 1, , 1, , , 1, , 1,, ,00, , 1, , 1 2, , 1, ,,0, , ,(1 4 ) ( )0 0,1,2,..., ; 0,1,2,...,1,2,..., 1, 1,..., 1; 1,2,3,...1,2,..., 1; 1,2,3,...i j k i j k i j k i j k i j k i j ki jj k j kN j k N j ki k i M ku u u u u uu i N j Mu u j M M M ku u j M ku u 0, , 1 20 0,1,2,...,1 ,...,j ki Nu j M M

　　(6) 可以证明,对于二维热传导方程,若满足 214 h

　　则差分格式式(4)或式(6)就就是稳定的差分格式,一般的讲,对于n维抛物线型微分方程差分格式稳定的充分条件就是: 212 h n

　　4、3解题步骤 给定 、 h 、 与 T 以及 XN 与 YM ,题目中已知 0.5 h ,14 , T 的值分别取0s,10s,100s,120s,150s,200s与1000s, XN 与 YM 取18与16; 计算XNNh 为36;YMMh 为32;2h 为0、05; k 的上界T; 计算初值与边值:, ,0( , )i ju ih jh ;0, , 1 (, )j ku g k jh ;, , 2 (, )N j ku g k jh ; ,0, 3 (, )i ku g k ih ;, , 4 (, )i N ku g k ih ;

　　用差分格式计算, , 1 i j ku; 4、4Matlab程序设计仿真 源程序见附录四 4、5Matlab仿真结果 通过Matlab画出0s 到1000s 之间的一些温度场的分布图,如下图4、1—图4、7分别为0s,10s,100s,120s,150s,200s,1000s的温度场分布图。

　　结论:很明显可以瞧出,温度呈整体下降的趋势。由于低温热源的范围比高温热源的更大,所以热量的流入大于流出。可以断定,只要时间足够长,整个温度场除高温热源外,其她地方的温度都要与低温热源相同(设为0)。1000s 时,如图4、7所示的场分布与无限长时间之后的场分布就已经很接近了。

　　图4、1

　　0s时的场分布

　　图4、2

　　10s时的场分布

　　图4、3

　　100s时的场分布

　　图4、4

　　120s时的场分布

　　图4、5

　　150s时的场分布

　　图4、6

　　200s时的场分布

　　图4、7

　　1000s时的场分布

　　第五章

　　矩量法在静电场边值问题计算中的应用 5、1问题描述 利用矩量法求无界空间中边长为2a的正方形导电薄板的电容。

　　5、2原理分析

　　一块正方形导体板,如上图所示。边长为 2a 米,位于 z=0 平面,中心坐标在原点,设 ( , ) x y 表示导电板上面电荷密度,板的厚度为零,则空间任意一点的静电位就是

　　( x , y , z )0( , )4a aa ad dR

　　(1) 式中2 2 2 1/2[(x ) (y ) z ] R , ( , ) 为待求的面电荷密度。

　　边界条件:0(x,y,0)

　　( , x a y a ) 导体板电容:1( , )a aa aqC d dV V 算子方程:00( , )( )4a aa ad d LR 算子:014a aa aL d dR （1）将导体板分为 N 个均匀小块nS ,并选基函数为分域脉冲函数。

　　1Nn nnp

　　其中1 S0nnnpS 在 上在其它 上

　　(2) 将式(2)代入式(1)得1Nmn nnV l

　　m=1,2,3,…N

　　(3) 式中2 214 ( ) ( )a amna al d dx y 据此电荷密度由逼近,平行板电容相应地近似为: 111Nn n mn nn mnC S l SV

　　(4)

　　若令 2 2 / b a N 表示的边长,由nS 本身面上的单位电荷密度在其中心处产生的电位就是: 2 21 2(0.8814)4b bmnb bbl dx dyx y （2）用点匹配法选权函数为 (x x ) (y y )m m mw , (x ,y )m m为mS

　　的中心点,求内积: , (p ) (x x ) (y )L(p )dxdymn m n m m nx ay al w L y

　　2 21(p )|4 (x ) (y )mn nmn n r rm ml L d d

　　(5) mnl 就是nS 处单位均匀电荷密度( 1np )在mS 处中心的电位。

　　0 0, (x x )(y y )m m m mx ay ag w g dxdy

　　01.1mg 式(5)适用于 m n 时mnl 的求解,当 m=n 时 2 2001 2ln(1 2)4b bmnb bbl d d (6) 其中22abN

　　(3)矩阵求逆解得: 1n mn ml g 1Nn nnp

　　5、3Matlab程序设计仿真 源程序见附录五 5、4Matlab仿真结果 当边长 2a=10 时,电容 C=7、9556e-010 由公式推导可知:C 的变化与 a 成正比; 有实验验证可知:C的变化也与a成正比。

　　结束语 经过这次计算物理学实验周的学习,我认识到自己对于以前学习过的一些课程掌握得还不够透彻,Matlab 编程语言的运用也不够熟练。通过这次实验也很好的巩固了以前学习的一些知识点,并且使我了解了如何利用计算机来模拟与计算一些物理问题。这次实验让我认识到数理方程的实用性,掌握了利用差分代替微分来求解波动方程、热传导方程、拉普拉斯方程等的基本原理与方法。

　　本次实践涉及到的二维拉普拉斯方程以及二维热传导方程的解题方法,都就是先将连续的方程以及边界条件离散化,再用计算机进行计算,因为计算机智能对离散的数值进行计算。

　　对于非线性方程的求解往往就是采用迭代的方法求解,本次实践主要涉及了Newton 迭代法的重要思想,也就是将连续的方程离散化后再进行计算。

　　矩量法主要分为三个步骤:(1)离散化;(2)取样检测;(2)矩阵求逆;适用于场源分布不确定的情况,用未知场的积分方程来计算给定媒质中的场的分布。

　　这次的实践,使我对 Matlab 的使用变得熟练了,并且在报告的写作过程中也熟练掌握了数学公式的录入,文章的排版等技能。

　　总的来说,这次实践带给了我很多的收获。

　　参考文献 [1]陈锺贤、计算物理学、哈尔滨工业大学出版社、2001、3 [2]杨振华,郦志新、数学实验、科学出版社、2010、2 [3]林亮,吴群英、数值分析方法与实验:基于MATLAB实现、高等教育出版社、2012、9

　　[4]李庆杨,王能超,易大义、数值分析、华中科技大学出版社、2006、7 [5]钟季康,鲍鸿吉、大学物理习题计算机解法—MATLAB编程应用、机械工业出版社、2008、1 [6]何红雨、电磁场数值计算法与MATLAB实现、华中科技大学出版社 附录一: close all; clear; clc; k=9e+9;

　　e_p=2e-6;

　　e_n=-e_p; d=-10:0、1:10; [x, y]=meshgrid(d);%产生格点矩阵 a=1、5,b=-1、5; x_n=-a; y_n=-b; x_p=

　　a; y_p=b;

　　V1=

　　k * e_n 、/ sqrt((x-x_n)、^2 + (y-y_n)、^2);

　　V2=

　　k * e_p 、/ sqrt((x-x_p)、^2 + (y-y_p)、^2);

　　V1_min=k * e_n /0、1; V2_max=k * e_p /0、1; V1(V1==-Inf)=V1_min;

　　V1(V1<V1_min)=V1_min; V2(V2==Inf)

　　=V2_max;

　　V2(V2>V2_max)=V2_max; V=

　　V1 + V2; [E_x, E_y]=gradient(-V); hold on; grid on; t=linspace(-pi, pi, 25); px=0、1 * cos(t) + x_p; py=0、1 * sin(t) + y_p; streamline(x, y, E_x, E_y, px, py);%画出电场线 sx=[min(d)/3*2,min(d),min(d),min(d),min(d)/3*2,min(d),min(d),min(d)/3*1,0,max(d)/3*1,max(d)/3*2]; sy=[min(d),min(d)/3*1, 0,max(d)/3*1, max(d),max(d)/3*2,max(d),max(d),max(d),max(d),max(d)]; streamline(x, y, E_x, E_y, sx, sy);%画出电场线 contour(x, y, V, linspace(min(V(:)), max(V(:)), 180));%画出等位线

　　plot(x_n, y_n, "ro",

　　x_n, y_n, "r-", "MarkerSize", 16);

　　plot(x_p, y_p, "ro",

　　x_p, y_p, "r+", "MarkerSize", 16);

　　axis([min(d), max(d), min(d), max(d)]);

　　title("电偶极子的场与等位线"); hold off; 附录二: function x=newton(fname,dfname,x0,e) if nargin<4,e=1e-4;end fname=inline("x*exp(x)-1"); dfname=inline("exp(x)+x*exp(x)"); x0=0、5; x=x0;x0=x+2*e; tic while abs(x0-x)>e

　　x0=x;x=x0-feval(fname,x0)/feval(dfname,x0); end toc 附录三: hx=17;hy=11;%设置网格 v1=ones(hy,hx);%设置二维数组 for j=1:hx%设置边界条件

　　v1(hy,j)=100*sin(pi*(2*(j-1)/(hx-1)));%假设恰好为一个周期

　　v1(1,j)=0; end

　　v1(:,1)=0; v2=v1;maxt=1;t=0;k=0;%初始化 while(maxt>0、00001) %迭代精度

　　k=k+1;%计算迭代总次数

　　maxt=0;

　　for i=2:hy-1

　　for j=2:hx-1

　　v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v2(i-1,j)+v2(i,j-1))/4;%拉普拉斯方程差分形式

　　t=abs(v2(i,j)-v1(i,j));

　　if(t>maxt) maxt=t;

　　end

　　end

　　end

　　v2(2:hy-1,hx)=v2(2:hy-1,hx-1);%右边界边界条件

　　v1=v2; end subplot(1,2,1),mesh(v2) %3D 网格图 axis([0,17,0,14,-20,100]) subplot(1,2,2),contour(v2,16)

　　hold on x=1:1:hx;y=1:1:hy; [xx,yy]=meshgrid(x,y); [Gx,Gy]=gradient(v2,0、6,0、6);%计算梯度 quiver(xx,yy,Gx,Gy,0、5,"r") %根据梯度画箭头 axis([-3、5,hx+6、5,-2,15]) plot([1,1,hx,hx,1],[1,hy,hy,1,1],"k")%画导体框 text(hx/2-2,hy+0、6,"\phi=100sin(\pix)","fontsize",11);%上标注 text(hx/2-1,0、5,"\phi=0","fontsize",11);%下标注 text(-1、8,hy/2,"\phi=0","fontsize",11);%左标注 text(hx+0、2,hy/2,"\partial\phi/\partialn=0","fontsize",11);% 右标注 title("静电场点位分布图 "); hold off 附录四: N=36;M=32;M1=12;M2=20;D=1;H=0、5;T=0、05;time=10;%初始参数定义 u=zeros(M+1,N+1);%定义场矩阵 u(M1+2:M2,1)=ones(M2-M1-1,1);%边界条件 for i=2:M for j=2:N u(i,j)=(i-1)*H*(j-1)*H;%初始条件 end end u2=u;%差分方程运算开始 for k=1:time/T%k为时间步数

　　for i=2:M

　　for j=2:N

　　u2(i,j)=(1-4*D*T/H/H)*u(i,j)+D*T/H/H*(u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1));

　　end

　　end

　　for i=1:M+1

　　for j=1:N+1

　　u(i,j)=u2(i,j);

　　end

　　u(i,N+1)=u2(i,N);

　　end

　　for i=1:M1+1

　　u(i,1)=u2(i,2);

　　end

　　for i=M2+1:M+1

　　u(i,1)=u2(i,2);

　　end end%差分方程运算结束 mesh(u)%画图 xlabel("X Axis"),ylabel("Y Axis"),zlabel("Temperature"),title("Thermal Field Distribution") 附录五: a=10; N=100; n1=sqrt(N); ltt=ones(N,N); b=a/n1; e1=1e-9; E=1/36/pi*e1; %介电常数 for i=1:n1

　　%获取各小块中心坐标

　　for j=1:n1

　　k=n1*(i-1)+j;

　　x(k)=(2*i-1)*b;

　　y(k)=(2*j-1)*b;

　　end end for m=1:N

　　for n=1:N

　　if m==n

　　ltt(m,n)=2*b/pi/E*0、8814;

　　else

　　ltt(m,n)=b^2/pi/E/sqrt((x(m)-x(n))^2+(y(m)-y(n))^2);

　　end

　　end end L1=ltt; L2=inv(L1); Lsum=sum(sum(L2)); C=4*b^2*Lsum

•  下载文档
•  收藏
•  分享 赏
•  0