8476阅读
33回复

[RFEDA原创]Python的实现

UID ：63797

注册：2010-07-23
登录：2010-10-19
发帖：28
等级：仿真二级

0楼发表于: 2010-07-26 18:10:28

— 本帖被 tensor 执行加亮操作(2010-09-21) —

正在学习《Electromagnetic Simulation Using the FDTD Method》，Sullivan。

基本的一维FDTD包括两个迭代公式：

ex[k] = ex[k]+0.5*(hy[k-1]-hy[k])
hy[k] = hy[k]+0.5*(ex[k]-ex[k+1])

还需要一个激励源，如高斯脉冲：

pulse = exp(-0.5*((t0-T)/spread)**2)

fd1d_1.1.c改写为Python如下：

from pylab import *

KE = 200
ex = zeros(KE)
hy = zeros(KE)
kc = KE/2
t0 = 40.0
spread = 12
T = 0
NSTEPS = 100

for n in range(NSTEPS):
    T += 1
    for k in range(1,KE):
        ex[k] = ex[k]+0.5*(hy[k-1]-hy[k])
    pulse = exp(-0.5*((t0-T)/spread)**2)
    ex[kc] = pulse
    for k in range(0,KE-1):
        hy[k] = hy[k]+0.5*(ex[k]-ex[k+1])

