1

Consider a differential equation, do be precise, the Polanyi-Wigner equation: 

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

This equation describes the desorption behavior of particles from a surface. An experimental technique called thermal desorption spectroscopy plots dy/dt thus making it possible to extract kinetic information such as the desorption Energy E, the desorption order n and the pre-exponential factor v from experimental data. It is oftentimes interesting to look upon the desorption energy E in dependence of dy/dt and y:

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

As we can see, a thermal desorption spectrum which is describable by eq (*) should, when eq (**) is applied, lead to a constant desorption energy E(y).

Its quite simple to solve eq (*) by applying 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]

However, if I try to solve eq (**) with this definition: 

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]

I run into problems, I only got a constant desorption energy when I assume a too large v (by a factor 10). Do I make some mistake?

--Edit: I have found a small mistake in the code of inv_pol_wig as I forgot to include the desorption order 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]

Unfortunately, this debug does not fix the problem.

-- Edit Nr. 2: Still no solution, though I made some troubleshooting and changed the 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]
Phil Ipp
  • 111
  • 4
  • 1
    Could we see the constant desortion value and when v is too large? – Katsu Sep 22 '13 at 20:23
  • @katsu: I have to say that I do not know when v is too large, but usually v has a value of about 10^13. Concerning the constant desorption value: do you mean what should come out or what comes out when applying inv_pol_wig? – Phil Ipp Sep 23 '13 at 18:43

0 Answers0