# Python和R热释光动能朗伯W函数解析方程

<details>

<summary>关键词</summary>

Python | R | 数值计算 | 绘图 | 热释光 | 离子 | 数学 | 方程 | 常微分方程 | 求解器 | 活化能 | 频率因子 | 捕获率 | 模拟 | 峰值 | 加热速率 | 材料 | 固体 | 动能 | 卷积 | 郎伯W函数 | 模型 | 离域电子 | 异常 | 衰落 | 辉光曲线 | 反卷积 | 量子 | 局域跃迁 | 积分 | 求和

</details>

{% embed url="https://viadean.notion.site/Python-R-W-1dd02b3c90d147fabb64d321dd1ea8bd" %}

{% content-ref url="../zhi-dian-mi-jin-brief" %}
[zhi-dian-mi-jin-brief](../zhi-dian-mi-jin-brief)
{% endcontent-ref %}

### 🎯要点

1. Python和R计算和绘图：
   1. 🎯一陷阱一复合中心模型计算：🖊常微分方程求解器求解 | 🖊不同活化能和频率因子动力学参数下热释放光强度 | 🖊改变重捕获率 | 🖊数值解光热激发一阶动能微分方程 | 🖊模拟可变初始陷阱浓度的一阶峰和二阶峰 | 🖊数值计算寻找活化能 | 🖊绘制不同加热速率的辉光曲线 | 🖊加热速率法求活化能和频率因子 | 🖊数值求解和绘制一般动能解析方程 | 🖊朗伯W函数求解和绘制一般离子陷阱方程 | 🖊数值计算和绘制混合级动能微分方程。🎯辉光曲线反卷积数值计算：🖊光热激发一阶动能微分方程反卷积 | 🖊一般动能方程反卷积 | 🖊离域电子跃迁方程反卷积 | 🖊离域电子跃迁方程氧化铝热释光反卷积 |🖊混合级动能微分方程氧化铝热释光反卷积 | 🖊离域电子跃迁方程两个峰值辉光曲线数据反卷积 | 🖊一般动能方程两个峰值辉光曲线数据反卷积 | 🖊离域电子跃迁方程九峰值辉光曲线反卷积。
   2. 🎯量子隧道局域跃迁模型计算：🖊评估立方体中电子和受体之间距离的最近邻分布 | 🖊评估异常衰落周期结束时距离的分布 | 🖊数值计算评估基态隧道效应模型 | 🖊数值计算和绘制亨特利方程：基态隧道效应模型剩余电子 | 🖊数值评估同时基态隧道模型 | 🖊数值计算和绘制同时基态隧道模型自然界中同时照射和异常衰落 | 🖊模拟未衰落取样辉光曲线求和。🎯长石和磷灰石等材料的模型实验数据分析。
   3. 🎯等温释光信号模型：🖊模型数据分析 | 🎯离域电子跃迁光激发光模型：🖊模型数据分析 | 🎯其他模型

### 🍇Python朗伯W函数一阶时滞微分方程

延时系统是指在系统输入输入和结果输出之间存在显着时间延迟的系统，这种延迟可能是固有的或故意引入的。延时系统可以使用延迟微分方程进行建模。

假设以下系统包含时滞微分方程：

\


$$
\begin{cases}x^{\prime}(t)=a x(t)+a_d x(t-h)+b u(t) & t>0 \\ x(t)=g(t) & t \in[-h, 0) \\ x(t)=x_0 & t=0\end{cases}
$$

其中，

* $$x(t)$$ 是未知函数
* $$a_r a_d$$ 和 $$b$$ 是标量常数 $$\in R$$
* $$h$$ 是常数 $$\in R ^{+}$$，因此严格为正数并表示延迟
* g(t) 是一个函数，当时间变量 t 包含在区间 \[-h, 0) 中时，该函数提供 x(t) 的值
* $$x(0)=x_0$$ 是柯西初始条件

根据论文，时滞微分方程的解如下：

