基于scipy.integrate.solve_ivp求解变刚度二自由度弹簧质量系统

你定义的状态向量y的第0位就是质量块m1的位移x1，可直接读取该值做条件判断，另外注意你预设的逻辑里刚度和返回值的对应关系写反了：规则为x1 <= -1时k1=14，对应你代码里的ydot2，x1 > -1时k1=7对应ydot1。

第一步：替换return段代码 #

将代码中return ???????的部分替换为以下内容：

# 读取m1的位移值
x1 = y[0]
if x1 <= -1:
    # 对应k1=14的刚度规则，返回ydot2，展平为一维数组适配solve_ivp要求
    return ydot2.flatten()
else:
    # 对应k1=7的刚度规则，返回ydot1，展平为一维数组适配solve_ivp要求
    return ydot1.flatten()

加.flatten()是因为你计算得到的ydot1/ydot2是二维列向量，而solve_ivp要求微分方程函数返回一维数组，否则会触发维度不匹配报错。

第二步（可选，提升计算精度） #

刚度突变属于常微分方程的不连续点，自适应步长求解器直接跨点计算可能出现精度损失，你可以给solve_ivp添加事件函数，让求解器在x1=-1的位置自动检测并调整步长：
首先在微分方程函数外定义事件函数：

def x1_cross_threshold(t, y):
    # 事件触发条件：y[0] + 1 = 0，即x1=-1
    return y[0] + 1
# 设置事件触发后不需要终止求解
x1_cross_threshold.terminal = False
# 检测两个方向的穿越（x1从大于-1变为小于等于-1，或者反向穿越）
x1_cross_threshold.direction = 0

调用solve_ivp时新增events参数即可：

sol = solve_ivp(F, time_interval, initial_conditions, t_eval = TE, vectorized=True, method = 'RK45', events=x1_cross_threshold)

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