矩量法实验报告.doc

1、矩量法实验报告姓名：学号：导师：班级：年月日题目一：用矩量法计算，边界条件为分析：显然，这是一个简单的边值问题，其精确解为（1）下面用矩量法求解这个问题，我们选择基函数为（2）则，原微分方程的解可以写成级数展开式为（3）对于检验函数我们选择（4）在这种情况下，就是伽略金法。由內积公式，（5）得，（6） （7）同时，由 （8）式中L是线性算子，g为已知函数，为未知函数。令在L的定义域中被展开为的组合，如（9）式中是系数。由于算子L是线性的，所以有（10）我们已经规定了一个合适的内积，由（6）、（7）式可把上式写成矩阵形式为（11）由此可求得（12）最后再把上式代入（3）式，即可得矩量法结果。因为

2、这是一个简单的微分方程，有精确解，所以为了体现N取不同值的时候矩量法的逼近程度，所以取N从13时矩量法的计算结果，并和解析解做比较。N=1时，由式（12）得。N=2时，得N=3时，得显然第三级解，即N=3时，矩量法所得的解和解析解是完全相同的。为了便于比较，把N取不同值的曲线画在同一张图里面，如图1。由图可以看出，当N=3的时候，用矩量法所得的解和解析解是完全相同的。源程序代码：clearclcx=linspace(0,1,100); %先画出解析结果以便和矩量法的结果相比较f0=5/6.*x-1/2.*x.2-1/3.*x.4;plot(x,f0,gp);grid onaxis(0 1 0

3、0.3)title(矩量法计算二次微分函数);hold on;for N=1:3 %N从1到3分别取不同的值，在此用循环分别计算之，更方便 f=0; l=zeros(N,N);g=zeros(N,1); %由于每次循环所用到的矩阵l、g的维数是不一样的，所以每次内循环之前都要先对矩阵初始化，这样可以加快运算的速度 for m=1:N g(m)=m*(3*m+8)/(2*(m+2)*(m+4); %与矩量法相应的激励向量 for n=1:N l(m,n)=m*n/(m+n+1); %与矩量法相应的阻抗矩阵 end end alpha=lg; %计算出每次的alpha for n=1:N %在上面

4、计算出一次alpha的值的时候，即马上画出相应的曲线 f=f+alpha(n)*(x-x.(n+1); end draw(x,f,N); %自定义的画图函数，在N取不同值时，赋予画图函数不同的线形 hold onendlegend(解析解,N取1时,N取2时,N取3时); %由画出的图可以看出，在N=3的时候，用矩量法实现的曲线与解析解画出的曲线完全重合function draw(x,y,n) %画图函数if n=1 plot(x,y,r);elseif n=2 plot(x,y,b);elseif n=3 plot(x,y,k);elseif n=4 plot(x,y,co)end 题目二：

5、Z一块正方形的导体板，边长为2a米，位于z=0的平面上，中心在坐标原点，如图2所示。设表示道题板上的面电荷密度，板的厚度为零。XYO图二则空间任一点的静电位是（1）式中，板上的边界条件是（常数）。此时方程是（在板上 z=0）（2）式中，待求的未知函数是电荷密度。一个有意义的参数是导体板的电容（3）假设将导体板划分为N个正方形小块。定义函数（4）则电荷密度就可表示为（5）将（5）代入（2），并且在每个的中点满足所得的方程，则有（6）式中（7）注意，是上单位振幅的均匀电荷密度在的中心处产生的电位。由求解式（6）得到。据此，电荷密度由式（5）逼近，对于式（3）的平板电容相应地近似为（8）此结果可以解

6、释为：物体的电容是其各小块电容的总和加上每一对小块间的互电容。为了将上述结果翻译成线性空间和矩量法的语言。令（9）（10）（板上电位）（11）于是与（2）式等效。取内积（12）为了应用矩量法，以函数式（4）为分域基并规定检验函数为（13）这是一个二维的狄拉克函数。为了得到数值结果，必须计算（7）式的。令表示每个的边长，由本身面上的单位电荷密度在其中心处产生的电位是（14）上单位电荷在中心处产生的电位可用同样的方法计算，但算式复杂。若将上的电荷视为点电荷，并应用（15）值得注意的是，在本题的编程和计算中要特别注意导体板的边长和2a和分块的边长2b的关系，同时注意a的赋值，a并不是边长，2a才是导

7、体板的边长。本题计算时，取分块数N=100。这是得到导体板的电容值为同时，沿导体板的电荷密度的分布的二维和三维图如下，图三和图四。图三图四程序代码如下：clearclc%=% %确定初始量%=%N=100; %确定导体板的分块数a=0.5;b=a/sqrt(N); %给定正方形导体板的边长ebslong=8.854e-12; %介电常数l=zeros(N); %给定l(m,n)的阶数，这样可以缩短循环的时间gm=ones(1,N); %给gm定值，同时确定阶数%确定每个分块的中心坐标x=-a+2*a/sqrt(N)/2:2*a/sqrt(N):a-2*a/sqrt(N)/2;%建立坐标系，以正

8、方形板的中心为坐标原点X=zeros(2,N);k=0; %矩阵初始化也即坐标生成,循环变量初始化for n=1:sqrt(N) for m=1:sqrt(N) k=k+1; X(1,k)=x(n); endendk=0;for n=1:sqrt(N) for m=1:sqrt(N) k=k+1; X(2,k)=x(m); endend %生成每个分块的中心坐标%求出面电荷密度for n=1:N for m=1:N if m=n l(m,n)=2*b*0.8814/(pi*ebslong); %计算m等于n的时候的元素 else l(m,n)=b2/(pi*ebslong*sqrt(X(1,m

9、)-X(1,n)2+(X(2,m)-X(2,n)2);%计算n不等于m的时候的矩阵元素 end endendalpha=lgm; %求出nC=(2*b)2*sum(alpha) %求出导体板的电容density=zeros(1,sqrt(N);m=0;for n=sqrt(N)/2:sqrt(N):N-sqrt(N)/2 m=m+1; density(m)=alpha(n); %沿导体板中间线的电荷密度end%为了画图的方便，在此画图的时候不再是以原来所建立的坐标系为参考，而是以每个分块的编号为参考figuresurface(reshape(alpha,sqrt(N),sqrt(N); %导体

10、板的三维电荷密度分布title(电荷分布三维图);xlabel(沿X轴的距离);ylabel(沿Y轴的距离);zlabel(电荷密度/电位);figureplot(density,r); grid on %电荷密度分布的横切面也即沿导体板中间线的电荷密度分布title(电荷密度);xlabel(沿板方向上的距离);ylabel(电荷密度/电位);figureplot(alpha);title(所有分块上的电荷密度);xlabel(分块的编号); %所有分块上的电荷密度的分布，可以看出在导体板的两边是有边沿效应的ylabel(电荷密度/电位); %以上图形的横坐标均没有归一化题目三：计算一个长度

11、为，直径为的直导线天线的方向图。分析：为了求解，可以用网络参数描述问题的解，把导线看成是N个小段连在一起的，每一个小段的终点确定了在空间的一对端点，这N对端点可以想像成一个N端口网络，而短路所有网络的端口就得到了线状物体。我们可以采用把电流源依次加在每个端口上，而在所有端口计算开路电压，就求得N端口的阻抗矩阵。此方法只包含空间的电流元。导纳矩阵是阻抗矩阵的逆矩阵，只要导纳矩阵已知，则对任一个特定激励的电压（外加电压），都可以用矩阵乘法来找出端口电流（在导线上的电流分布）。在已知外加场的作用下，在导体S上的电荷密度和电流密度J的方程可用下述方程求得。用滞后位的积分来表示由和J产生的散射场，并应用

12、在s上的边界条件，这些公式归纳如下： （1） （2） （3） （4）在S上 （5）图1表示一根任意的细导线，对于它可作如下的近似：（1）假定电流只是沿着导线轴的方向流动；（2）电流和电荷密度可以近似地认为是线电流I及在导线轴上的；（3）只对导线表面上E的轴向分量使用边界条件式（5）做了这些近似之后，式（1）至式（5）就变成导线轴细导线 图一图二 在S上（6）（7）（8）（9）式中是沿导线轴的长度变量，R是从轴上源点指向导线表面的场点之间的距离。想要用计算机对上述方程求解，则需对上面的方程离散化。积分可近似为沿N个小段积分的总和，此时，在每个小段上视I和q为常数。在积分所处的相同区间上，导数可由

13、有限差分来近似。图2表明将导线轴划分为N个小段。第n小段由始点，中点n和终点组成，增量表明是在和之间，和分别表示增量沿上移动负的和正的二分之一增量。所需要的式（6）至（9）的近似为（10）（11）（12） （13）同时还有类似于式（12）和式（13）的和的方程式。由式（13）可以看出，各个可以用各个I来表示。因此，式（10）可以写成只包含的形式。我们可以认为由式（10）表示的N个方程是一个带有端子对的N端口网络方程。外加到每个端口的电压近似为。因此，定义矩阵 （14）我们便可以将（10）式写成矩阵形式如下： （15）将式（11）至（13）代入式（10）并重新排成式（15）的形式，便可以得到矩阵

14、Z的元素。另一方面，可将式（10）至（13）用于两个孤立元素而直接得出阻抗元素。我们要用的是后一种方式，因为它比较容易做到。现在研究如图3所示的两个代表性的导线散射体元。式（11）和式（12）具有相同的积分形式，即可以表示为（16）图三 导线的两小段式中是从上一点到m点的距离。符号+与在适当的时候加在m与n上。令图3的元素n由一个线电流I(n)和静电荷为（17）的两个线电荷组成。式中。由式（11），在m点由I(n)产生的矢位为（18）在和点由式（17）的电荷产生的标位为（19）将式（18）和式（19）代入式（10），且形成，则可得到（20）此结果适用于互阻抗，也同样适用于自阻抗(m=n)。若两

15、电流元相隔很远时，则可使用一个比较简单的公式，此公式可以根据电流元产生的辐射场得出。在所含的近似条件下，可以用其阻抗矩阵完全表现线状物体的特性。物体由在导线轴上的2N个点加上导线半径来确定。阻抗元素由式（20）计算，而电压矩阵则由式（14）决定的外加场来求得。在散射体N个点上的电流则由电流矩阵从式（15）的逆矩阵求得I=YV （21）只要电流分布已知，则各种有用参数，如场方向图、输入阻抗、散射面积等都可以用相应公式以数值计算方法算出来。按照矩量法，以上的解相当于利用脉冲函数同时作电流和电荷的展开函数，而以点选配作为检验函数。为了避免微分，将这个方法用于一个有限差分代替微分而得的近似算子。应当指

16、出，导线的终点可以看成带有零电流的一小段的中点，从线端留出一个小段的一半然后才开始分段，它在数学上等于在线端得边界条件I=0。注意在线端得电荷并不是零，这与延伸出以外半个区段以代表电荷是一致的。只要计算出来，阻抗元素式（20）就是已知的。为此，我们建立一个局部坐标系，使原点在n上，z轴沿着，如图4所示，则（22）式中 （23）是线的半径。将指数展开为马格劳林级数，得到的近似式：（24）式中第一项相当于线电荷的静电电位，第二项与无关。当m=n时，这两项给出的精确度是令人满意的，而XYZ 图四 （25）若时，粗略的近似是将在积分式（22）中作为常数，则 （26）式中是从n到m的距离。因为式（26）

17、在极限时是不精确的，正如在原来的讨论中所讨论过的那样，它有一个残留误差。一根导线在其沿线一点上或者更多点上加以一个集总电压源来激励，就是一个线天线。如果导线在第i个区间被激励，则外加的电压矩阵式变为 （27）这就是说，除了等于电源电压的第i项之外，其余各项为零。电流分布由式（21）变为 （28）因此，导纳矩阵的第i列就是在第i区间加上单位电压源时的电流分布。这样，阻抗矩阵的求逆运算同时给出了沿任意点激励的天线电流分布。导纳矩阵中的对角线元素是导线在第i区间馈电的输入导纳，是在第i区间的端口与第j区间的端口之间的转移导纳。辐射方向图的推导可以根据互易定理获得。图5表明一个远区电流元，调整到使其在

18、天线附近产生一个单位振幅的平面波： （29）时所需之值。式中是规定波的极化方向的单位矢量，是指向波的传播方向的波数矢量，而是指向天线上n点的矢径。根据互易定理 （30）式中是天线产生的E的分量，I是天线电流，常数是为在原点产生单位振幅的平面波所必需之值，即（31）式（30）的数值近似式可以由定义一个电压矩阵而得到 （32）其中由式（29）给出，并将式（30）近似为矩阵相乘（33）其中是V的转置矩阵。注意，是与式（14）同样类型的矩阵，即为导线上平面波激励的电压矩阵。式（33）对于一个任意激励保持有效。辐射场分量的功率方向图为 （34）式中为空间波阻抗，是天线的输入功率：（35）对于单个源的特殊

19、情况，即方程式（27）的情况，变为简单的。将式（33）和式（35）代入式（34），则得（36）其中是由式（32）在不同的入射角和下得出的。式（36）给出了只有单一极化辐射场时的增益方向图。如果要求总功率增益方向图时，正交极化的各个g值应加在一起。图五源程序代码：clearlamda=1;%波长ra=0.005*lamda;%振子的半径me=8.85e-12;%介电常数mu=4*pi*(1e-7);%磁导率c=3e+8;f=c/lamda;arg=2*pi*f;%角频率tl=0.5*lamda;%振子的总长nm=21;%匹配点数目pi=3.14159265;rad=pi/180;beta=2.0

20、*pi;eta=120*pi;hl=tl/2;nmh=0.5*(nm+1);dz=2*hl/nm;zm=hl-0.5*dz;b=0.5*dz;a=-0.5*dz;n=79;hd=(b-a)/(n+1);lzm=-hl+dz/2:dz:hl-dz/2;%匹配点for I = 1: nmzn=hl-(I-0.5)*dz;za1=zn-zm+a;recgp=sqrt(ra*ra+za1*za1);cgp1=exp(-1i*beta*recgp)*(1.0+1i*beta*recgp)*(2.0*recgp*recgp-3.0*ra*ra)+(beta*ra*recgp)2)/(2.0*beta*re

21、cgp5);zb1=zn-zm+b;roc=sqrt(ra*ra+zb1*zb1);cgp2=exp(-1i*beta*roc)*(1.0+1i*beta*roc)*(2.0*roc*roc-3.0*ra*ra)+(beta*ra*roc)2)/(2.0*beta*roc5);crt=cgp1+cgp2;for k = 1: nxk=a+k*hd;zx1=zn-zm+xk;r=sqrt(ra*ra+zx1*zx1);cgp3=exp(-1i*beta*r)*(1.0+1i*beta*r)*(2.0*r*r-3.0*ra*ra)+(beta*ra*r)2)/(2.0*beta*r5);if mo

22、d(k,2)=0crt=crt+4.0*cgp3;elsecrt=crt+2.0*cgp3;endendcrt=crt*hd*0.33333;zmn(I)=crt;if I=1zmn(nm+I-1)=crt;endendfor n=1:nmp(1,n)=zmn(n);endfor m=1:nmp(m,1)=zmn(m);endfor m=2:nmfor n=2:nmp(m,n)=p(m-1,n-1);endendV=zeros(nm,1);fedp=(nm+1)/2;%馈电点的位置V(fedp)=-1i*beta/(120*pi*dz);%馈电的电位I=inv(p)*V;Z=1/I(fedp)%输入阻抗figuresubplot(2,1,1);plot(lzm,abs(I),grid on;xlabel(L/lamda);ylabel(电流幅值);title(电流分布);subplot(2,1,2);plot(lzm,180*angle(I)/pi),grid on;xlabel(L/lamda);ylabel(电流相位);