【你购物我买单】Arduino UNO Q 数值计算

感谢 DigiKey 和电源网共同举办的【你购物我买单】活动，我采购的物料是 Arduino UNO Q 开发板，料号 ABX00162 .

本文介绍了 Arduino UNO Q 开发板结合 scipy 实现微分方程组的数值求解的项目设计。

项目介绍

Arduino UNO Q 开发板结合 scipy 与 matplotlib 软件库实现微分方程组的数值求解。

准备工作：硬件连接、系统安装、软件更新等；

环境搭建：科学计算所需软件库的安装、测试等；

数值计算：结合实际物理问题进行数值求解和计算，包括流程图、代码、效果演示等；

原子系统：考虑全量子 JC 模型，数值计算 Rabi 振荡动力学及其参量依赖关系。

准备工作

包括硬件连接、系统安装、软件更新等。

硬件连接

这里采用 SSH 远程控制，使用 Type-C 数据线供电并 WiFi 联网即可。

软件更新

更新软件包

sudo apt update
sudo apt upgrade

环境搭建

安装 scipy 和 matplotlib 软件包，便于调用 RK45 算法和科学绘图；

sudo apt install python3-scipy
sudo apt install python3-matplotlib

Runge-Kutta 算法

为了获得更为精确的数值解，常用方案是采用四阶 Runge-Kutta 方法。

使用 Python 编程，通常采用 odeint 和 solve_ivp 函数实现；

solve_ivp 是 Python 中 SciPy 库提供的一个函数，用于求解初值问题 (IVP, Initial Value Problem) 的常微分方程组 (ODEs)。

使用方法：

scipy.integrate.solve_ivp(fun, t_span, y0, method='RK45', t_eval=None, vectorized=False, args=None, **options)

数值计算

通过调用 solve_ivp 函数实现常微分方程组的数值求解。

Lorenz 吸引子

Lorenz attractor 是混沌理论中的经典模型。

该模型源于对大气对流方程的简化，旨在研究流体力学中的混沌现象，常用于描述 “蝴蝶效应” 。

该模型的简化数学方程为

可使用 4 阶 Runge-Kutta 算法进行数值求解。

代码

终端执行 touch ode_demo.py 指令新建文件，并添加如下代码

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# 洛伦兹方程组
def lorenz(t, X, sigma, rho, beta):
    x, y, z = X
    return np.array([
        sigma * (y - x),          # dx/dt
        x * (rho - z) - y,        # dy/dt
        x * y - beta * z          # dz/dt
    ])

# 参数与初值
sigma = 1
beta  = 3
rho   = 0
x0    = [0, 0, 0]
t_span = (0, 10)                    # 积分区间
t_eval = np.arange(0, 10, 0.01)  # 固定步长 0.01 输出

# 数值计算
sol = solve_ivp(lorenz, t_span, x0, args=(sigma, rho, beta),
                t_eval=t_eval, rtol=1e-4, atol=1e-4)

# 绘图
fig = plt.figure(figsize=(6, 5))
ax = fig.add_subplot(111, projection='3d')
ax.plot(sol.y[0], sol.y[1], sol.y[2], lw=0.5)
ax.set_xlabel('x'); ax.set_ylabel('y'); ax.set_zlabel('z')
ax.set_title('Lorenz Attractor')
plt.show()

终端运行 python lorentz_attractor.py 弹窗显示结果

二能级系统

考虑二能级系统光学 Bloch 方程组的数值求解。

Hamiltonian

系统示意图如下

系统 Hamiltonian 为

密度矩阵方程组

结合 Lindblad 主方程，给出矩阵元对角元的运动方程

详见：电磁诱导透明机制下基于四波混频过程的光学参量放大动力学研究 .

代码

终端执行指令 touch twolevel_system.py 新建程序文件，添加如下代码

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# ---------- 参数 ----------
Delta = 1.0
Omega = 1
gamma = 1
R0 = np.array([0.0, 0.0])      # 初值 [Re(ρ), Im(ρ)]
t_span = (0, 10)                # 时间取值范围

