EEPW论坛

【树莓派5】科学计算应用：微分方程组的求解

本文介绍了树莓派 5 结合 Python 软件包 scipy 实现常微分方程组的数值求解，给出该设备在科学计算领域的应用解决方案。

项目介绍

树莓派 5 部署 Pyhton 科学计算相关软件包，实现数值求解常微分方程组。

• 准备工作：Python 安装、scipy 和 matplotlib 等软件包的安装、单个微分方程求解测试等；

• 微分方程组的求解：结合 solve_ivp 函数及 Runge-Kutta 算法实现半经典三能级系统的 Bloch 方程组数值求解，并给出介质极化率与失谐的函数关系；

• 二能级原子系统动力学演化：数值求解全量子 Bloch 方程组，给出 Rabi 振荡的三维动力学演化结果；

准备工作

• 树莓派 5 安装 Raspberry Pi OS 官方操作系统；

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

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

测试

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

代码

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

 from scipy.integrate import odeint
 import numpy as np
 import matplotlib.pyplot as plt
 #定义ODE方程
 def model(y, t):
     dydt = -2 * y + t
     return dydt
 # 初始条件
 y0 = 1
 # 时间点
 t = np.linspace(0, 5, 100)
 # 使用odeint求解ODE
 y = odeint(model, y0, t)
 # 绘制结果
 plt.plot(t, y, 'ro-')
 plt.xlabel('time')
 plt.ylabel('y(t)')
 plt.show()

效果

终端运行 python sci_demo.py 弹窗显示数值结果

MATLAB 对比

使用 MATLAB 的 ode45 函数验证上述结果

 ode = @(t, y) -2*y + t;
 [t, y] = ode45(ode, [0, 5], 1);
 
 % 显示结果
 plot(t, y, 'b-o');
 xlabel('t'); ylabel('y');
 title('ode45求解: dy/dt = -2y + t');
 grid on;

运行代码后输出数值求解曲线

对比 MATLAB 和 Python 的数值求解结果，使用树莓派 5 完全可以胜任常微分方程的求解计算，并获得精确可靠的数值结果。

微分方程组

在前面测试常微分方程数值求解的基础上，进一步考虑更为复杂的实际科学计算。

这里以量子光学中常见的 Lambda 型三能级系统 Bloch 方程组数值求解为例，在树莓派 5 平台结合 Python 工具包函数实现。

Hamiltonian

考虑半经典（经典光和量子化原子）条件下的光与三能级原子相互作用体系。

该系统的示意图如下

信号光耦合原子能级 |1> 和 |3> 之间的跃迁，控制光耦合|2> 和 |3> 之间的跃迁。

该光与原子耦合系统 Hamiltonian 可表示为

Bloch 方程组

结合 Maxwell-Liouville 方程（亦称光学 Bloch 方程或主方程），

给出 Hamiltonian 各个矩阵元的运动方程

该方程组可采用数值或解析的方法求解。

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

Runge-Kutta 算法

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

详见：Runge-Kutta method .

• 对于 MATLAB 编程，通常采用 ode45 函数实现常微分方程组的数值求解；

• 对于 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 | scipy .

代码

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

 import numpy as np
 from scipy.integrate import solve_ivp
 import matplotlib.pyplot as plt
 
 def rk_ode(t, R, Deltas, Deltac, Omegas, Omegac, gamma21, gamma31):
     """定义微分方程系统"""
     dR = np.zeros(4)
 
     dR[0] = -(gamma21 * R[0] + (Deltas - Deltac) * R[1]) - Omegac * R[3]
     dR[1] = -(gamma21 * R[1] - (Deltas - Deltac) * R[0]) + Omegac * R[2]
     dR[2] = -(gamma31 * R[2] + Deltas * R[3]) - Omegac * R[1]
     dR[3] = -(gamma31 * R[3] - Deltas * R[2]) + Omegac * R[0] + Omegas
 
     return dR
 
 def main():
     # 参数设置
     Hp = 0.1        # step
     Deltas = np.arange(-10, 10 + Hp, Hp)  # range [-10, 10]
     Nd = len(Deltas) # number of delta_s
 
     Deltac = 0      # delta of control
     Omegac = 3      # Rabi frequency of control
     Omegas = 0.01   # Rabi frequency of signal
     gamma21 = 0.001 # decay of transition between |2> and |1> 
     gamma31 = 3     # decay of transition between |3> and |1> 
 
     # initial condition [rho12_re, rho12_im, rho13_re, rho13_im]
     R0 = np.array([0.0, 0.0, 0.0, 0.0])
 
     # 时间区间
     tspan = [0, 10]
 
     # 保存演化结果
     r31R = np.zeros(Nd)
     r31I = np.zeros(Nd)
 
     # 遍历 delta_s
     print("开始计算...")
     for m in range(Nd):
         # 使用 solve_ivp 求解 ODE（对应 MATLAB 的 ode45）
         sol = solve_ivp(
             lambda t, R: rk_ode(t, R, Deltas[m], Deltac, Omegas, Omegac, gamma21, gamma31),
             tspan,
             R0,
             method='RK45',
             dense_output=True,
             rtol=1e-6,
             atol=1e-8
         )
 
         # 获取演化最终时刻的值
         Nt = len(sol.t)
         r31R[m] = sol.y[2, -1]  # Rho(Nt,3)
         r31I[m] = sol.y[3, -1]  # Rho(Nt,4)
 
         # 显示计算进度
         if (m+1) % 20 == 0:
             print(f"进度: {m+1}/{Nd} (Delta = {Deltas[m]:.1f})")
 
     print("Calculation Complete! ")
 
     # Draw figure
     plt.figure(figsize=(5, 4))
 
     plt.plot(Deltas, r31R, 'b-', linewidth=2, label='Re[rho31]')
     plt.plot(Deltas, r31I, 'r-', linewidth=2, label='Im[rho31]')
 
     plt.xlabel('delta_s', fontsize=12)
     plt.ylabel('rho_31', fontsize=12)
     plt.title('EIT', fontsize=14)
     plt.grid(True, alpha=0.3)
     plt.legend(fontsize=10)
     plt.tight_layout()
 
     # Show figure
     try:
         plt.show()
     except KeyboardInterrupt:
         print("\n\n程序被用户中断")
     except Exception as e:
         print(f"\n程序运行出错: {e}")
     finally:
         print("\n程序结束")
 
 if __name__ == "__main__":
     main()

