【“星睿O6”AI PC开发套件评测】科学计算

本文介绍了瑞莎星睿 O6 (Radxa Orion O6) 开发板结合 scipy 实现微分方程组的数值求解的项目设计，结合具体实际物理问题，给出相应的数值解决方案。

项目介绍

硬件连接：联网、显示、供电、SSH 登录等；
环境搭建：科学计算所需软件库的安装、测试等；
数值计算：结合实际物理问题进行数值求解和计算，包括流程图、代码、效果演示等；

硬件连接

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

环境搭建

为了避免影响系统 Python，创建并激活虚拟环境；
安装 scipy 和 matplotlib 软件包，便于调用 RK45 算法和科学绘图；

mkdir ~/SCI          # 创建文件夹
cd ~/SCI             # 进入文件夹
python3 -m venv venv # 新建虚拟环境
source venv/bin/activate  # 激活环境
python -m pip install -U pip # 升级 pip（避免旧 pip 装包失败）
python -m pip install numpy matplotlib scipy  # 安装科学计算三件套

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

# -------- parameters ---------
sigma = 1
beta  = 1
rho   = 1.0
x0    = [1.0, 0.0, 0.0]
t_span = (0, 100)                    # time span
t_eval = np.arange(0, 100, 0.01)  # fixed step 0.01 output

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

# ---------- Draw -----------
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_delta.py 新建程序文件，并添加如下代码

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

# ---------- parameters ----------
Delta   = np.arange(-5, 5, 0.05)    # delta range
Omega   = 0.1
gamma21 = 0.7
t_span  = (0, 10)                   # integral range
R0      = np.array([0.0, 0.0])      # initialize 
Nd      = len(Delta)                # number of delta
# ---------- ODE ----------
def rk(t, R, Delta, Omega, gamma21):
    r1, r2 = R
    dr1 = -(gamma2*r1 - Delta*r1)
    dr2 = -(gamma1*r2 + Delta*r2) 
    return np.array([dr1, dr2])

# ---------- scan delta ----------
r21R = np.zeros(Nd)   # real 
r21I = np.zeros(Nd)   # image 

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)
    # get terminal value
    r21R[m] = sol.y[0, -1]
    r21I[m] = sol.y[1, -1]
    # show progress
    if (m+1) % 10 == 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("Program Terminate.")

保存代码。