【VisionFive 2 Lite 单板计算机】科学计算应用方案：微分方程的数值求解

lijinlei · December 17, 2025, 6:34am

【VisionFive 2 Lite 单板计算机】科学计算解决方案：微分方程组的数值求解

本文介绍了昉·星光 VisionFive2 Lite 单板计算机结合 scipy 实现微分方程组数值求解的项目设计。

项目介绍

VisionFive 2 Lite 结合 scipy 与 matplotlib 软件库实现微分方程组的数值求解。

准备工作：硬件连接、系统安装、软件更新等；
环境搭建：科学计算所需软件库的安装、测试等；
微分方程组：结合实际物理问题进行数值求解和计算，包括流程图、代码、效果演示等；
二能级系统：考虑全量子 JC 模型，数值计算 Rabi 振荡动力学及其参量依赖关系。

准备工作

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

硬件连接

包括供电、联网、显示、鼠标键盘输入等。

系统安装

安装和烧录 VisionFive2 Lite 官方 Ubuntu 操作系统。

详见：VF2 Lite 系统安装。

软件更新

更新软件包

wget https://github.com/starfive-tech/Debian/releases/download/v0.15.0-engineering-release-wayland/install_full.sh
sudo chmod +x install_full.sh
sudo ./install_full.sh

详见： StarFive-Tech | Debian .

环境搭建

安装 scipy 和 matplotlib 软件包，便于调用 RK45 算法和科学绘图；
```
sudo apt install python3-scipy
sudo apt install python3-matplotlib
```

Runge-Kutta 算法

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

详见：Runge-Kutta method .

使用 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 .

数值计算

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

Lorenz 吸引子

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

最早由美国气象学家爱德华·诺顿·洛伦兹在其论文《Deterministic Nonperiodic Flow》中提出。

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

Lorenz 吸引子具有三维、非线性和确定性特征，其相轨线呈现出复杂的双纽线形状，状态随时间以非重复模式演变。

其简化的方程形式为

where σ is Prandtl number, ρ is Rayleigh number, β is system parameter.

对于该常微分方程组，可使用 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  = 1
rho   = 1
x0    = [1.0, 0.0, 0.0]
t_span = (0, 100)                    # 积分区间
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 ode_demo.py 弹窗显示数值结果

使用 VisionFive2 Lite 完全可以胜任常微分方程的求解计算，并获得精确可靠的数值结果。

二能级系统

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

这里以量子光学中常见的二能级系统 Bloch 方程组数值求解为例，在 VisionFive 2 Lite 平台结合 Python 工具包函数实现。

Hamiltonian

考虑半经典（经典光和量子化原子）条件下的光与两能级原子相互作用体系，该系统的示意图如下

探测光耦合原子能级 |1> 和 |2> 之间的跃迁。

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

Bloch 方程组

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

liouville_equation

给出 Hamiltonian 矩阵元对角元的运动方程

two-level_bloch_equations

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

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

代码

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

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

# ---------- 参数 ----------
Delta   = np.arange(-1, 1, 0.01) # range  -10:.1:10
Omega   = 0.1
gamma21 = 1.0
t_span  = (0, 10)                   # 积分区间
R0      = np.array([0.0, 0.0])      # 初值 ρ=[ρ1, ρ2]
Nd      = len(Delta)                # number of delta
# ---------- ODE ----------
def rk(t, R, Delta, Omega, gamma21):
    r1, r2 = R
    dr1 = -(gamma21*r2 + Delta*r1)
    dr2 = -(gamma21*r1 - Delta*r2) 
    return np.array([dr1, dr2])

# ---------- 扫频 ----------
r21R = np.zeros(Nd)   # 实部 χ′
r21I = np.zeros(Nd)   # 虚部 χ″

print("Start calculation...")
for m, d in enumerate(Delta):
    sol = solve_ivp(rk, t_span, R0, args=(d, Omega, gamma21),
                    method='RK45', rtol=1e-6, atol=1e-6)
    # 取终值
    r21R[m] = sol.y[0, -1]
    r21I[m] = sol.y[1, -1]
    # 显示计算进度
    if (m+1) % 20 == 0:
        print(f"Process: {m+1}/{Nd} (Delta = {Delta[m]:.1f})")
print("Calculation complete!")

# ---------- Plot ----------
print("Start Drawing...")
plt.plot(Delta, r21R, label=r"$\chi'$", linewidth=4)
plt.plot(Delta, r21I, label=r"$\chi''$", linewidth=4)
plt.xlabel(r'$\Delta$ ($\gamma$)')
plt.ylabel(r'$\chi$')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()
print("Terminate Program.")

保存代码。

效果

终端运行 python sci_two-level_system.py ，输出计算进度

弹窗显示绘图结果

获得介质极化率的实部（色散）和虚部（吸收）对应的光谱。在探测光共振时吸收最强，色散最小。

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

可知，矩阵对角元 ρ21 的时间演化快速趋于稳定。

全量子模型

在全量子（光和原子均量子化）领域，基于二能级系统的 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

# ---------- Parameters ----------
omega1 = 0
Lambda1 = 1
Lambda2 = 1
gamma1 = 1
gamma2 = 1
F = 1
Gamma = (gamma1 + gamma2) / 2 + F
A = 1
Dn = np.array([0, 1, 2, 3])          # n = 0,1,2,3
t_eval = np.arange(0, 80.1, 0.1)     # time range

# ---------- ODE defination ----------
def rk(t, R, delta, omega1, Lambda1, Lambda2, Gamma, gamma1, gamma2, A):
    r11, r22, re21, im21 = R
    dr11 = Lambda1 + gamma1 * r11 + A * r22 + 2 * im21
    dr22 = Lambda2 + gamma2 * r22 - 2 * omega1 
    dre21 = -Gamma * re21 - delta * im21
    dim21 = delta * re21 - Gamma * im21 - omega1 * (r22 - r11) / 2
    return np.array([dr11, dr22, dre21, dim21])

# ---------- Scan N ----------
P = np.empty((len(t_eval), 4))   # Population ρ₁₁−ρ₂₂
for k, n in enumerate(Dn):
    delta = n * omega1
    R0 = np.array([1.0, 0.0, 0.0, 0.0])  # [ρ11, ρ22, Reρ21, Imρ21]
    sol = solve_ivp(rk, [0, 80], R0, args=(delta, omega1, Lambda1, Lambda2,
                                           Gamma, gamma1, gamma2, A),
                    t_eval=t_eval, rtol=1e-6, atol=1e-6)
    P[:, k] = sol.y[0, :] - sol.y[1, :]  # ρ₁₁−ρ₂₂

print("Calculation complete!")

# ---------- Draw ----------
print("Start Drawing...")
plt.plot(t_eval, P[:, 0], '-k',  label=r'$\delta=0$', linewidth=3)
plt.plot(t_eval, P[:, 1], '--r', label=r'$\delta=0.2\gamma$', linewidth=3)
plt.plot(t_eval, P[:, 2], '-.b', label=r'$\delta=0.4\gamma$', linewidth=3)
plt.plot(t_eval, P[:, 3], ':m',  label=r'$\delta=0.6\gamma$', linewidth=3)
plt.xlabel(r'Time ($1/\gamma$)')
plt.ylabel(r'$\rho_{11}-\rho_{22}$')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()
print("Terminate Program.")

保存代码。

效果

终端运行 python jc_p-t.py ，弹窗显示绘图结果；

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

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

给出对应的三维渲染结果

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

总结

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