Python中使用odeint求解含线性方程的常微分方程组

你需要求解的是含线性代数约束的常微分方程组，其中u(t)不是独立的外部输入，是由状态x(t)、y(t)直接决定的代数变量，不需要单独迭代更新。

原代码的错误点 #

• 把u错误当成了分段常数的外部输入，每个积分区间内传入固定的u值，没有在积分过程中根据当前时刻的x、y实时计算u，完全破坏了u = x + 3*y的代数约束关系。
• 逐段循环调用odeint属于冗余实现，对于光滑的常微分系统，一次性传入全时间区间求解精度更高、运行效率更高。
• 初始条件存在矛盾：u作为非独立的代数变量没有独立初始值，当x(0)=0、y(0)=0时，按代数关系u(0)必然为0，和你给出的u(0)=1冲突。如果需要满足u(0)=1，需要调整x、y的初始值满足x(0)+3*y(0)=1。

对于这类可以直接消去代数变量的系统，最简单的方法是把u的表达式直接代入微分方程，将系统转化为纯常微分方程后直接调用odeint求解，不需要额外循环：

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

def ode_system(z, t):
    x, y = z
    # 积分过程中直接根据当前状态计算u，不需要外部传入
    u = x + 3 * y
    dxdt = (-x + u) / 2.0
    dydt = (-y + x) / 5.0
    return [dxdt, dydt]

# 初始条件：若保持x(0)=0,y(0)=0则u(0)自然为0，若需要u(0)=1可调整为比如z0=[1,0]
z0 = [0, 0]
t = np.linspace(0, 15, 150)
# 一次性求解全时段状态
z_sol = odeint(ode_system, z0, t)
x = z_sol[:, 0]
y = z_sol[:, 1]
u = x + 3 * y # 求解完成后直接计算全时段u值

# 绘图
plt.plot(t, u, label='u(t)')
plt.plot(t, x, label='x(t)')
plt.plot(t, y, label='y(t)')
plt.xlabel('t')
plt.legend()
plt.show()

• 对于无法直接消去代数变量的复杂微分代数方程组（DAE），可以使用scipy.integrate.solve_ivp通过定义质量矩阵的方式求解，也可以使用专门的DAE求解库。你这个案例属于最简单的线性DAE，直接消元是效率最高的方案。
• 该系统是线性时不变系统，可以写成状态空间形式$\dot{\boldsymbol{z}} = \boldsymbol{A}\boldsymbol{z}$，通过计算矩阵指数得到解析解，结果和数值积分完全一致。

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