矩量法实验报告

矩量法实验报告姓名：学号：导师：班级：年月日题目一：用矩量法计算22214dfxdx，边界条件为(0)(1)0ff分析：显然，这是一个简单的边值问题，其精确解为24511()623fxxxx（1）下面用矩量法求解这个问题，我们选择基函数为1,1,2,3,,nnfxxnN（2）则，原微分方程的解可以写成级数展开式为11()Nnnnfxx（3）对于检验函数我们选择1nnnfxx（4）在这种情况下，就是伽略金法。由內积公式，10,()()fgfxgxdx（5）得，1110,()(1)1mnmnmnmnlLfxxnnxdxmn（6）1120(38),()(14)2(2)(4)mmmmmggxxxdxmm（7）同时，由()Lfg（8）式中L是线性算子，g为已知函数，f为未知函数。令f在L的定义域中被展开为123,,,fff的组合，如nnnff（9）式中n是系数。由于算子L是线性的，所以有()nnnLfg（10）我们已经规定了一个合适的内积，由（6）、（7）式可把上式写成矩阵形式为mnnmlg（11）由此可求得1nmnmlg（12）最后再把上式代入（3）式，即可得矩量法结果。因为这是一个简单的微分方程，有精确解，所以为了体现N取不同值的时候矩量法的逼近程度，所以取N从1~3时矩量法的计算结果，并和解析解做比较。N=1时，111111,330lg，由式（12）得11110。N=2时，得N=3时，得121102312312013显然第三级解，即N=3时，矩量法所得的解和解析解是完全相同的。为了便于比较，把N取不同值的曲线画在同一张图里面，如图1。00.10.20.30.40.50.60.70.80.9100.050.10.150.20.25矩量法计算二次微分函数图一解析解N取1时N取2时N取3时由图可以看出，当N=3的时候，用矩量法所得的解和解析解是完全相同的。源程序代码：clearclcx=linspace(0,1,100);%先画出解析结果以便和矩量法的结果相比较f0=5/6.*x-1/2.*x.^2-1/3.*x.^4;plot(x,f0,'gp');gridonaxis([0100.3])title('矩量法计算二次微分函数');holdon;forN=1:3%N从1到3分别取不同的值，在此用循环分别计算之，更方便f=0;l=zeros(N,N);g=zeros(N,1);%%由于每次循环所用到的矩阵l、g的维数是不一样的，所以每次内循环之前都要先对矩阵初始化，这样可以加快运算的速度form=1:Ng(m)=m*(3*m+8)/(2*(m+2)*(m+4));%与矩量法相应的激励向量forn=1:Nl(m,n)=m*n/(m+n+1);%与矩量法相应的阻抗矩阵endendalpha=l\g;%计算出每次的alphaforn=1:N%在上面计算出一次alpha的值的时候，即马上画出相应的曲线f=f+alpha(n)*(x-x.^(n+1));enddraw(x,f,N);%自定义的画图函数，在N取不同值时，赋予画图函数不同的线形holdonendlegend('解析解','N取1时','N取2时','N取3时');%由画出的图可以看出，在N=3的时候，用矩量法实现的曲线与解析解画出的曲线完全重合functiondraw(x,y,n)%画图函数ifn==1plot(x,y,'r');elseifn==2plot(x,y,'b');elseifn==3plot(x,y,'k');elseifn==4plot(x,y,'co')end题目二：一块正方形的导体板，边长为2a米，位于z=0的平面上，中心在坐标原点，如图2所示。设(,)xy表示道题板上的面电荷密度，板的厚度为零。图二则空间任一点的静电位是XYO(,)xyZ1(,)(',')(,,)''44aaSaaxyxyxyzdSdxdyRR（1）式中222(')(')Rxxyyz，板上的边界条件是V（常数）。此时方程是22(',')''4(')(')aaaaxyVdxdyxxyy（在板上z=0）（2）式中,xaya，待求的未知函数是电荷密度(,)xy。一个有意义的参数是导体板的电容11(,)aaSaaqCdSdxdyxyVVV（3）假设将导体板划分为N个正方形小块。定义函数10nnmfs在s上在所有其他上（4）则电荷密度就可表示为1(,)Nnnnxyf（5）将（5）代入（2），并且在每个ms的中点(,)mmxy满足所得的方程，则有1,1,2,,NmnnnVlmN（6）式中221''4(')(')nnmnxyldxdyxxyy（7）注意，mnl是ns上单位振幅的均匀电荷密度在ms的中心处产生的电位。由求解式（6）得到n。据此，电荷密度由式（5）逼近，对于式（3）的平板电容相应地近似为111NnnmnnnmnCslsV（8）此结果可以解释为：物体的电容是其各小块电容的总和加上每一对小块间的互电容。为了将上述结果翻译成线性空间和矩量法的语言。令(,)(,)fxyxy（9）(,),,gxyVxaya（10）22(',')()''4(')(')aaaafxyLfdxdyxxyy（板上电位）（11）于是()Lfg与（2）式等效。取内积,(,)(,)aaaafgdxdyfxygxy（12）为了应用矩量法，以函数式（4）为分域基并规定检验函数为()()mmmxxyy（13）这是一个二维的狄拉克函数。为了得到数值结果，必须计算（7）式的mnl。令22abN表示每个ns的边长，由ns本身面上的单位电荷密度在其中心处产生的电位是22122ln(12)(0.8814)4bbmnbbbbldxdyxy（14）ns上单位电荷在ms中心处产生的电位可用同样的方法计算，但算式复杂。若将ns上的电荷视为点电荷，并应用222,4()()nmnmnmnmnsblmnRxxyy（15）值得注意的是，在本题的编程和计算中要特别注意导体板的边长和2a和分块的边长2b的关系，同时注意a的赋值，a并不是边长，2a才是导体板的边长。本题计算时，取分块数N=100。这是得到导体板的电容值为39.832/CFm同时，沿导体板的电荷密度的分布的二维和三维图如下，图三和图四。1234567891023456x10-11电荷密度沿板方向上的距离电荷密度/电位图三02468100510234567891011x10-11沿X轴的距离电荷分布三维图沿Y轴的距离电荷密度/电位图四程序代码如下：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;%建立坐标系，以正方形板的中心为坐标原点X=zeros(2,N);k=0;%矩阵初始化也即坐标生成,循环变量初始化forn=1:sqrt(N)form=1:sqrt(N)k=k+1;X(1,k)=x(n);endendk=0;forn=1:sqrt(N)form=1:sqrt(N)k=k+1;X(2,k)=x(m);endend%生成每个分块的中心坐标%%求出面电荷密度forn=1:Nform=1:Nifm==nl(m,n)=2*b*0.8814/(pi*ebslong);%计算m等于n的时候的元素elsel(m,n)=b^2/(pi*ebslong*sqrt((X(1,m)-X(1,n))^2+(X(2,m)-X(2,n))^2));%计算n不等于m的时候的矩阵元素endendendalpha=l\gm';%求出αnC=(2*b)^2*sum(alpha)%求出导体板的电容density=zeros(1,sqrt(N));m=0;forn=sqrt(N)/2:sqrt(N):N-sqrt(N)/2m=m+1;density(m)=alpha(n);%沿导体板中间线的电荷密度end%%为了画图的方便，在此画图的时候不再是以原来所建立的坐标系为参考，而是以每个分块的编号为参考figuresurface(reshape(alpha,sqrt(N),sqrt(N)));%导体板的三维电荷密度分布title('电荷分布三维图');xlabel('沿X轴的距离');ylabel('沿Y轴的距离');zlabel('电荷密度/电位');figureplot(density,'r');gridon%电荷密度分布的横切面也即沿导体板中间线的电荷密度分布title('电荷密度');xlabel('沿板方向上的距离');ylabel('电荷密度/电位');figureplot(alpha);title('所有分块上的电荷密度');xlabel('分块的编号');%所有分块上的电荷密度的分布，可以看出在导体板的两边是有边沿效应的ylabel('电荷密度/电位');%以上图形的横坐标均没有归一化题目三：计算一个长度为274.210l，直径为21a的直导线天线的方向图。分析：为了求解，可以用网络参数描述问题的解，把导线看成是N个小段连在一起的，每一个小段的终点确定了在空间的一对端点，这N对端点可以想像成一个N端口网络，而短路所有网络的端口就得到了线状物体。我们可以采用把电流源依次加在每个端口上，而在所有端口计算开路电压，就求得N端口的阻抗矩阵。此方法只包含空间的电流元。导纳矩阵是阻抗矩阵的逆矩阵，只要导纳矩阵已知，则对任一个特定激励的电压（外加电压），都可以用矩阵乘法来找出端口电流（在导线上的电流分布）。在已知外加场iE的作用下，在导体S上的电荷密度和电流密度J的方程可用下述方程求得。用滞后位的积分来表示由和J产生的散射场sE，并应用在s上的边界条件()0isnEE，这些公式归纳如下：sEjA（1）4jkRSeAJdSR（2）14jkRSedSR（3）1Jj（4）sinEnE在S上（5）图1表示一根任意的细导线，对于它可作如下的近似：（1）假定电流只是沿着导线轴的方向流动；（2）电流和电荷密度可以近似地认为是线电流I及在导线轴上的；（3）只对导线表面上E的轴向分量使用边界条件式（5）做了这些近似之后，式（1）至式（5）就变成图一图二illEjAl在S上（6）()4jkReAIldlR轴上（7）1()4jkReldlR轴上（8）1=Ijl（9）式中l是沿导线轴的长度变量，R是从轴上源点指向导线表面的场点之间的距离。想要用计算机对上述方程求解，则需对上面的方程离散化。积分可近似为沿N个小段积分的总和，此时，在每个小段上视I和q为常数。在积分所处的相同区间上，导数可由有限差分来近似。图2表明将导线轴划分为N个小段。第n小段由始点n，中点n导线轴细导线和终点n组成，增量nl表明是在n和n之间，nl和nl分别表示