$$
x(t)=\sum_{k=-\infty}^{+\infty} e^{s_k t} C_k^I+\int_0^t \sum_{k=-\infty}^{+\infty} e^{s_k(t-\eta)} C_k^N b u(\eta) d \eta
$$

其中，



$$
C_k^I=\frac{x_0+a_d e^{-s_k h} \int_0^h e^{-s_k t} g(t-h) d t}{1+a_d h e^{-s_k h}}
$$

和

$$
C_k^N=\frac{1}{1+a_d h e^{-s_k h}}
$$

$$
s_k=\frac{1}{h} W_k\left(a_d h e^{-a h}\right)+a
$$

其中 $$W_k$$ 是索引 $$k$$ 的朗伯函数 $$W$$。朗伯函数 $$W$$ 是在随着索引 k 变化而获得的复数域中定义的函数族。

在开始Python实现之前需要注意两点：

* 函数W不能用初等函数表示，因此我们将使用SciPy提供的数值实现`scipy.special.lambertw`；
* 对于积分的计算，我们将始终使用 SciPy，特别是 `scipy.integrate.quad`，但请记住，我们在复数域中操作，我们必须小心分别对实部和虚部进行积分，因为 `scipy.integrate.quad` 不支持复数域中的积分。

```python
 import numpy as np
 from scipy import real, imag
 from scipy.integrate import quad
 from scipy.special import lambertw
 import matplotlib.pyplot as plt
```

本示例中使用的设置是：

```python
 t_begin=0.
 t_end=10.
 t_nsamples=101
 t_space, t_step = np.linspace(t_begin, t_end, t_nsamples, retstep=True)
 ​
 k_range=9
 a = 0.5
 ad = -2.5
 b = 1.75
 h = 1.
 g = lambda t : 1. - 0.1 * t
 u = lambda t : 0.2 * t
 x0 = 1.5
```

然后，

```python
 sk_fn = lambda k :  (1./h) * lambertw(ad * h * np.exp(-a * h), k) + a
 SK = [sk_fn(k) for k in range (-k_range, k_range+1)]
```

执行该片段后，变量 SK 包含所有各个 sk 的列表 $$𝑠_𝑘$$ 一劳永逸地预先计算。$$x(t)$$的实现在Python中如下：

```python
 def x(t):
     def integrand_for_cki(t_, sk):
         return np.exp(-sk * t_) * g(t_ - h)
 ​
     def integral_for_cki(sk):
         def real_func(t_, sk):
             return np.real(integrand_for_cki(t_, sk))
         def imag_func(t_, sk):
             return np.imag(integrand_for_cki(t_, sk))
         real_integral = quad(real_func, 0., h, args=(sk))
         imag_integral = quad(imag_func, 0., h, args=(sk))
         return real_integral[0] + 1.j*imag_integral[0]
 ​
     def integrand_for_x_t(eta):
         tot = 0.
         for k in range (-k_range, k_range+1):
             sk = SK[k + k_range]
             ck_denom = (1. + ad * h * np.exp(-sk * h))
             ckn = 1. / ck_denom
             tot += np.exp(sk * (t - eta)) * ckn * b * u(eta)
         return tot
 ​
     def integral_for_x_t():
         def real_func(eta):
             return np.real(integrand_for_x_t(eta))
         def imag_func(eta):
             return np.imag(integrand_for_x_t(eta))
         real_integral = quad(real_func, 0., t)
         imag_integral = quad(imag_func, 0., t)
         return real_integral[0] + 1.j*imag_integral[0]
 ​
     tot = 0.
     for k in range (-k_range, k_range+1):
         sk = SK[k + k_range]
         int_for_cki = integral_for_cki(sk)
         ck_denom = (1. + ad * h * np.exp(-sk * h))
         cki = (x0 + ad * np.exp(-sk * h) * int_for_cki) / ck_denom
         tot += np.exp(sk * t) * cki
     tot += integral_for_x_t()
     return tot
```

\


