Python求解复值耦合ODE异常：氨微波激射器径向波近似数值解问题

核心问题 #

1. complex_ode调用完全错误：它的API和odeint不同，不能直接传入初始条件和时间数组，需要通过对象方法分步初始化和求解。
2. 初始条件未声明为复数：微分方程的解是复值，初始条件需显式设为复数类型，否则会被强制转为实数，导致后续计算异常。
3. 频率项过大：代码中2*100的频率太高，在0到3π的时间范围内振荡极快，采样点无法捕捉变化，看起来像是水平线。
4. 布居数计算错误：布居数应该是复振幅的模平方（|A|²、|B|²），而非取实部。

修正后的代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint, complex_ode

# 非旋波近似下的微分方程组
def odes1(x, t, Omega):
    A, B = x
    # 调整频率为Omega，方便观察振荡
    dgp_dt = (1/1j) * 0.5 * B * np.exp(-1j * t * 2*Omega)
    dgm_dt = (1/1j) * 0.5 * A * np.exp(1j * t * 2*Omega)
    return [dgp_dt, dgm_dt]

# 初始条件：显式声明为复数
t0 = np.array([0+0j, 1+0j])
# 时间数组
t = np.linspace(0, 3*np.pi, 1000)
# 设置频率为1，便于观察sin²/cos²曲线
Omega = 1

# 方法1：用odeint直接求解（支持复数）
sol_odeint = odeint(odes1, t0, t, args=(Omega,))
# 计算布居数：模平方
pop_plus = np.abs(sol_odeint[:, 0])**2
pop_minus = np.abs(sol_odeint[:, 1])**2

# 绘图
plt.figure(figsize=(8,4))
plt.plot(t, pop_plus, label=r'$\gamma_{+}$ (No RWA)', color='blue')
plt.plot(t, pop_minus, label=r'$\gamma_{-}$ (No RWA)', color='red')
plt.legend()
plt.xlabel(r"Time ($\frac{1}{\Omega_{0}}$)")
plt.ylabel("Population")
plt.grid(alpha=0.3)
plt.show()

# 方法2：正确使用complex_ode的方式
# 创建complex_ode对象
solver = complex_ode(odes1)
# 设置初始条件和参数
solver.set_initial_value(t0, t[0])
solver.set_f_params(Omega)
# 分步求解
sol_complex = np.zeros((len(t), 2), dtype=np.complex128)
sol_complex[0] = t0
for i in range(1, len(t)):
    sol_complex[i] = solver.integrate(t[i])
# 计算布居数
pop_plus_c = np.abs(sol_complex[:,0])**2
pop_minus_c = np.abs(sol_complex[:,1])**2

# 验证两种方法结果一致
plt.figure(figsize=(8,4))
plt.plot(t, pop_plus_c, label=r'$\gamma_{+}$ (complex_ode)', color='blue', linestyle='--')
plt.plot(t, pop_minus_c, label=r'$\gamma_{-}$ (complex_ode)', color='red', linestyle='--')
plt.legend()
plt.xlabel(r"Time ($\frac{1}{\Omega_{0}}$)")
plt.ylabel("Population")
plt.grid(alpha=0.3)
plt.show()

关键修正点说明 #

• complex_ode正确用法：需先创建solver对象，通过set_initial_value设置初始条件，set_f_params传递额外参数，再循环调用integrate完成求解。
• 复数初始条件：用np.array([0+0j,1+0j])替代实数数组，确保微分方程计算全程保持复值。
• 频率调整：将2*100改为可调节的2*Omega，并设Omega=1，这样在0到3π时间内可以清晰看到布居数在0和1之间以sin²/cos²形式振荡。
• 布居数计算：复振幅的模平方才是实际的粒子布居数，这是量子力学中常见的处理方式，而非直接取实部。

内容的提问来源于stack exchange，提问作者George