[求助]急寻有限差分法计算波导内电场分布的matlab程序

离线天不空

UID ：4999

注册：2007-09-18
登录：2011-07-20
发帖：78
等级：仿真一级

0楼发表于: 2008-05-29 10:47:24

哪位高手可提供一些有关利用有限差分法计算波导内电场分布的matlab程序的例子？？？
在书上只是看到计算的原理，具体到程序代码的部分就没有了，不知道该怎么编写
..

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

离线gwzhao

方恨少

UID ：17098

注册：2008-08-24
登录：2019-01-09
发帖：1374
等级：荣誉管理员

1楼发表于: 2008-09-01 12:24:43

function Fdtd_J()
clear;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%常数
mu0=4e-7*pi;
c0=3e8;
Z0=mu0*c0;
eps0=1/(Z0*c0);

%%%%%%%%%%%%%%%%%%%'
nb=10;
nx=30;
ny=30;
nz=30;
nt=120;
smax=4e8;
tpulse=30;

dx=1;
dy=1;
dz=1;
dt=1/(c0*sqrt(1/dx^2+1/dy^2+1/dz^2));

exs=zeros(nx+2*nb,ny+2*nb,nz+2*nb,2);
eys=zeros(nx+2*nb,ny+2*nb,nz+2*nb,2);
ezs=zeros(nx+2*nb,ny+2*nb,nz+2*nb,2);
hxs=zeros(nx+2*nb,ny+2*nb,nz+2*nb,2);
hys=zeros(nx+2*nb,ny+2*nb,nz+2*nb,2);
hzs=zeros(nx+2*nb,ny+2*nb,nz+2*nb,2);

sxe=zeros(1,nx+2*nb);
sye=zeros(1,ny+2*nb);
sze=zeros(1,nz+2*nb);
sxh=zeros(1,nx+2*nb);
syh=zeros(1,ny+2*nb);
szh=zeros(1,nz+2*nb);

for i=1:nx+2*nb
    sxe=smax*rho(i+1,nx,nb)^2;
    sxh=smax*rho(i+0.5,nx,nb)^2;
end

for i=1:ny+2*nb
    sye=smax*rho(i+1,ny,nb)^2;
    syh=smax*rho(i+0.5,ny,nb)^2;
end

for i=1:nz+2*nb
    sze=smax*rho(i+1,nz,nb)^2;
    szh=smax*rho(i+0.5,nz,nb)^2;
end

for n=1:nt
    if n<=tpulse
        F_in=-fun(n,tpulse);
        ezs(nx/2,ny/2,nz/2,:)=ezs(nx/2,ny/2,nz/2,:)-F_in/2;
    end

Ex=exs;
Ey=eys;
Ez=ezs;
Hx=hxs;
Hy=hys;
Hz=hzs;

[v]=add_diff_y(syh,-mu0,Hx(2:nx+2*nb,:,:,1),Ez(2:nx+2*nb,:,:,:));
Hx(2:nx+2*nb,:,:,1)=v;
[v]=add_diff_z(szh,mu0,Hx(2:nx+2*nb,:,:,2),Ey(2:nx+2*nb,:,:,:));
Hx(2:nx+2*nb,:,:,2)=v;

[v]=add_diff_z(szh,-mu0,Hy(:,2:ny+2*nb,:,1),Ex(:,2:ny+2*nb,:,:));
Hy(:,2:ny+2*nb,:,1)=v;
[v]=add_diff_x(sxh,mu0,Hy(:,2:ny+2*nb,:,2),Ez(:,2:ny+2*nb,:,:));
Hy(:,2:ny+2*nb,:,2)=v;

[v]=add_diff_x(sxh,-mu0,Hz(:,:,2:nz+2*nb,1),Ey(:,:,2:nz+2*nb,:));
Hz(:,:,2:nz+2*nb,1)=v;
[v]=add_diff_y(syh,mu0,Hz(:,:,2:nz+2*nb,2),Ex(:,:,2:nz+2*nb,:));
Hz(:,:,2:nz+2*nb,2)=v;

[v]=add_diff_y(sye,eps0,Ex(:,2:ny+2*nb,2:nz+2*nb,1),Hz(:,:,2:nz+2*nb,:));
Ex(:,2:ny+2*nb,2:nz+2*nb,1)=v;
[v]=add_diff_z(sze,-eps0,Ex(:,2:ny+2*nb,2:nz+2*nb,2),Hy(:,2:ny+2*nb,:,:));
Ex(:,2:ny+2*nb,2:nz+2*nb,2)=v;

[v]=add_diff_z(sze,eps0,Ey(2:nx+2*nb,:,2:nz+2*nb,1),Hx(2:nx+2*nb,:,:,:));
Ey(2:nx+2*nb,:,2:nz+2*nb,1)=v;
[v]=add_diff_x(sxe,-eps0,Ey(2:nx+2*nb,:,2:nz+2*nb,2),Hz(:,:,2:nz+2*nb,:));
Ey(2:nx+2*nb,:,2:nz+2*nb,2)=v;

