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

本科实验报告实验名称： 计算电磁学实验（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);grid onsubplot(212);w=1;x = -2*pi/w:0.01:2*pi/w;x1=w*x;plot(x,sin(x1);title(sin(wx);grid on2.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);grid on2.1.3k=10;w=1;u=2*pi/k;T=2*pi/w;x=-u:0.01:u;xlabel(x);ylabel(y(x);for t=0:0.1:2*T y=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);grid on 2.1.5n0=10;J0=2.35;dt=0.01;w=pi;i=1;for n=-10*n0:0.01:10*n0 %(U(ndt) if n0 u(i)=0; %i,i1 elseif n=n0; u(i)=n/n0; else u(i)=1; end i=i+1;endsubplot(211);n=-10*n0:0.01:10*n0;plot(n,u);xlabel(n);ylabel(U(nt);axis(-50 50 -1.5 1.5);grid on clear;%n0=10n0=10;J0=2.35;dt=0.01;w=pi;i=1;for n=-100*n0:dt:100*n0; if n0 Jz(i)=0; else if n*dt0 Jz(i)=0; elseif n*dt=n0 Jz(i)=J0*(n0-n)/n0*sin(w*n*dt); else Jz(i)=J0*sin(w*n*dt); end end i=i+1;endsubplot(212);n=-100*n0:0.01:100*n0;plot(n,Jz);xlabel(n);ylabel(Jz(n);grid on;2.1.6dt=0.01;x=-10:dt:10;y=sin(x);i=1;for x=-10:dt:10 dy=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;for x=-0:dt:pi-dt; s=sin(x+dt/2)*dt; y1=y1+s i=i+1;endy1 clear;dt=pi/120;x=-10:dt:10;y2=0;i=1;s=0;for x=-0:dt:2*pi-dt; s=sin(x+dt/2)*dt y2=y2+s; i=i+1;endy2得y1=2，y2=4.0804e-16clear;dt=pi/60;i=1;j=1for x=0:dt:4*pi-dt; yt=0; for t=0:dt:x-dt; s=sin(t+dt/2); yt=yt+s; i=i+1; end y(j)=yt; j=j+1;endx=0:dt:4*pi-dt;plot(x,y);xlabel(x);ylabel(y(x);grid on;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);for i=1:n u(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; %删双层循环的时候忘记移出来了for x=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;for y=6:dt:3-dt; y0=y+dt/2; s=(15-y0.2).*dt y1=y1+s; i=i+1;end%y1y2=0;s=0;for x=5:dt:3-dt; x0=x+dt/2; s=(3.*x0).*dt y2=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);grid off;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,-5E-7,-4E-7,-3E-7,-2E-7,-E-7,0,E-7,2E-7,3E-7,4E-7,5E-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;for w=-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;for n=-10:dt:400; i=i+1; if n0 F(i)=0; else F(i)=F0.*exp(-(n.*dt-t0).2/(2.*d0.2); endendn=-10:dt:400;plot(n,F);axis(-10 500 0 2.1);title(F(n); clear;F0=2;t0=4;d0=1;dt=0.1;i=0;for n=-10:1*dt:100 %老师讲n的步长要和dt相同 i=i+1; if n0 F(i)=0; else F(i)=F0.*exp(-(n.*dt-t0).2/(2.*d0.2); endendn=-10:1*dt:100;l=0;for w=-20:0.5*dt:20 l=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;for n=-10:1*dt:100 i=i+1; if n0 F(i)=0; else F(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;for w=-40:0.5*dt:40 l=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,xtickLabel,-5E-7,-4E-7,-3E-7,-2E-7,-E-7,0,E-7,2E-7,3E-7,4E-7,5E-7);xlabel(t);ylabel(y(t); subplot(3,2,3,4); %【以傅里叶法计算频谱】i=0;for w=-5e7:1e5:5e7; i=i+1; dF(i,:)=y.*exp(-1j.*w.*t).*dt;endF=sum(dF,2);w=-5e7:1e5:5e7;plot(w,abs(F);xlabel();ylabel(|F()|);% n=-N/Fs:1/Fs:N/Fs;% s=sin(w0.*n)+2.*cos(2.*w0.*n);subplot(3,2,5); %【以FFT函数算频谱】Y=fft(y,N);A=abs(Y);A=A/(N/2);%xinhaozhenfuA(1)=A(1)/2;%zhiliuchuxinhaozhenfuF=(1:N/2-1)*Fs / ( 2 * times ) / 10;% Fs/(2*times)就是f轴间隔,所以不需要除以N,除以10因为时域点数是10个周期的 /plot(F(1:N/2),A(1:N/2);xlim(0 1e7);xlabel(f / Hz);subplot(3,2,6); %【换成角频率的】FF = F * 2 * pi;plot(FF,A(1:N/2);xlim(0 5e7);xlabel();一维FDTDclear;cc = 3e8;m0 = 4*pi*1e-7;eps0 = 1/cc/cc/m0;freq = 3e6;%频率wavelength = cc/freq;%波长dx = wavelength/20;dt = dx/cc/2;maxz = 200;nmax = 300;%演示时间ex = zeros(maxz,1);hy = zeros(maxz,1);vc = 2:(maxz-1)for n = 1:nmax % 边界条件 ex1=ex(1); ex2=ex(2); exmaxz_1=ex(maxz-1); exmaxz=ex(maxz); ex(vc)=ex(vc)-dt/eps0/dx*(hy(vc)-hy(vc-1); ex(1)=ex2+(cc*dt-dx)/(cc*dt+dx)*(ex(2)-ex1); ex(maxz)=exmaxz_1+(cc*dt-dx)/(cc*dt+dx)*(ex(maxz-1)-exmaxz); ex(src)=cos(2*pi*freq*n*dt); hy1=hy(1); hy2=hy(2); hymaxz_1=hy(maxz-1); hymaxz=hy(maxz); hy(vc)=hy(vc)-dt/m0/dx*(ex(vc+1)-ex(vc); hy(1)=hy2+(cc*dt-dx)/(cc*dt+dx)*(hy(2)-hy1); hy(maxz)=hymaxz_1+(cc*dt-dx)/(cc*dt+dx)*(hy(maxz-1)-hymaxz); figure(1) plot(ex); title(ex); axis(1 maxz -2 2) pause(0.05);end二维FDTDZ=200; c=3*108; dx=0.005; dt=dx/(2*c); u0=(4*pi)*1e-7; s=(1/(36*pi)*1e-9; n0=100; Ca=1;Cb=dt/(s*dx); Da=1;Db=dt/(u0*dx); Ex=zeros(Z+1,Z+1);Ey=zeros(Z+1,Z+1);Hz=zeros(Z+1,Z+1); Hzzz=zeros(Z+1,Z+1); for n=0:400 for y=2:Z+1 %更新Ex for x=1:Z+1 Ex(x,y)=Ca*Ex(x,y)+Cb*(Hz(x,y)-Hz(x,y-1); end end for y=1:Z+1 %更新Ey for x=2:Z+1 Ey(x,y)=Ca*Ey(x,y)-Cb*(Hz(x,y)-Hz(x-1,y); end end for y=1:Z %更新Hz for x=1:Z Hz(n0,n0)=sin(1.20*pi*1011*n*dt); Hz(x,y)=Da*Hz(x,y)-Db*(Ey(x+1,y)-Ey(x,y)-Ex(x,y+1)+Ex(x,y); end end for y=2:Z %边界条件 Hz(1,y)=Hzzz(2,y)+(c*dt-dx)/(c*dt+dx)*(Hz(2,y)-Hzzz(1,y);