1916阅读
4回复

[已解决]Sullivan书中的傅里叶变化程序总运行不对，版主帮忙看看！

UID ：33555

注册：2009-05-25
登录：2010-04-25
发帖：56
等级：仿真一级

0楼发表于: 2010-03-02 18:17:21

Sullivan书中2.2.c程序中傅里叶变化中入射波的变换计算部分，代码如下：
    /* Fourier Transform of the input pulse */
    if (T<100)
    {
     for (m=0; m<=2; m++)
     {
       real_in[m] = real_in[m] + cos (arg[m]*T)*ex[10];
       imag_in[m] = imag_in[m] - sin (arg[m]*T)*ex[10];
     }
    }入射波的位置是100，这里ex[10]是不是应该改成100？个人对书中程序觉得不正确的地方。

另：我是用MATLAB的，所以将傅里叶变换部分改成MATLAB程序运算，得不到正确的结果，各位帮忙看看。
我的代码是在Sullivan书中一维有耗介质传播程序的中加入傅里叶变化部分，如下：
clc,clear
% 空间时间步长
delta_x =0.01;
delta_t = delta_x/(2*3e8);
% 空间时间网格数
max_time = 400;
max_space = 100;
% 相关参数设置
sigma = 0.04;
epszero = 8.85419e-12;
muz=4.0*pi*1.0e-7;
frequency=2.45e9;
omega = 2*pi*frequency;
dielectric_constant = 4;
% 初始化电磁场变量
E = zeros(max_space,1);
H = E;
%传播波的傅里叶变化变量
real_pt = zeros(max_space,1); %实部变量
imag_pt = zeros(max_space,1); %虚部变量
ampn = zeros(max_space,1);   %振幅标量
phasen = zeros(max_space,1); %相位变量
%入射波的傅里叶变化变量
real_in = 0;
imag_in = 0;
amp_in =0;
phase_in =0;
% 介质位置设置
dielectric_begin =50;
dielectric_end = max_space;
dielectic_space_a = ones(max_space,1);
dielectic_space_b = ones(max_space,1);
% 方程参数设置
eta = 0.5;
for i = 1:max_space
    if (i>=dielectric_begin & i <= dielectric_end)
      loss_term = (delta_t*sigma) / (2*epszero*dielectric_constant);
      dielectric_space_a(i)=(1-loss_term)/(1+loss_term);
      dielectric_space_b(i)=0.5/(dielectric_constant *(1+loss_term));
    else
        dielectric_space_a(i)=1;
        dielectric_space_b(i)=eta;
    end
end
%ABS needed for LHS of problem space do not work with Lossy medium
E_low_m2 = 0;
E_low_m1 = 0;
E_high_m2 = 0;
E_high_m1 = 0;
%Initialize Pulse Variables（初始化脉冲变量）
place_pulse = 2;
% spread = 12;
% center_problem_space = max_space/2;
% t0 = 40.0;
%Set up dielectric medium within the problem space
for n = 1:max_time
    for k = 2:max_space
        E(k) = dielectric_space_a(k)*E(k) + dielectric_space_b(k)*(H(k-1)-H(k));
    end
% 傅里叶计算
    for k = 1:max_space
        real_pt(k) = real_pt(k) + cos(2*pi*frequency*delta_t*n)*E(k); %对空间每一点计算实部值
        imag_pt(k) = imag_pt(k) - sin(2*pi*frequency*delta_t*n)*E(k); %对空间每一点计算虚部值

    end
% 入射波傅里叶计算
    if n<100
        real_in = real_in+cos(2*pi*frequency*delta_t*n)*E(place_pulse);
        imag_in = imag_in-sin(2*pi*frequency*delta_t*n)*E(place_pulse);
    end

    pulse = sin(2*pi*frequency*delta_t*n);
    E(place_pulse) = E(place_pulse) + pulse;

    E(1) = E_low_m2;
    E_low_m2 = E_low_m1;
    E_low_m1 = E(2);
    E(max_space) = E_high_m2;
    E_high_m2 = E_high_m1;
    E_high_m1 = E(max_space-1);
    for h = 1:max_space-1
        H(h) = H(h) +eta*(E(h)-E(h+1));
    end
    plot(E)
  &n ..

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

离线cem-uestc

UID ：9061

注册：2008-03-07
登录：2019-01-05
发帖：2575
等级：荣誉管理员

1楼发表于: 2010-03-02 20:13:01

for k = 1:max_space ??
for (m=0; m<=2; m++) ??

看看比两个的循环

欢迎光临
http://www.mwtee.com/home.php?mod=space&uid=13535

离线lizi0908

UID ：33555

注册：2009-05-25
登录：2010-04-25
发帖：56
等级：仿真一级

2楼发表于: 2010-03-02 21:19:09

版主可能没有看过Sullivan书中的C程序，Sullivan书中的C程序激励源是脉冲信号，然后对100MHZ，200MHZ，500MHZ，三个频率进行分析，m变量代表的是这个频率。我的程序只有一个频率，所以没有用m变量。Sullivan书中的C程序放在附件中，辛苦版主看看。我对书中的理论也看了很多次了，还是没找出问题。

附件:

FDTD.doc (43 K) 下载次数:10

离线lizi0908

UID ：33555

注册：2009-05-25
登录：2010-04-25
发帖：56
等级：仿真一级

3楼发表于: 2010-03-02 23:29:18

刚才把论坛的帖子翻了一遍，看到也有网友出现过这样的问题，帖子地址：
http://bbs.rfeda.cn/read-htm-tid-27893.html
http://bbs.rfeda.cn/read-htm-tid-19840-page-1.html
后面是怎么解决的没说啊。
大家帮帮忙吧！

离线lizi0908

UID ：33555

注册：2009-05-25
登录：2010-04-25
发帖：56
等级：仿真一级

4楼发表于: 2010-03-03 18:16:35

自己认为ampn变量即幅度值最后应该除以总时间步长数即max_time变量。不除的话时间越长，幅度值越大，这样应该不合理吧！

发帖回复