电磁场数值分析期末

1《电磁场数值分析》（期末作业）---2019学年---学院：学号：姓名：联系方式：任课教师：2019年5月2作业1模拟真空中二维TM电磁波的传播，边界设置为一阶Mur吸收边界，观察电磁波的传播过程。波源为正弦函数：sin()sin(2)25ztcEtnt代码：clcclearcloseallxmesh=150;ymesh=150;mu0=4*pi*1.0E-7;eps0=8.85E-12;C=3.0E8;dx=1.0;dt=0.7*dx/C;timestep=150;ez(1:xmesh+1,1:ymesh+1)=0.0;hx(1:xmesh+1,1:ymesh)=0.0;hy(1:xmesh,1:ymesh+1)=0.0;coef1=dt/(mu0*dx);coef2=dt/(eps0*dx);coef3=(C*dt-dx)/(C*dt+dx);ez1=ez;fornow=1:timestephx=hx-coef1*(ez(:,2:ymesh+1)-ez(:,1:ymesh));hy=hy+coef1*(ez(2:xmesh+1,:)-ez(1:xmesh,:));ez(2:xmesh,2:ymesh)=ez(2:xmesh,2:ymesh)-...coef2*(hx(2:xmesh,2:ymesh)-hx(2:xmesh,1:150波源150xy3ymesh-1))+...coef2*(hy(2:xmesh,2:ymesh)-hy(1:xmesh-1,2:ymesh));ez(1,:)=ez1(2,:)+coef3*(ez(2,:)-ez1(1,:));ez(xmesh+1,:)=ez1(xmesh,:)+coef3*(ez(xmesh,:)-ez1(xmesh+1,:));ez(:,1)=ez1(:,2)+coef3*(ez(:,2)-ez1(:,1));ez(:,ymesh+1)=ez1(:,ymesh)+coef3*(ez(:,ymesh)-ez1(:,ymesh+1));ez(xmesh/2+1,ymesh/2+1)=sin(now*dt*2*pi*C/25.0);mesh(ez);pause(0.05)ez1=ez;end结果与分析：第10时间步第100时间步4第150时间步5作业2基于Pocklington方程用MoM分析半波对称振子天线：观察天线线径和分段数目分别取不同值对天线阻抗和辐射特性的影响（半径分别取0.001λ，0.0001λ，0.00001λ，分段数取11,21,31，可列表说明）代码：clearall;closeall;clc;%初始化参数c=3e8;%光速r=1%波长f=c/r;%频率w=2*pi*f;%角频率e0=8.85e-12;%介电常数u0=4*pi*1e-7;%磁导率a=0.00001*r;%半径L=0.5*r;%振子长度k=2*pi/r;%波数N=11;%分段数（奇数段）dl=L/(N+1);%每段长度（分母中+1为两头半段之和）l=L/2-dl/2;%两头空出半段，满足电流为0的边界条件lz=-l:dl:l;lzs=lz(1:N);%每一小段的起点坐标lzm=lz(1:N)+dl/2;%每一小段的中点坐标lze=lz(2:N+1);%每一小段的终点坐标6%阻抗矩阵元素求解fi=log(dl/a)/(2*pi*dl)-k/(4*pi)*1i;fi_1=exp(-k*dl*1i)/(4*pi*dl);fi_2=exp(-k*2*dl*1i)/(8*pi*dl);z=ones(N,N);form=1:Nforn=1:Nifm==nfi1=fi;fi2=fi_1;fi3=fi_1;z(m,n)=((k^2*dl^2-2)*fi1+fi2+fi3);elseifabs(m-n)==1fi1=fi_1;fi2=fi;fi3=fi_2;z(m,n)=((k^2*dl^2-2)*fi1+fi2+fi3);elsefi1=exp(-k*abs(m-n)*dl*1i)/(4*pi*abs(m-n)*dl);fi2=exp(-k*abs(m+1-n)*dl*1i)/(4*pi*abs(m+1-n)*dl);fi3=exp(-k*abs(n+1-m)*dl*1i)/(4*pi*abs(n+1-m)*dl);z(m,n)=((k^2*dl^2-2)*fi1+fi2+fi3);endendend%电压矩阵求解V=zeros(N,1);V((N+1)/2)=-1*(1i*w*e0);%计算电流系数矩阵I=z\V;%计算输入阻抗Z_in=1/I((N+1)/2);7disp(['输入阻抗=',num2str(Z_in)]);%计算振子上归一化电流分布I_amp=abs(I);Max=max(I_amp);Iunit2=[0;I_amp/Max(1);0];%两端零电流figure(1);h=0:dl/r:L/r;Ithe=sin(pi*h*r/L);%半波振子电流解析值plot(h,Iunit2,'b',h,Ithe,'r','linewidth',2);legend('pocklinton','解析值');gridon;xlabel('电长度L/\lambda');ylabel('归一化电流');%画方向图theta=0:0.01:2*pi;abs_f=zeros(1,length(theta));forn=1:1:Nabs_f=abs_f+I(n)*exp(k*(n*dl-L/2)*cos(theta)*1i);endabs_f=abs(sin(theta)*dl.*abs_f);Max_f=abs(sum(I)*dl);Far_patten2=abs_f/Max_f(1);theta_2=0:0.1:2*pi;Far_theory=abs((cos(k*(L/2)*cos(theta_2))-cos(k*L/2))./sin(theta_2));figure(2);polar(theta,Far_patten2,'-b');holdon;polar(theta_2,Far_theory,'or');holdoff;legend('pocklinton','解析值');title('半波振子天线E面方向图');8figure(3);polar(theta,ones(1,length(theta)),'-b');title('半波振子天线H面方向图');%半波振子增益I_in=I((N+1)/2);A=(w*u0)^2/(4*pi*sqrt(u0/e0)*real(Z_in)*(abs(I_in))^2);G_theta=A*abs_f.^2;Max_gain=max(G_theta)Max_gain_dB=10*log10(Max_gain);disp(['半波振子增益=',sprintf('%.4fdBi',Max_gain_dB)]);结果与分析：分段数波长1121310.001输入阻抗=81.0029+36.1494iMax_gain=1.6446半波振子增益=2.1607dBi输入阻抗=83.0507+39.2545iMax_gain=1.6503半波振子增益=2.1757dBi输入阻抗=83.9979+41.066iMax_gain=1.6520半波振子增益=2.1802dBi0.0001输入阻抗=77.7374+38.0429iMax_gain=1.6420半波振子增益=2.1538dBi输入阻抗=78.8657+39.7265iMax_gain=1.6469半波振子增益=2.1666dBi输入阻抗=79.3079+40.9582iMax_gain=1.6482半波振子增益=2.1700dBi0.00001输入阻抗=76.3172+39.4969iMax_gain=1.6408半波振子增益=2.1506dBi输入阻抗=77.1417+40.0554iMax_gain=1.6454半波振子增益=2.1626dBi输入阻抗=77.44+40.936iMax_gain=1.6465半波振子增益=2.1657dBi9作业3基于电场积分方程用MoM分析对称振子天线：计算振子总长度分别为0.25λ，0.5λ，λ，1.5λ时，振子的输入阻抗和E面方向图。代码：clearall;closeall;clc;%初始化参数c=3e8;%光速r=1;%波长f=c/r;%频率w=2*pi*f;%角频率e0=8.85e-12;%介电常数u0=4*pi*1e-7;%磁导率a=0.001*r;%半径L=1.5*r;%振子长度k=2*pi/r;%波数N=11;%分段数（奇数段）dl=L/(N+1);%每段长度（分母中+1为两头半段之和）l=L/2-dl/2;%两头空出半段，满足电流为0的边界条件lz=-l:dl:l;lzs=lz(1:N);%每一小段的起点坐标lzm=lz(1:N)+dl/2;%每一小段的中点坐标lze=lz(2:N+1);%每一小段的终点坐标10%阻抗矩阵元素求解fi=2*log(dl/a)/dl-k*1i;fi_1=exp(-k*dl*1i)/dl;fi_2=exp(-k*2*dl*1i)/(2*dl);form=1:Nforn=1:Nifm==nfi11=fi;fi12=fi;fi13=fi;fi2=fi_1;fi3=fi_1;elseifabs(m-n)==1fi11=fi_1;fi12=fi_1;fi13=fi_1;ifnmfi2=fi_2;fi3=fi;elsefi3=fi_2;fi2=fi;endelsefi11=exp(-k*abs(m-n)*dl*1i)/(abs(m-n)*dl);fi12=exp(-k*abs(m-n)*dl*1i)/(abs(m-n)*dl);fi13=exp(-k*abs(m-n)*dl*1i)/(abs(m-n)*dl);ifnmfi2=exp(-k*(abs(n-m+1)*dl*1i))/((abs(n-m+1)*dl));fi3=exp(-k*(abs(n-m-1)*dl*1i))/((abs(n-m-1)*dl));elsefi2=exp(-k*(abs(m-n-1)*dl*1i))/((abs(m-n-1)*dl));fi3=exp(-k*(abs(m-n+1)*dl*1i))/((abs(m-n+1)*dl));endendz(m,n)=1i*w*u0/(4*pi)*dl*dl*fi11+(1/(1i*4*pi*w*e0))*(fi12-fi3-fi2+fi13);11endend%电压矩阵求解V=zeros(N,1);V((N+1)/2)=1;%计算电流系数矩阵I=z\V;%计算输入阻抗Z_in=1/I((N+1)/2);disp(['输入阻抗=',num2str(Z_in)]);%画方向图theta=0:0.01:2*pi;abs_f=zeros(1,length(theta));forn=1:1:Nabs_f=abs_f+I(n)*exp(k*(n*dl-L/2)*cos(theta)*1i);endabs_f=abs(sin(theta)*dl.*abs_f);Max_f=abs(sum(I)*dl);Far_patten2=abs_f/Max_f(1);polar(theta,Far_patten2,'-b');title('振子天线E面方向图');结果与分析：0.25λ输入阻抗==12.76051-440.6278i12E面方向图0.5λ输入阻抗=81.1696+36.6291iE面方向图1λ输入阻抗=1076.2926-1067.7693iE面方向图131.5λ输入阻抗=121.0432+101.4787iE面方向图14作业4请你谈谈对课程的建议、自己的收获等等。通过电磁场数值分析这门课程的学习，我了解掌握了许多运用matlab编程解决数值计算的方法，发现了这是一门与数学和计算机编程联系非常紧密的课程