矩量法实验报告.doc
《矩量法实验报告.doc》由会员分享，可在线阅读，更多相关《矩量法实验报告.doc（21页珍藏版）》请在冰点文库上搜索。
矩量法实验报告
姓名：
学号：
导师：
班级：
年月日
题目一：
用矩量法计算，边界条件为
分析：
显然，这是一个简单的边值问题，其精确解为
（1）
下面用矩量法求解这个问题，我们选择基函数为
（2）
则，原微分方程的解可以写成级数展开式为
（3）
对于检验函数我们选择
（4）
在这种情况下，就是伽略金法。
由內积公式，
（5）
得，
（6）
（7）
同时，由
（8）
式中L是线性算子，g为已知函数，为未知函数。
令在L的定义域中被展开为的组合，如
（9）
式中是系数。
由于算子L是线性的，所以有
（10）
我们已经规定了一个合适的内积，由（6）、（7）式可把上式写成矩阵形式为
（11）
由此可求得
（12）
最后再把上式代入（3）式，即可得矩量法结果。
因为这是一个简单的微分方程，有精确解，所以为了体现N取不同值的时候矩量法的逼近程度，所以取N从1~3时矩量法的计算结果，并和解析解做比较。
N=1时，，由式（12）得。
N=2时，得 N=3时，得
显然第三级解，即N=3时，矩量法所得的解和解析解是完全相同的。
为了便于比较，把N取不同值的曲线画在同一张图里面，如图1。
由图可以看出，当N=3的时候，用矩量法所得的解和解析解是完全相同的。
源程序代码：
clear
clc
x=linspace（0,1,100）;%先画出解析结果以便和矩量法的结果相比较
f0=5/6.*x-1/2.*x.^2-1/3.*x.^4;
plot（x,f0,'gp'）;
gridon
axis（[0100.3]）
title（'矩量法计算二次微分函数'）;
holdon;
forN=1:
3%N从1到3分别取不同的值，在此用循环分别计算之，更方便
f=0;
l=zeros（N,N）;g=zeros（N,1）;
%%由于每次循环所用到的矩阵l、g的维数是不一样的，所以每次内循环之前都要先对矩阵初始化，这样可以加快运算的速度
form=1:
N
g（m）=m*（3*m+8）/（2*（m+2）*（m+4））;%与矩量法相应的激励向量
forn=1:
N
l（m,n）=m*n/（m+n+1）;%与矩量法相应的阻抗矩阵
end
end
alpha=l\g;%计算出每次的alpha
forn=1:
N%在上面计算出一次alpha的值的时候，即马上画出相应的曲线
f=f+alpha（n）*（x-x.^（n+1））;
end
draw（x,f,N）;%自定义的画图函数，在N取不同值时，赋予画图函数不同的线形
holdon
end
legend（'解析解','N取1时','N取2时','N取3时'）;%由画出的图可以看出，在N=3的时候，用矩量法实现的曲线与解析解画出的曲线完全重合
functiondraw（x,y,n）%画图函数
ifn==1
plot（x,y,'r'）;
elseifn==2
plot（x,y,'b'）;
elseifn==3
plot（x,y,'k'）;
elseifn==4
plot（x,y,'co'）
end
题目二：
Z
一块正方形的导体板，边长为2a米，位于z=0的平面上，中心在坐标原点，如图2所示。
设表示道题板上的面电荷密度，板的厚度为零。
X
Y
O
图二
则空间任一点的静电位是
（1）
式中，板上的边界条件是（常数）。
此时方程是
（在板上z=0）
（2）
式中，待求的未知函数是电荷密度。
一个有意义的参数是导体板的电容
（3）
假设将导体板划分为N个正方形小块。
定义函数
（4）
则电荷密度就可表示为
（5）
将（5）代入
（2），并且在每个的中点满足所得的方程，则有
（6）
式中 （7）
注意，是上单位振幅的均匀电荷密度在的中心处产生的电位。
由求解式（6）得到。
据此，电荷密度由式（5）逼近，对于式（3）的平板电容相应地近似为
（8）
此结果可以解释为：
物体的电容是其各小块电容的总和加上每一对小块间的互电容。
为了将上述结果翻译成线性空间和矩量法的语言。
令
（9）
（10）
（板上电位） （11）
于是与
（2）式等效。
取内积
（12）
为了应用矩量法，以函数式（4）为分域基并规定检验函数为
（13）
这是一个二维的狄拉克函数。
为了得到数值结果，必须计算（7）式的。
令表示每个的边长，由本身面上的单位电荷密度在其中心处产生的电位是
（14）
上单位电荷在中心处产生的电位可用同样的方法计算，但算式复杂。
若将上的电荷视为点电荷，并应用
（15）
值得注意的是，在本题的编程和计算中要特别注意导体板的边长和2a和分块的边长2b的关系，同时注意a的赋值，a并不是边长，2a才是导体板的边长。
本题计算时，取分块数N=100。
这是得到导体板的电容值为
同时，沿导体板的电荷密度的分布的二维和三维图如下，图三和图四。
图三
图四
程序代码如下：
clear
clc
%===========================================%
%确定初始量
%===========================================%
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;%建立坐标系，以正方形板的中心为坐标原点
X=zeros（2,N）;
k=0;%矩阵初始化也即坐标生成,循环变量初始化
forn=1:
sqrt（N）
form=1:
sqrt（N）
k=k+1;
X（1,k）=x（n）;
end
end
k=0;
forn=1:
sqrt（N）
form=1:
sqrt（N）
k=k+1;
X（2,k）=x（m）;
end
end%生成每个分块的中心坐标
%%求出面电荷密度
forn=1:
N
form=1:
N
ifm==n
l（m,n）=2*b*0.8814/（pi*ebslong）;%计算m等于n的时候的元素
else
l（m,n）=b^2/（pi*ebslong*sqrt（（X（1,m）-X（1,n））^2+（X（2,m）-X（2,n））^2））;%计算n不等于m的时候的矩阵元素
end
end
end
alpha=l\gm';%求出αn
C=（2*b）^2*sum（alpha）%求出导体板的电容
density=zeros（1,sqrt（N））;
m=0;
forn=sqrt（N）/2:
sqrt（N）:
N-sqrt（N）/2
m=m+1;
density（m）=alpha（n）;%沿导体板中间线的电荷密度
end
%%为了画图的方便，在此画图的时候不再是以原来所建立的坐标系为参考，而是以每个分块的编号为参考
figure
surface（reshape（alpha,sqrt（N）,sqrt（N）））;%导体板的三维电荷密度分布
title（'电荷分布三维图'）;xlabel（'沿X轴的距离'）;
ylabel（'沿Y轴的距离'）;zlabel（'电荷密度/电位'）;
figure
plot（density,'r'）;gridon%电荷密度分布的横切面也即沿导体板中间线的电荷密度分布
title（'电荷密度'）;xlabel（'沿板方向上的距离'）;
ylabel（'电荷密度/电位'）;
figure
plot（alpha）;
title（'所有分块上的电荷密度'）;xlabel（'分块的编号'）;%所有分块上的电荷密度的分布，可以看出在导体板的两边是有边沿效应的
ylabel（'电荷密度/电位'）;%以上图形的横坐标均没有归一化
题目三：
计算一个长度为，直径为的直导线天线的方向图。
分析：
为了求解，可以用网络参数描述问题的解，把导线看成是N个小段连在一起的，每一个小段的终点确定了在空间的一对端点，这N对端点可以想像成一个N端口网络，而短路所有网络的端口就得到了线状物体。
我们可以采用把电流源依次加在每个端口上，而在所有端口计算开路电压，就求得N端口的阻抗矩阵。
此方法只包含空间的电流元。
导纳矩阵是阻抗矩阵的逆矩阵，只要导纳矩阵已知，则对任一个特定激励的电压（外加电压），都可以用矩阵乘法来找出端口电流（在导线上的电流分布）。
在已知外加场的作用下，在导体S上的电荷密度和电流密度J的方程可用下述方程求得。
用滞后位的积分来表示由和J产生的散射场，并应用在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为常数。
在积分所处的相同区间上，导数可由有限差分来近似。
图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）的形式，便可以得到矩阵[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）。
若两电流元相隔很远时，则可使用一个比较简单的公式，此公式可以根据电流元产生的辐射场得出。
在所含的近似条件下，可以用其阻抗矩阵完全表现线状物体的特性。
物体由在导线轴上的2N个点加上导线半径来确定。
阻抗元素由式（20）计算，而电压矩阵则由式（14）决定的外加场来求得。
在散射体N个点上的电流则由电流矩阵从式（15）的逆矩阵求得
[I]=[Y][V] （21）
只要电流分布已知，则各种有用参数，如场方向图、输入阻抗、散射面积等都可以用相应公式以数值计算方法算出来。
按照矩量法，以上的解相当于利用脉冲函数同时作电流和电荷的展开函数，而以点选配作为检验函数。
为了避免微分，将这个方法用于一个有限差分代替微分而得的近似算子。
应当指出，导线的终点可以看成带有零电流的一小段的中点，从线端留出一个小段的一半然后才开始分段，它在数学上等于在线端得边界条件I=0。
注意在线端得电荷并不是零，这与延伸出以外半个区段以代表电荷是一致的。
只要计算出来，阻抗元素式（20）就是已知的。
为此，我们建立一个局部坐标系，使原点在n上，z轴沿着，如图4所示，则
（22）
式中
（23）
是线的半径。
将指数展开为马格劳林级数，得到的近似式：
（24）
式中第一项相当于线电荷的静电电位，第二项与无关。
当m=n时，这两项给出的精确度是令人满意的，而
X
Y
Z
图四
（25）
若时，粗略的近似是将在积分式（22）中作为常数，则
（26）
式中是从n到m的距离。
因为式（26）在极限时是不精确的，正如在原来的讨论中所讨论过的那样，它有一个残留误差。
一根导线在其沿线一点上或者更多点上加以一个集总电压源来激励，就是一个线天线。
如果导线在第i个区间被激励，则外加的电压矩阵式变为
（27）
这就是说，除了等于电源电压的第i项之外，其余各项为零。
电流分布由式（21）变为
（28）
因此，导纳矩阵的第i列就是在第i区间加上单位电压源时的电流分布。
这样，阻抗矩阵的求逆运算同时给出了沿任意点激励的天线电流分布。
导纳矩阵中的对角线元素是导线在第i区间馈电的输入导纳，是在第i区间的端口与第j区间的端口之间的转移导纳。
辐射方向图的推导可以根据互易定理获得。
图5表明一个远区电流元，调整到使其在天线附近产生一个单位振幅的平面波：
（29）
时所需之值。
式中是规定波的极化方向的单位矢量，是指向波的传播方向的波数矢量，而是指向天线上n点的矢径。
根据互易定理
（30）
式中是天线产生的E的分量，I是天线电流，常数是为在原点产生单位振幅的平面波所必需之值，即
（31）
式（30）的数值近似式可以由定义一个电压矩阵而得到
（32）
其中由式（29）给出，并将式（30）近似为矩阵相乘
（33）
其中是[V]的转置矩阵。
注意，是与式（14）同样类型的矩阵，即为导线上平面波激励的电压矩阵。
式（33）对于一个任意激励保持有效。
辐射场分量的功率方向图为
（34）
式中为空间波阻抗，是天线的输入功率：
（35）
对于单个源的特殊情况，即方程式（27）的情况，变为简单的。
将式（33）和式（35）代入式（34），则得
（36）
其中是由式（32）在不同的入射角和下得出的。
式（36）给出了只有单一极化辐射场时的增益方向图。
如果要求总功率增益方向图时，正交极化的各个g值应加在一起。
图五
源程序代码：
clear
lamda=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*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;%匹配点
forI=1:
nm
zn=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*recgp^5）;
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*roc^5）;
crt=cgp1+cgp2;
fork=1:
n
xk=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*r^5）;
ifmod（k,2）~=0
crt=crt+4.0*cgp3;
else
crt=crt+2.0*cgp3;
end
end
crt=crt*hd*0.33333;
zmn（I）=crt;
ifI~=1
zmn（nm+I-1）=crt;
end
end
forn=1:
nm
p（1,n）=zmn（n）;
end
form=1:
nm
p（m,1）=zmn（m）;
end
form=2:
nm
forn=2:
nm
p（m,n）=p（m-1,n-1）;
end
end
V=zeros（nm,1）;
fedp=（nm+1）/2;%馈电点的位置
V（fedp）=-1i*beta/（120*pi*dz）;%馈电的电位
I=inv（p）*V;
Z=1/I（fedp）%输入阻抗
figure
subplot（2,1,1）;
plot（lzm,abs（I））,gridon;
xlabel（'L/lamda'）;
ylabel（'电流幅值'）;
title（'电流分布'）;
subplot（2,1,2）;
plot（lzm,180*angle（I）/pi）,gridon;
xlabel（'L/lamda'）;
ylabel（'电流相位'）;
升级会员 