北理工-计算电磁学实验报告(2017.6)

本科实验报告实验名称：计算电磁学实验（MATLAB实现）课程名称：计算电磁学实验时间：周三上午8:00-9:35任课教师：宋巍实验地点：信教2004实验教师：宋巍实验类型：√原理验证□综合设计□自主创新学生姓名：代明东学号/班级：组号：学院：同组搭档：专业：成绩：2.1.1subplot(211);w=2;x=-2*pi/w:0.01:2*pi/w;x1=w*x;plot(x,cos(x1));title('cos(wx)');gridonsubplot(212);w=1;x=-2*pi/w:0.01:2*pi/w;x1=w*x;plot(x,sin(x1));title('sin(wx)');gridon2.1.2k=1;w=1;u=2*pi/k;T=2*pi/w;x=-u:0.01:u;plot(x,cos(x));title('cos(-kx)');gridon2.1.3k=10;w=1;u=2*pi/k;T=2*pi/w;x=-u:0.01:u;xlabel('x');ylabel('y(x)');fort=0:0.1:2*Ty=cos(w*t-k*x);plot(x,y);pause(0.1);end;2.1.4w=2;x=-20*pi/w:pi/2:20*pi/w;x1=w*x;plot(x,cos(x1),'.');ylim([-1.5,1.5]);title('cos(wx),¼ä¸ôÈ¡¦Ð/2');gridon2.1.5n0=10;J0=2.35;dt=0.01;w=pi;i=1;forn=-10*n0:0.01:10*n0%(U(ndt)ifn0u(i)=0;%iÊÇÊý×éµÄË÷Òý£¬Ö»ÄÜÊÇÕýÕûÊý,iÒªÉè³õʼֵ£¬³õʼֵֻÄÜΪ1elseifn=n0;u(i)=n/n0;elseu(i)=1;endi=i+1;endsubplot(211);n=-10*n0:0.01:10*n0;plot(n,u);xlabel('n');ylabel('U(n¡÷t)');axis([-5050-1.51.5]);gridonclear;%È¡n0=10µÄʱºòÇ°ÃæµÄÊý¾ÝÒªÇåÁ㣬²»È»Êý×é¹æÄ£»á´ÓÉÏÃæ¼Ì³ÐÏÂÀ´n0=10;J0=2.35;dt=0.01;w=pi;i=1;forn=-100*n0:dt:100*n0;ifn0Jz(i)=0;elseifn*dt0Jz(i)=0;elseifn*dt=n0Jz(i)=J0*(n0-n)/n0*sin(w*n*dt);elseJz(i)=J0*sin(w*n*dt);endendi=i+1;endsubplot(212);n=-100*n0:0.01:100*n0;plot(n,Jz);xlabel('n');ylabel('Jz(n)');gridon;2.1.6dt=0.01;x=-10:dt:10;y=sin(x);i=1;forx=-10:dt:10dy=sin(x+dt/2)-sin(x-dt/2);d(i)=dy/dt;i=i+1;endx=-10:dt:10;plot(x,d);title('sin(x)的导数');xlabel('x');ylabel('dy/dx');2.1.7dt=0.01;x=-10:dt:10;i=1;y1=0;forx=-0:dt:pi-dt;s=sin(x+dt/2)*dt;y1=y1+si=i+1;endy1clear;dt=pi/120;x=-10:dt:10;y2=0;i=1;s=0;forx=-0:dt:2*pi-dt;s=sin(x+dt/2)*dty2=y2+s;i=i+1;endy2得y1=2，y2=4.0804e-16clear;dt=pi/60;i=1;j=1forx=0:dt:4*pi-dt;yt=0;fort=0:dt:x-dt;s=sin(t+dt/2);yt=yt+s;i=i+1;endy(j)=yt;j=j+1;endx=0:dt:4*pi-dt;plot(x,y);xlabel('x');ylabel('y(x)');gridon;2.1.8【用meshgrid】[x,y]=meshgrid(-1:0.4:5,-1:0.4:3);u=x.*y;v=3.*x-y.^2;figurequiver(x,y,u,v)【不用meshgrid】x=-1:0.4:5;y=-1:0.4:3;m=length(x);%列n=length(y);%行u=zeros(n,m);fori=1:nu(i,:)=x*y(i);%':'表所有，u（列，行）v(i,:)=3*x-y(i)^2;endquiver(u,v)2.1.9（1）dt=-0.01;i=1;j=1;y1=0;%删双层循环的时候忘记移出来了forx=5:dt:3-dt;x0=x+dt/2;s=(1.5.*x0.*(x0-1)+(3.*x0-2.25.*(x0-1)^2).*1.5).*dt%大中小括号在Matl里含义不同，数学表达式一律用小括号y1=y1+s;i=i+1;endy1y1=-10.0000（2）dt=-0.01;i=1;y1=0;s=0;fory=6:dt:3-dt;y0=y+dt/2;s=(15-y0.^2).*dty1=y1+s;i=i+1;end%y1y2=0;s=0;forx=5:dt:3-dt;x0=x+dt/2;s=(3.*x0).*dty2=y2+s;i=i+1;end%y2y=y1+y2;yy=-6.00002.1.10dt=0.1;[x,y]=meshgrid(-1-dt:dt:5+dt,-1-dt:dt:3+dt);%meshgrid就是构建二维数组的函数dux=((x+dt).*y-(x-dt).*y)/(2*dt);dvy=((3.*x-(y+dt).^2)-(3.*x-(y-dt).^2))/(2.*dt);divF=dux+dvy;contour3(x,y,divF,250);title('É¢¶È');xlabel('x');ylabel('y');zlabel('z');figure(2);%单开一个网格dt=3.*dt;[x,y,z]=meshgrid(-1-dt:dt:5+dt,-1-dt:dt:3+dt,-10:dt:10);Size=size(x);%计算x矩阵每维的长度，将长度记在一个矩阵里empty=zeros(Size(1),Size(2),Size(3));%建立一个和Size矩阵一样大的空矩阵p=empty;q=empty;r=((3.*(x+dt)-y.^2)-(3.*(x-dt)-y.^2))/(2.*dt)-(x.*(y+dt)-x.*(y-dt))/(2.*dt);quiver3(x,y,z,p,q,r,'color',[0,0,0]./255);%quiver3是建立向量三维坐标，contour3是普通三维网格或者标量图title('旋度');xlabel('x');ylabel('y');zlabel('z');gridoff;2.1.11clear;f0=3e6;w0=2*pi*f0;T=1/f0;dt=1e-9*pi;t=-5*T-dt:dt:5*T+dt;y=sin(w0*t)+2*cos(2*w0*t);subplot(211);plot(t,y);set(gca,'XTick',[-5e-7*pi:1e-7*pi:5e-7*pi]);set(gca,'xtickLabel',{'-5¦ÐE-7','-4¦ÐE-7','-3¦ÐE-7','-2¦ÐE-7','-¦ÐE-7','0','¦ÐE-7','2¦ÐE-7','3¦ÐE-7','4¦ÐE-7','5¦ÐE-7'});%GCA用法示例set(gca,'box','on','xlim',[02*pi],'YDir','reverse')后，图形变成下图所示，出现坐标边界（box），x轴显示坐标范围缩小（xlim）,y轴方向反转（ydir）xlabel('t');ylabel('y(t)');subplot(212);i=0;forw=-5e7:1e5:5e7;%取值和之前的y计算的取值点必须相同，不然没法用上一问算出的y矩阵i=i+1;dF(i,:)=y.*exp(-1j.*w.*t).*dt;%t和y都是上一问已经得到的矩阵，此处用上可以省略很多运算endF=sum(dF,2);%sum(x)列求和,sum(x,2)行求和,sum(x(:))矩阵求和w=-5e7:1e5:5e7;plot(w,abs(F));xlabel('w');ylabel('|F(w)|');2.1.12clear;F0=2;t0=0.1;d0=1;dt=0.01;i=0;forn=-10:dt:400;i=i+1;ifn0F(i)=0;elseF(i)=F0.*exp(-(n.*dt-t0).^2/(2.*d0.^2));endendn=-10:dt:400;plot(n,F);axis([-1050002.1]);title('F(n)');clear;F0=2;t0=4;d0=1;dt=0.1;i=0;forn=-10:1*dt:100%老师讲n的步长要和dt相同i=i+1;ifn0F(i)=0;elseF(i)=F0.*exp(-(n.*dt-t0).^2/(2.*d0.^2));endendn=-10:1*dt:100;l=0;forw=-20:0.5*dt:20l=l+1;dF(l,:)=F.*exp(-1j.*w.*n.*dt).*dt;endFw=sum(dF,2);w=-20:0.5*dt:20;plot(w,abs(Fw));xlabel('w');ylabel('F(w)');2.1.13clear;subplot(211);F0=2;t0=4;d0=1;dt=0.1;i=0;w0=20;forn=-10:1*dt:100i=i+1;ifn0F(i)=0;elseF(i)=F0.*exp(-(n.*dt-t0).^2/(2.*d0.^2)).*cos(w0.*n.*dt);endendplot(F);title('时域');subplot(212);n=-10:1*dt:100;l=0;forw=-40:0.5*dt:40l=l+1;dF(l,:)=F.*exp(-1j.*w.*n.*dt).*dt;endFw=sum(dF,2);w=-40:0.5*dt:40;plot(w,abs(Fw));xlabel('w');ylabel('F(w)');title('频域');2.1.11FFTw0=128;Fs=256;%采样频率N=256;%采样点数n=[0:1/Fs:N/Fs];%采样时间间隔作为坐标s=sin(w0.*n)+2.*cos(2.*w0.*n);Y=fft(s,N);Ayy=abs(Y);%取模Ayy=Ayy/(N/2);%换算实际幅度Ayy(1)=Ayy(1)/2;F=([1:N]-1)*Fs/N;%换算实际频率值stem(F(1:N/2),Ayy(1:N/2));%显示换算后模值结果title('FFT算频谱');clear;f0=3e6;w0=2*pi*f0;T=1/f0;dt=1e-9*pi;t=-5*T-dt:dt:5*T+dt;Fs=f0;%这个就是保证FFT出来每个点的间隔代表1Fs/(2*times)//%所以Fs要关联时域点数即关联T因为length(t)是通过T限定的//times=4;N=times*length(t)*2;%最小取样点数1064*2//y=sin(w0*t)+2*cos(2*w0*t);%【时域图】subplot(3,2,[1,2]);plot(t,y);set(gca,'XTick',[-5e-7*pi:1e-7*pi:5e-7*pi]);set(gca,'