自己按照教材上编了一个2维ADI FDTD程序,当把时间步Time设为200时不出结果,不知道哪里出了问题?刚学编程,望多指教!
%二维TE波ADI
clc;clear;
c=3e8;
eb0=8.854e-12;
miu0=4*3.1415926536e-7;
f=3e8;
dz=1/20;
dt=dz/2/c;
len=40;
ex=zeros(len,len+1);
ey=zeros(len+1,len);
hz=zeros(len,len);
%n+1/2时刻ex系数值
a=-dt*dt/(miu0*eb0*dz*dz);
b=1+2*a;
c=a;
%n+1/2时刻ex系数矩阵
A=zeros(len+1,len+1);
B=zeros(len,len+1);
A(1,1)=1;
A(len+1,len+1)=1;
for i=2:len
A(i,i-1)=a;
A(i,i)=b;
A(i,i+1)=c;
end
%n+1时刻ey系数矩阵
C=zeros(len+1,len);
Time=200;
for n=1:Time
%========================================过程一==================================
ex(len/2,len/2)=sin(2*pi*f*n*dt);
for i=1:len
for j=2:len
B(i,j)=ex(i,j)+dt/eb0/dz*(hz(i,j)-hz(i,j-1))-dt*dt/miu0/eb0/dz/dz*(ey(i+1,j)-ey(i,j)-ey(i+1,j-1)+ey(i,j-1));
end
end
ex=(inv(A)*B')';
for i=2:len
for j=1:len
ey(i,j)=ey(i,j)+dt/eb0/dz*(hz(i-1,j)-hz(i,j));
end
end
for i=1:len
for j=1:len
hz(i,j)=hz(i,j)-dt/miu0/dz*(ex(i,j)-ex(i,j+1)+ey(i+1,j)-ey(i,j));
end
end
%========================================过程二==================================
for i=2:len
for j=1:len
C(i,j)=ey(i,j)+dt/eb0/dz*(hz(i-1,j)-hz(i,j))-dt*dt/(miu0*eb0*dz*dz)*(ex(i,j+1)-ex(i,j)-ex(i-1,j+1)+ex(i-1,j));
end
end
for i=1:len
for j=2:len
ex(i,j)=ex(i,j)+dt/eb0/dz*(hz(i,j)-hz(i,j-1));
end
end
ey=inv(A)*C;
for i=1:len
for j=1:len
hz(i,j)=hz(i,j)-dt/miu0/dz*(ex(i,j)-ex(i,j+1)+ey(i+1,j)-ey(i,j));
end
end
end
mesh(ex);
%二维TE波ADI
clc;clear;
c=3e8;
eb0=8.854e-12;
miu0=4*3.1415926536e-7;
f=3e8;
dz=1/20;
dt=dz/2/c;
len=40;
ex=zeros(len,len+1);
ey=zeros(len+1,len);
hz=zeros(len,len);
%n+1/2时刻ex系数值
a=-dt*dt/(miu0*eb0*dz*dz);
b=1+2*a;
c=a;
%n+1/2时刻ex系数矩阵
A=zeros(len+1,len+1);
B=zeros(len,len+1);
A(1,1)=1;
A(len+1,len+1)=1;
for i=2:len
A(i,i-1)=a;
A(i,i)=b;
A(i,i+1)=c;
end
%n+1时刻ey系数矩阵
C=zeros(len+1,len);
Time=200;
for n=1:Time
%========================================过程一==================================
ex(len/2,len/2)=sin(2*pi*f*n*dt);
for i=1:len
for j=2:len
B(i,j)=ex(i,j)+dt/eb0/dz*(hz(i,j)-hz(i,j-1))-dt*dt/miu0/eb0/dz/dz*(ey(i+1,j)-ey(i,j)-ey(i+1,j-1)+ey(i,j-1));
end
end
ex=(inv(A)*B')';
for i=2:len
for j=1:len
ey(i,j)=ey(i,j)+dt/eb0/dz*(hz(i-1,j)-hz(i,j));
end
end
for i=1:len
for j=1:len
hz(i,j)=hz(i,j)-dt/miu0/dz*(ex(i,j)-ex(i,j+1)+ey(i+1,j)-ey(i,j));
end
end
%========================================过程二==================================
for i=2:len
for j=1:len
C(i,j)=ey(i,j)+dt/eb0/dz*(hz(i-1,j)-hz(i,j))-dt*dt/(miu0*eb0*dz*dz)*(ex(i,j+1)-ex(i,j)-ex(i-1,j+1)+ex(i-1,j));
end
end
for i=1:len
for j=2:len
ex(i,j)=ex(i,j)+dt/eb0/dz*(hz(i,j)-hz(i,j-1));
end
end
ey=inv(A)*C;
for i=1:len
for j=1:len
hz(i,j)=hz(i,j)-dt/miu0/dz*(ex(i,j)-ex(i,j+1)+ey(i+1,j)-ey(i,j));
end
end
end
mesh(ex);
