第三个程序：1流程图：2输出图形：在命令行输出时间步T为30,60,90,115，可分别得到：3Matlab程序：closeall;clear;clc;%设置网格数为60*60IE=60;JE=60;%确定源的位置ic=IE/2;jc=JE/2;%设置总场和散射场的边界ia=7;ib=IE-ia-1;ja=7;jb=JE-ja-1;ddx=0.01;%空间步长dt=ddx/6e8;%时间步长，其中6e8=2*cepsz=8.8e-12;%真空中的介电常数pi=3.14159;%进行初始化dz=zeros(IE,JE);ez=zeros(IE,JE);hx=zeros(IE,JE);hy=zeros(IE,JE);ihx=zeros(IE,JE);ihy=zeros(IE,JE);ga=ones(IE,JE);%自由空间中，其值为1ez_inc=zeros(1,JE);hx_inc=zeros(1,JE);ez_inc_low_m1=0;ez_inc_low_m2=0;ez_inc_high_m1=0;ez_inc_high_m2=0;%对PML中一位数组进行初始化gi2=ones(1,IE);gi3=ones(1,IE);fi1=zeros(1,IE);fi2=ones(1,IE);fi3=ones(1,IE);gj2=ones(1,JE);gj3=ones(1,JE);fj1=zeros(1,JE);fj2=ones(1,JE);fj3=ones(1,JE);%设置PML的层数npml=8;%计算PML的参数值，注意c语言与matlab中数组的区别，做适当修改fori=1:npml+1xnum=npml-(i-1);xd=npml;xxn=xnum/xd;xn=0.33*(xxn^3);gi2(i)=1.0/(1.0+xn);gi2(IE-i+1)=1.0/(1.0+xn);gi3(i)=(1.0-xn)/(1.0+xn);gi3(IE-i+1)=(1.0-xn)/(1.0+xn);xxn=(xnum-0.5)/xd;xn=0.25*(xxn^3);fi1(i)=xn;fi1(IE-i)=xn;fi2(i)=1.0/(1.0+xn);fi2(IE-i)=1.0/(1.0+xn);fi3(i)=(1.0-xn)/(1.0+xn);fi3(IE-i)=(1.0-xn)/(1.0+xn);endforj=1:npml+1xnum=npml-(j-1);xd=npml;xxn=xnum/xd;xn=0.33*(xxn^3);gj2(j)=1.0/(1.0+xn);gj2(IE-j+1)=1.0/(1.0+xn);gj3(j)=(1.0-xn)/(1.0+xn);gj3(IE-j+1)=(1.0-xn)/(1.0+xn);xxn=(xnum-0.5)/xd;xn=0.25*(xxn^3);fj1(j)=xn;fj1(JE-j)=xn;fj2(j)=1.0/(1.0+xn);fj2(JE-j)=1.0/(1.0+xn);fj3(j)=(1.0-xn)/(1.0+xn);fj3(JE-j)=(1.0-xn)/(1.0+xn);endt0=20;%脉冲的中心spread=8;%脉冲的宽度T=0;%时间初始值nsteps=1;%时间步初始值Tnd=str2num(input('Timesteps=','s'));%输入所需要的时间步数while(nsteps>0&&T<Tnd)fornsteps=1:TndT=T+1;%记录循环次数，得到时间步数forj=2:JEez_inc(j)=ez_inc(j)+0.5*(hx_inc(j-1)-hx_inc(j));end%入射缓冲吸收边界条件ez_inc(1)=ez_inc_low_m2;ez_inc_low_m2=ez_inc_low_m1;ez_inc_low_m1=ez_inc(2);ez_inc(JE)=ez_inc_high_m2;ez_inc_high_m2=ez_inc_high_m1;ez_inc_high_m1=ez_inc(JE-1);%计算电位移值forj=2:JEfori=2:IEdz(i,j)=gi3(i)*gj3(j)*dz(i,j)+gi2(i)*gj2(j)*0.5*(hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1));endend%加入高斯脉冲激励源pulse=e