odeint求解含离散函数的ODE结果异常，请求排查修正

首先，我对比了你的Python代码和Modelica模型的数学表达式，确认两者的核心方程是完全一致的：

从Modelica的方程：

m * a + b * v + k * x + M * sign(v) = Fmax;

整理得到加速度的表达式：
$$a = \frac{F_{max} - bv - kx - M \cdot sign(v)}{m}$$

而你的Python代码中dXdt函数返回的第二个元素（加速度）：

- b * X[1] / m - k * X[0] / m - M * np.sign(X[1]) / m + Fmax / m

和上面的数学表达式完全等价，所以模型本身没有错误。

结果不一致的原因 #

问题出在不连续函数的积分处理上：

• 你的模型中包含sign(v)这个不连续函数，当速度v=0时，函数值会从-1突变到1（或者反过来）。
• Python的odeint（基于LSODA积分器）默认设置下，对这种不连续点的检测能力有限，积分过程中可能会跳过这些突变点，导致误差累积，最终和Modelica的结果产生偏差。
• 而OpenModelica这类Modelica仿真器会自动识别不连续事件（比如v=0就是一个事件触发条件），在事件发生时会暂停积分，更新系统状态，然后重新启动积分，这样能更精准地处理不连续带来的状态变化。

解决方法：改用solve_ivp并添加事件检测 #

Scipy的solve_ivp函数比odeint更适合处理带不连续项的ODE，它支持自定义事件函数来检测v=0的时刻，让积分器在这些点处停止并重新初始化，从而得到和Modelica一致的结果。

修改后的Python代码如下：

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

m = 1
k = 1
M = 0.1
b = 1
Fmax = 1

def dXdt(t, X):
    # 注意solve_ivp的函数签名是(t, X)，和odeint的(X,t)相反
    return [X[1], (Fmax - b * X[1] - k * X[0] - M * np.sign(X[1])) / m]

# 定义事件函数：当v=0时触发事件，设置terminate=False表示不终止积分，仅记录事件点
def event(t, X):
    return X[1]  # 速度v=X[1]，当返回值为0时触发事件
event.terminal = False
event.direction = 0  # 检测所有穿越0的情况（从正到负或负到正）

X0 = [1, 2]
t_span = [0, 10]
t_eval = np.linspace(0, 10, 200)

# 调用solve_ivp，添加事件参数
sol = solve_ivp(dXdt, t_span, X0, t_eval=t_eval, events=event, rtol=1e-6, atol=1e-9)

plt.plot(sol.t, sol.y[0])
plt.xlabel('Time')
plt.ylabel('x(t)')
plt.show()

额外说明 #

• 如果你坚持要用odeint，可以尝试增大mxstep参数（比如mxstep=10000），让积分器更频繁地检查状态，一定程度上改善结果，但这种方法不如solve_ivp的事件处理可靠。
• 另外，np.sign(0)返回0，和Modelica的sign(0)行为一致，所以这部分没有差异。

内容的提问来源于stack exchange，提问作者Foad S. Farimani