# ---------- ODE define ----------
def rk(t, R, Delta, Omega, gamma):
    r1, r2 = R
    dr1 = -(gamma*r2 + Delta*r1)
    dr2 = -(gamma*r1 - Delta*r2)
    return np.array([dr1, dr2])

# ---------- Solve ----------
sol = solve_ivp(rk, t_span, R0, args=(Delta, Omega, gamma),
                method='RK45', rtol=1e-6, atol=1e-6, dense_output=True)

t = sol.t          # 时间轴
Rho = sol.y.T      # 数组 (N,2)

# ---------- 画图 ----------
plt.plot(t, R[:, 0], '-k', label=r"$\chi'$", linewidth=3)
plt.plot(t, R[:, 1], '--b', label=r"$\chi''$", linewidth=3)
plt.xlabel(r'Time ($1/\gamma$)')
plt.ylabel(r'$\chi$')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

效果

终端运行 python twolevel_system.py ，输出计算进度

弹窗显示绘图结果

获得介质极化率的实部（色散）和虚部（吸收）对应的光谱。

对应矩阵元的时间演化曲线如下

矩阵元随时间的演化快速趋于稳定。

全量子模型

进一步考虑全量子的 Jaynes-Cummings 模型。

系统示意图如下

主方程

结合二能级系统的全量子 Hamiltonian 和密度矩阵的 Lindblad 主方程，可得

同样使用 Runge-Kutta 算法，结合 solve_ivp 函数，给出全数值解决方案。

代码

终端执行指令 touch jc_p-t.py 新建程序文件，添加如下代码

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# ---------------- 参数 ----------------
omega   = 1
Lambda1 = 1
Lambda2 = 1
delta   = 0
F       = 1
A       = 1
Gm      = np.array([0, 2, 5, 1])   # Γ 取值
t_span  = (0, 100)                       # 积分区间
t_eval  = np.arange(0, 10, 0.1)       # 时间
# ---------------------------------------

# 微分方程组
def rhs(t, R, delta, omega, Lambda1, Lambda2, Gamma, gamma1, gamma2, A):
    dR = np.zeros(4)
    dR[0] = Lambda1 - gamma1*R[1] + A*R[2]    # drho11/dt
    dR[1] = Lambda2 + 2*omega*R[0]            # drho22/dt
    dR[2] = -Gamma*R[1] - delta*R[3]                        # dRe(rho21)/dt
    dR[3] = delta*R[1] - Gamma*R[0]   # dIm(rho21)/dt
    return dR

# 存储结果
T = []   # 时间轴
P = []   # rho11 - rho22

# 对四种 Gamma 分别积分
for k, gm in enumerate(Gm):
    gamma1 = gm * omega
    gamma2 = gamma1
    Gamma  = (gamma1 + gamma2)/2 + F
    R0 = [0, 0, 0, 0]      # 初始条件 [rho11, rho22, Re(rho21), Im(rho21)]

    sol = solve_ivp(rhs, t_span, R0, t_eval=t_eval, args=(
                    delta, omega, Lambda1, Lambda2, Gamma, gamma1, gamma2, A))
    t = sol.t
    rho = sol.y                       # 4×N 的数组
    pop_diff = rho[0, :] - rho[1, :]  # rho11 - rho22

    T.append(t)
    P.append(pop_diff)

# ---------------- 绘图 ----------------
plt.figure(figsize=(6, 4))
styles = ['--b', ':r', '-.m', '-k']
labels = [r'$\Gamma=0$', r'$\Gamma=0.2\omega$',
          r'$\Gamma=0.5\omega$', r'$\Gamma=\omega$']

for tk, pk, st, lb in zip(T, P, styles, labels):
    plt.plot(t, p, s, label=lb)

plt.xlabel(r'Time (1/$\gamma$)')
plt.ylabel(r'$\rho_{11}-\rho_{22}$')
plt.legend()
plt.tight_layout()
plt.show()