保存代码。

效果

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

弹窗显示绘图结果

相应的矩阵对角元随时间的演化曲线

由此可知，当信号光失谐一定时，矩阵对角元（光与原子的耦合项）的演化最终趋于稳定。

二能级系统

在全量子（光和原子均量子化）领域，基于二能级系统的 Jaynes-Cummings 模型扮演着基础且重要的角色。

Bloch 方程

光与二能级原子相互作用的示意图如下

二能级系统的 Hamiltonian 可表示为

结合光学 Bloch 方程，可得

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

代码

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

 import numpy as np
 from scipy.integrate import solve_ivp
 import matplotlib.pyplot as plt
 
 # 定义ODE系统
 def rk_ode(t, R, delta, omega1):
     dR = np.zeros(3)
     dR[0] = 2 * omega1 * R[2]
     dR[1] = -delta * R[2]
     dR[2] = delta * R[1] - omega1 * R[0] / 2
     return dR
 
 # 设置参数
 omega1 = 0.2
 D = np.arange(0, 5.1, 0.1)  # 0:.1:5
 num_D = len(D)
 tspan = np.arange(0, 101)  # 0:100
 num_t = len(tspan)
 
 # 预分配数组
 P = np.zeros((num_t, num_D))           # ρ22-ρ11
 sigma21_Re = np.zeros((num_t, num_D))  # Re[ρ21]
 sigma21_Im = np.zeros((num_t, num_D))  # Im[ρ21]
 Dta = np.zeros((num_t, num_D))         # delta矩阵
 T = np.zeros((num_t, num_D))           # 时间矩阵
 
 # 循环求解
 print("正在计算...")
 for k in range(num_D):
     delta = D[k] * omega1
     R0 = np.array([-1.0, 0.0, 0.0])  # ρ22-ρ11, Re[ρ21], Im[ρ21]
 
     # 求解 ODE（使用 RK45 方法）
     sol = solve_ivp(
         lambda t, R: rk_ode(t, R, delta, omega1),
         [tspan[0], tspan[-1]],
         R0,
         method='RK45',
         t_eval=tspan,
         rtol=1e-6,
         atol=1e-8
     )
 
     # 存储结果
     T[:, k] = sol.t
     P[:, k] = sol.y[0, :]           # ρ22-ρ11
     sigma21_Re[:, k] = sol.y[1, :]  # Re[ρ21]
     sigma21_Im[:, k] = sol.y[2, :]  # Im[ρ21]
     Dta[:, k] = delta
 
     # 显示进度
     if (k+1) % 10 == 0:
         print(f"进度: {k+1}/{num_D}")
 
 print("计算完成!")
 
 # 创建网格用于surf绘图
 T_mesh, Dta_mesh = np.meshgrid(tspan, Dta[0, :]/omega1, indexing='ij')
 
 # 绘图1：俯视图
 plt.figure(figsize=(5, 4))
 plt.contourf(T_mesh, Dta_mesh, P, 50, cmap='viridis')
 plt.colorbar(label='ρ₁₁-ρ₂₂')
 plt.xlabel('Time (1/γ)', fontsize=12)
 plt.ylabel('Δ (γ)', fontsize=12)
 plt.title('Population Difference ρ₁₁-ρ₂₂', fontsize=14)
 plt.tight_layout()
 plt.show()
 
 # 绘图2：3D曲面图
 fig = plt.figure(figsize=(5, 4))
 ax = fig.add_subplot(111, projection='3d')
 surf = ax.plot_surface(T_mesh, Dta_mesh, P, cmap='viridis', alpha=0.8)
 plt.colorbar(surf, ax=ax, label='ρ₁₁-ρ₂₂')
 ax.set_xlabel('Time (1/γ)', fontsize=12)
 ax.set_ylabel('Δ (γ)', fontsize=12)
 ax.set_zlabel('ρ₁₁-ρ₂₂', fontsize=12)
 ax.set_title('Population Difference ρ₁₁-ρ₂₂', fontsize=14)
 plt.show()

保存代码。

效果

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

弹窗显示绘图结果

• 不同失谐（Delta）条件下，粒子布居振荡（ρ11-ρ22）随时间的演化如下

• 给出布居振荡随时间和失谐变化的伪彩图

• 给出对应的三维渲染结果

该结果表明原子在基态和激发态之间的 Rabi 振荡 rho11 - rho22 随失谐的增加而减弱。

总结

本文介绍了树莓派 5 结合 Python 软件包 scipy 实现常微分方程组的数值求解，给出了该设备在科学计算领域的应用解决方案，包括半经典和全量子体系的 Bloch 方程组的数值求解、原理、方程、代码与结果，为相关产品在工业、科研、科学计算领域的快速应用提供了参考。