660阅读
0回复

[求助]能帮忙看看我这easy的程序么？为什么出来的图总是发散的？

离线lovie

lovie

UID ：74467

注册：2011-03-25
登录：2011-05-30
发帖：11
等级：仿真新人

0楼发表于: 2011-04-30 10:09:52

200步是达到10^57了散的厉害

%定义几个常量
c=3*10^8;
e0=8.854*10^(-12);
u0=4*pi*10^(-7);

deta_x=0.001;%定义空间网格长度，为1/20波长
deta_y=0.001;
deta_z=0.001;
x_max=50;%定义总的网格数
y_max=50;
z_max=50;
%定义时间步长
deta_t=1/c/sqrt((1/deta_x)^2+(1/deta_y)^2+(1/deta_z)^2);
%定义迭代的次数
step=100;
%场量置零
hx1=zeros(x_max+1,y_max,z_max);
hx2=zeros(x_max+1,y_max,z_max);
hy1=zeros(x_max,y_max+1,z_max);
hy2=zeros(x_max,y_max+1,z_max);
hz1=zeros(x_max,y_max,z_max+1);
hz2=zeros(x_max,y_max,z_max+1);
ez1=zeros(x_max+1,y_max+1,z_max);
ez2=zeros(x_max+1,y_max+1,z_max);
ex1=zeros(x_max,y_max+1,z_max+1);
ex2=zeros(x_max,y_max+1,z_max+1);
ey1=zeros(x_max+1,y_max,z_max+1);
ey2=zeros(x_max+1,y_max,z_max+1);
%定义电流源的位置
source_x=25;
source_y=25;
source_z=25;
%定义观察点的位置
field_x=23;
field_y=23;
field_z=23;
%定义激励源
%时谐源
frequency=15*10^9;%波长为0.02
omiga=2*pi*frequency;
%高斯源
t0=0.0667*10^(-9);
tt=0.0667*10^(-9);
%定义保存结果的变量
result=zeros(step,1,1);
%迭代开始
for t=1:1:step
    t
   % current(t)=sin(omiga*t*deta_t);
current(t)=exp(-4*pi*(t*deta_t-t0)^2/tt^2);
%ex的计算；
for i=1:1:x_max
    for j=2:1:y_max
        for k=2:1:z_max
            if i==source_x&j==source_y&k==source_z%判断是不是源点
               ex2(i,j,k)=ex1(i,j,k)+deta_t/e0*((hz2(i,j,k)-hz2(i,j-1,k))/deta_y-(hy2(i,j,k)-hy2(i,j,k-1))/deta_z)-deta_t/e0*current(t);
            else
               ex2(i,j,k)=ex1(i,j,k)+deta_t/e0*((hz2(i,j,k)-hz2(i,j-1,k))/deta_y-(hy2(i,j,k)-hy2(i,j,k-1))/deta_z);
            end
        end
    end
end

%ey的计算；
for i=2:1:x_max
    for j=1:1:y_max
        for k=2:1:z_max
            if i==source_x&j==source_y&k==source_z%判断是不是源点
               ey2(i,j,k)=ey1(i,j,k)+deta_t/e0*((hx2(i,j,k)-hx2(i,j,k-1))/deta_z-(hz2(i,j,k)-hz2(i-1,j,k))/deta_x)-deta_t/e0*current(t);
            else
               ey2(i,j,k)=ey1(i,j,k)+deta_t/e0*((hx2(i,j,k)-hx2(i,j,k-1))/deta_z-(hz2(i,j,k)-hz2(i-1,j,k))/deta_x);
            end
        end
    end
end
%ez的计算；
for i=2:1:x_max
    for j=2:1:y_max
        for k=1:1:z_max
            if i==source_x&j==source_y&k==source_z%判断是不是源点
                ez2(i,j,k)=ez1(i,j,k)+deta_t/e0*((hy2(i,j,k)-hy2(i-1,j,k))/deta_x-(hx2(i,j,k)-hx2(i,j-1,k))/deta_y)-deta_t/e0*current(t);
            else
                ez2(i,j,k)=ez1(i,j,k)+deta_t/e0*((hy2(i,j,k)-hy2(i-1,j,k))/deta_x-(hx2(i,j,k)-hx2(i,j-1,k))/deta_y);
            end
        end
    end
end
%hx的计算；
for i=1:1:x_max+1
    for j=1:1:y_max
        for k=1:1:z_max
            hx2(i,j,k)=hx1(i,j,k)-deta_t/u0*((ez1(i,j+1,k)-ez1(i,j,k))/deta_y-(ey1(i,j,k+1)-ey1(i,j,k))/deta_z);
        end
    end
end
%hy的计算；
for i=1:1:x_max
    for j=1:1:y_max+1
        for k=1:1:z_max
            hy2(i,j,k)=hy1(i,j,k)-deta_t/u0*((ex1(i,j,k+1)-ex1(i,j,k))/deta_z-(ez1(i+1,j,k)-ez1(i,j,k))/deta_x);
        end
    end
end
%hz的计算；
for i=1:1:x_max
    for j=1:1:y_max
        for k=1:1:z_max+1
            hz2(i,j,k)=hz1(i,j,k)-deta_t/u0*((ey1(i+1,j,k)-ey1(i,j,k))/deta_x-(ex1(i,j+1,k)-ey1(i,j,k))/deta_y) ..

未注册仅能浏览部分内容，查看全部内容及附件请先登录或注册

发帖回复