手动实现函数与基于odeint的矩阵乘法微分方程求解结果差异排查

首先要明确：你看到的结果差异，核心原因是原函数f的符号错误，导致两个函数求解的是完全不同的微分方程系统，而你的f_new函数其实是正确实现了手动定义的阻尼振荡系统。

1. 原函数f的符号错误 #

你手动定义的微分方程是：

def f(x, t): return [ -0.1 * x[0] + 2 * x[1], -2 * x[0] - 0.1 * x[1] ]

对应的系统矩阵是：
$$
A = \begin{bmatrix} -0.1 & 2 \\ -2 & -0.1 \end{bmatrix}
$$
这是一个阻尼振荡系统，特征值为$-0.1 \pm 2j$，解会随时间指数衰减到0。

但你适配矩阵的f函数第二个式子用了错误的符号：

def f(x, t, a): return [ a[0] * x[0] + a[1] * x[1], a[2] * x[0] - a[3] * x[1] ]

结合你的matrix_a = np.array([-0.09999975, 1.999999, -1.999999, -0.09999974])，第二个式子实际变成：
$$-1.999999x_0 - (-0.09999974)x_1 = -1.999999x_0 + 0.09999974x_1$$
对应的系统矩阵是：
$$
A_{\text{wrong}} = \begin{bmatrix} -0.09999975 & 1.999999 \\ -1.999999 & 0.09999974 \end{bmatrix}
$$
这个矩阵的迹几乎为0，行列式为正，特征值是纯虚数，系统是无阻尼振荡，解不会随时间衰减，所以t=1000时仍有明显振幅——这和你手动定义的系统完全不是一回事。

2. f_new函数的正确性与优化点 #

你的f_new函数是正确的：从new_matrix_a的第二行非零元素可以看出，第二个微分方程是$-1.999999x_0 -0.09999974x_1$，和手动定义的系统完全一致，所以t=1000时解衰减到几乎为0是符合预期的正确结果。

不过f_new有几个可以优化的地方，提升效率和可读性：

• 避免重复创建PolynomialFeatures实例：每次调用f_new都重新初始化并拟合特征转换器完全没必要，因为输入特征的结构（2维，degree=5）是固定的，应该在函数外提前创建好：

  from sklearn.preprocessing import PolynomialFeatures
  import numpy as np
  from scipy.integrate import odeint
  
  # 提前创建并拟合特征转换器（用任意2维样本即可，只是确定特征结构）
  polynomials = PolynomialFeatures(degree=5)
  polynomials.fit(np.array([[0, 0]]))
  
  def f_new(x, t, parameters):
      x = np.array(x).reshape(1, -1)
      # 直接转换特征，不需要重复拟合
      polynomial_features = polynomials.transform(x)[0]  # 得到(21,)的特征数组
      # 用点积替代矩阵乘法，更简洁高效
      x_ode = parameters[0] @ polynomial_features
      y_ode = parameters[1] @ polynomial_features
      return [x_ode, y_ode]
  

• 简化维度处理：不需要对特征数组转置，直接用一维数组的点积就能得到标量结果，代码更简洁。

总结 #

你看到的结果差异，本质是原函数f的符号错误导致求解了错误的微分方程系统。f_new的实现是正确的，完全匹配你手动定义的阻尼振荡系统，末尾结果趋近于0是符合物理规律的正确表现。

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