1

考虑一个微分方程，准确地说，Polanyi-Wigner 方程：

dy/dt = - v*(y(t))^n * exp(-E/(kt))/b        (*)

该方程描述了颗粒从表面的解吸行为。一种称为热解吸光谱图的实验技术，因此可以从实验数据dy/dt中提取动力学信息，例如解吸能 E、解吸顺序n和指前因子。v根据 dy/dt 和 y 查看解吸能 E 通常很有趣：

E(y) = -kt*ln((dy/dt * b)/(v*(y(t))^n))    (**)

正如我们所看到的，由 eq ( *) 描述的热解吸光谱应该在应用 eq ( **) 时导致恒定的解吸能E(y)。

*通过应用来解决 eq ( ) 非常简单odeint：

def polwig_single(coverage,E,n,v, beta, Tmin,Tmax):
    k=8.6173/100000.
    def f(y, t):
        Si = y[0]
        # the Polanyi-Wigner Equation
        f0 = -v*np.power(Si,n)*np.exp(-E/(k*t))/beta
        return [f0]
    S0 = coverage        
    y0 = [S0]    
    t  = np.linspace(Tmin, Tmax, 10*(Tmax-Tmin))  
    soln = odeint(f, y0, t)
    S = soln[:, 0]
    S1=[0 if x<0 else x for x in S]
    pw=-np.gradient(S1)
    return[S1,pw]

但是，如果我尝试**用这个定义解决 eq ( )：

def inv_pol_wig(pw,S,n,v,beta,Tmin,Tmax):
    t  = np.linspace(Tmin, Tmax, 10*(Tmax-Tmin))   # time grid
    k=8.6173/100000# Boltzmann constant
    list1=[x / y for x,y in zip(pw,S)]
    list2=[beta * elem0 for elem0 in list1]
    list3=[elem1 / v for elem1 in list2]
    list4=[np.log(elem2) for elem2 in list3]
    Edes=[-(k * t) * elem3 for t,elem3 in zip(t,list4)]
    return[Edes]

我遇到了问题，当我假设一个太大v（10 倍）时，我只有一个恒定的解吸能量。我犯了一些错误吗？

--编辑：我在 inv_pol_wig 的代码中发现了一个小错误，因为我忘记包含解吸顺序n：

 def inv_pol_wig(pw,S,n,v,beta,Tmin,Tmax):
    t  = np.linspace(Tmin, Tmax, 10*(Tmax-Tmin))   # time grid
    k=8.6173/100000# Boltzmann constant
    list1=[x / np.power(y,n) for x,y in zip(pw,S)]
    list2=[beta * elem0 for elem0 in list1]
    list3=[elem1 / v for elem1 in list2]
    list4=[np.log(elem2) for elem2 in list3]
    Edes=[-(k * t) * elem3 for t,elem3 in zip(t,list4)]
    return[Edes]

不幸的是，此调试不能解决问题。

- 编辑编号。2：仍然没有解决方案，尽管我进行了一些故障排除并更改了 polwig_single-function：

def polwig_single(coverage,E,n,v, beta, Tmin,Tmax):
    k=8.6173/100000.
    def f(y, t):
        Si = y[0]
        # the Polanyi-Wigner Equation
        f0 = -v*np.power(Si,n)*np.exp(-E/(k*t))/beta
        return [f0]
    S0 = coverage        
    y0 = [S0]    
    t  = np.linspace(Tmin, Tmax, 10*(Tmax-Tmin))  
    soln = odeint(f, y0, t)
    S = soln[:, 0]
    S1=S
    pw=-np.gradient(S1)
    return[S1,pw]
4

0 回答 0