plot(ex,label='Ex')
plot(hy,label='Hy') R( FQ+h

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

描述:fd1d_1.1

图片:fd1d_1.1.JPG

共2条评分

tensor	会员等级 +1	十分期待您分享下一贴！	2010-09-21
tensor	威望 +1	十分期待您分享下一贴！	2010-09-21

离线caocheng82

UID ：10116

注册：2008-03-28
登录：2025-05-26
发帖：697
等级：积极交流六级

1楼发表于: 2010-07-26 19:28:53

想请问一下python的代码相比与C的速度或者说效率如何？

解释型的语言的循环的效率是公认的差。

离线samuelp

UID ：63797

注册：2010-07-23
登录：2010-10-19
发帖：28
等级：仿真二级

2楼发表于: 2010-07-26 22:29:29

python肯定比c运行慢，也肯定比c开发快。python是胶水语言，很适合混合语言编程，可以充分利用各种语言的长处。

共1条评分

hefang

技术分 +1

积极参与讨论+技术分论坛感谢您的参与

2010-07-27

离线samuelp

UID ：63797

注册：2010-07-23
登录：2010-10-19
发帖：28
等级：仿真二级

3楼发表于: 2010-07-27 12:50:00

第二步是加入吸收边界条件:

ex_n[0] = ex_n-2[1]
ex_n[KE-1] = ex_n-2[KE-2]

预想的结果应该是电磁波到达边界后不会反射回来。

fd1d_1.2.c改写如下：

from pylab import *

KE = 200
ex = zeros(KE)
hy = zeros(KE)
kc = KE/2
t0 = 40.0
spread = 12
T = 0
NSTEPS = 250
ex_low_m1 = 0.0
ex_low_m2 = 0.0
ex_high_m1 = 0.0
ex_high_m2 = 0.0

for n in range(NSTEPS):
    T += 1

    for k in range(1,KE):
        ex[k] = ex[k]+0.5*(hy[k-1]-hy[k])

    pulse = exp(-0.5*((t0-T)/spread)**2)
    ex[kc] = pulse

    ex[0] = ex_low_m2
    ex_low_m2 = ex_low_m1
    ex_low_m1 = ex[1]
    ex[-1] = ex_high_m2
    ex_high_m2 = ex_high_m1
    ex_high_m1 = ex[-2]

    for k in range(0,KE-1):
        hy[k] = hy[k]+0.5*(ex[k]-ex[k+1])

plot(ex,label='T=250')
legend()
show()

图片:fd1d_1.2.JPG

离线samuelp

UID ：63797

注册：2010-07-23
登录：2010-10-19
发帖：28
等级：仿真二级

4楼发表于: 2010-07-27 14:22:19

如果存在对应于自由空间的相对介电常数，迭代公式中电场的系数由0.5修订为cb：

ex[k] = ex[k] + cb[k]*(hy[k-1]-hy[k])
hy[k] = hy[k] + 0.5 *(ex[k]-ex[k+1])

式中，

cb[k] = 0.5/epsilon

fd1d_1.3.c如下（右半边的相对介电常数为4）：

from pylab import *

KE = 200
ex = zeros(KE)
hy = zeros(KE)
kc = KE/2
t0 = 40.0
spread = 12
T = 0
NSTEPS = 320

ex_low_m1 = 0.0
ex_low_m2 = 0.0
ex_high_m1 = 0.0
ex_high_m2 = 0.0

epsilon = 4
kstart = 100
cb = 0.5*ones(KE)
cb[kstart:KE] = cb[kstart:KE]/epsilon

for n in range(NSTEPS):
    T += 1

    for k in range(1,KE):
        ex[k] = ex[k]+cb[k]*(hy[k-1]-hy[k])

    pulse = exp(-0.5*((t0-T)/spread)**2)
    ex[5] = ex[5]+pulse

    ex[0] = ex_low_m2
    ex_low_m2 = ex_low_m1
    ex_low_m1 = ex[1]
    ex[-1] = ex_high_m2
    ex_high_m2 = ex_high_m1
    ex_high_m1 = ex[-2]

    for k in range(0,KE-1):
        hy[k] = hy[k]+0.5*(ex[k]-ex[k+1])

plot(ex,label='T=320')
legend()
show()

在激励源中，有时为：

ex[kc] = pulse

有时为：

ex[kc] = ex[kc]+pulse

应该是大家讨论的硬源与软源的区别吧。

图片:fd1d_1.3.JPG

离线samuelp

UID ：63797

注册：2010-07-23
登录：2010-10-19
发帖：28
等级：仿真二级

5楼发表于: 2010-07-27 15:01:32

高斯脉冲源可以用其他激励源代替，如：

pulse = sin(2*pi*freq_in*dt*T)

fd1d_1.4.c改写如下：

from pylab import *

KE = 200
ex = zeros(KE)
hy = zeros(KE)
kc = KE/2
t0 = 40.0
spread = 12
T = 0
NSTEPS = 425

ex_low_m1 = 0.0
ex_low_m2 = 0.0
ex_high_m1 = 0.0
ex_high_m2 = 0.0

epsilon = 4
kstart = 100
cb = 0.5*ones(KE)
cb[kstart:KE] = cb[kstart:KE]/epsilon

ddx = 0.01
dt = ddx/(2*3e8)
freq_in = 700e6

for n in range(NSTEPS):
    T += 1

    for k in range(1,KE):
        ex[k] = ex[k]+cb[k]*(hy[k-1]-hy[k])

    pulse = sin(2*pi*freq_in*dt*T)
    ex[5] = ex[5]+pulse

    ex[0] = ex_low_m2
    ex_low_m2 = ex_low_m1
    ex_low_m1 = ex[1]
    ex[-1] = ex_high_m2
    ex_high_m2 = ex_high_m1
    ex_high_m1 = ex[-2]

    for k in range(0,KE-1):
        hy[k] = hy[k]+0.5*(ex[k]-ex[k+1])

plot(ex,label='T=425')
legend()
show()

图片:fd1d_1.4.JPG

离线samuelp

UID ：63797

注册：2010-07-23
登录：2010-10-19
发帖：28
等级：仿真二级

6楼发表于: 2010-07-27 15:51:21

如果介质是有损的，迭代公式修正为：

ex[k] = ca[k]*ex[k] + cb[k]*(hy[k-1]-hy[k])
hy[k] = hy[k] + 0.5*(ex[k]-ex[k+1])

式中：

eaf = dt*sigma/(2*epsz*epsilon)
ca[k] = (1-eaf)/(1+eaf)
cb[k] = 0.5/(epsilon*(1+eaf))

fd1d_1.5.c：

from pylab import *

KE = 200
ex = zeros(KE)
hy = zeros(KE)
kc = KE/2
t0 = 40.0
spread = 12
T = 0
NSTEPS = 500

ex_low_m1 = 0.0
ex_low_m2 = 0.0
ex_high_m1 = 0.0
ex_high_m2 = 0.0

epsilon = 4
kstart = 100

ddx = 0.01
dt = ddx/(2*3e8)
freq_in = 700e6

epsz = 8.85419e-12
sigma = 0.04
eaf = dt*sigma/(2*epsz*epsilon)
ca = ones(KE)
ca[kstart:KE] = (1-eaf)/(1+eaf)
cb = 0.5*ones(KE)
cb[kstart:KE] = 0.5/(epsilon*(1+eaf))

for n in range(NSTEPS):
    T += 1

    for k in range(1,KE):
        ex[k] = ca[k]*ex[k]+cb[k]*(hy[k-1]-hy[k])

    pulse = sin(2*pi*freq_in*dt*T)
    ex[5] = ex[5]+pulse

    ex[0] = ex_low_m2
    ex_low_m2 = ex_low_m1
    ex_low_m1 = ex[1]
    ex[-1] = ex_high_m2
    ex_high_m2 = ex_high_m1
    ex_high_m1 = ex[-2]

    for k in range(0,KE-1):
        hy[k] = hy[k]+0.5*(ex[k]-ex[k+1])

plot(ex,label='T=500')
legend()
show()

图片:fd1d_1.5.JPG

离线samuelp

UID ：63797

注册：2010-07-23
登录：2010-10-19
发帖：28
等级：仿真二级

7楼发表于: 2010-07-27 15:57:02

第一章的练习做完了，奖励自己一瓶冰红茶。

离线caocheng82

UID ：10116

注册：2008-03-28
登录：2025-05-26
发帖：697
等级：积极交流六级

8楼发表于: 2010-07-27 18:31:12

楼主能否介绍一下混合编程？比如C++和python，据我所知，这两种语言的配置异常复杂，对于没有LINUX使用经验的人来说，真得不好配置。而且还对python版本以及C++的编译器版本有限制。

离线samuelp

UID ：63797

注册：2010-07-23
登录：2010-10-19
发帖：28
等级：仿真二级

9楼发表于: 2010-07-27 23:23:03

安装pythonxy，所有配置一次自动搞定，windows,linux版本都行。混合编程可参考www。ibm。com/developerworks/cn/linux/l-pythc/

发帖回复