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

🎯要点

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

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

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

假设以下系统包含时滞微分方程：
$$\{\begin{array}{ll}{x}^{\prime }(t)=ax(t)+a_{d}x(t-h)+bu(t)& t>0\\ x(t)=g(t)& t\in [-h,0)\\ x(t)=x_0& t=0\end{array}$$
其中，

• $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}bu(\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)dt}{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}^{-ah}\right)+a$$

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

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

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

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

本示例中使用的设置是：

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

然后，

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)]

执行该片段后。$x(t)$的实现在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

参阅一：计算思维

参阅二：亚图跨际