[v]=add_diff_x(sxe,eps0,Ez(2:nx+2*nb,2:ny+2*nb,:,1),Hy(:,2:ny+2*nb,:,:));
Ez(2:nx+2*nb,2:ny+2*nb,:,1)=v;
[v]=add_diff_y(sye,-eps0,Ez(2:nx+2*nb,2:ny+2*nb,:,2),Hx(2:nx+2*nb,:,:,:));
Ez(2:nx+2*nb,2:ny+2*nb,:,2)=v;

end

%%%%%%%%%%%%%%%%%%%%
function  rho=rho(x,nmax,nb)
   rho=max(max(-x,x-nmax),0)/max(nb,1);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%
function g=g(t)
g=sin(2*t)*sin(t)^5;
end
%%%%%%%%%%%%%%%%%%%%%%
function fun=fun(t,tpulse)
fun=g(min(t/tpulse,1)*pi);
end
%%%%%%%%%%%%%%%%%%%%%%
function [v]=add_diff_x(s,pfac,v,w)
nsize2=size(v(1,:,1));
nsize2=nsize2(2);
alpha=zeros(1,nsize2);
beta=zeros(1,nsize2);
[alpha,beta]=coeff(pfac,dy,s,alpha,beta,nsize1);
nsize1=size(v(:,1,1));
nsize1=nsize1(2);
nsize3=size(v(1,1,:));
nsize3=nsize3(2);
for k=1:nsize3-1
    for j=1:nsize2-1
        for i=1:nzise1-1
            v(i,j,k)=alpha(i)*v(i,j,k)+beta(i)*(w(i+1,j,k,1)+w(i+1,j,k,2)-w(i,j,k,1)-w(i,j,k,2));
        end
    end
end

end
%%%%%%%%%%%%%%%%%%%%%%
function [v]=add_diff_y(s,pfac,v,w)
nsize2=size(v(1,:,1));
nsize2=nsize2(2);
alpha=zeros(1,nsize2);
beta=zeros(1,nsize2);
[alpha,beta]=coeff(pfac,dy,s,alpha,beta,nsize2);
nsize1=size(v(:,1,1));
nsize1=nsize1(2);
nsize3=size(v(1,1,:));
nsize3=nsize3(2);
for k=1:nsize3-1
    for j=1:nsize2-1
        for i=1:nzise1-1
            v(i,j,k)=alpha(j)*v(i,j,k)+beta(j)*(w(i,j+1,k,1)+w(i,j+1,k,2)-w(i,j,k,1)-w(i,j,k,2));
        end
    end
end

end
%%%%%%%%%%%%%%%%%%%%%%
function [v]=add_diff_z(s,pfac,v,w)
nsize1=size(v(:,1,1));
nsize1=nsize1(2);
nsize3=size(v(1,1,:));
nsize3=nsize3(2);
nsize2=size(v(1,:,1));
nsize2=nsize2(2);
alpha=zeros(1,nsize3);
beta=zeros(1,nsize3);
[alpha,beta]=coeff(pfac,dz,s,alpha,beta,nsize3);

for k=1:nsize3-1
    for j=1:nsize2-1
        for i=1:nzise1-1
            v(i,j,k)=alpha(k)*v(i,j,k)+beta(k)*(w(i,j,k+1,1)+w(i,j,k+1,2)-w(i,j,k,1)-w(i,j,k,2));
        end
    end
end

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [alpha,beta]=coeff(pfac,d,s,alpha,beta,n)
for i=1:n
    t=exp(-s(i)*dt);
    alpha(i)=t;
    if t>1-1e-7
        beta(i)=dt/(pfac*d);
    else
        beta(i)=(1-t)/(pfac*s(i)*d);
    end
end
end

逆流而上

离线gwzhao

方恨少

UID ：17098

注册：2008-08-24
登录：2019-01-09
发帖：1374
等级：荣誉管理员

2楼发表于: 2008-09-01 12:28:47

这个只是极其简单的例子，实际上数值模拟最难的是边界的处理。对三维而言，内存使用，速度等等。

其实我建议你不要向别人要源码，因为这样对你自己帮助很小很小。
而且简单的波导，数值计算就可以搞定了，没有必要用数值模拟。况且微波书上对波导讲了比较详细了，对简单情况可以直接得到的。

数值模拟主要针对边界比较复杂，无法解析求解的情况下，使用才有意义。

逆流而上

离线qiupan23

UID ：60223

注册：2010-05-23
登录：2010-06-28
发帖：17
等级：仿真新人

3楼发表于: 2010-06-28 16:13:14

这是时域有限差分法FDTD，不是FDM，FDTD是应用于时域的，FDM是应用于频域的

发帖回复

[求助]急寻 有限差分法计算波导内电场分布的matlab程序

[求助]急寻有限差分法计算波导内电场分布的matlab程序