如何用Python的ivp_solver事件选项实现双摆碰撞弹跳模拟

问题说明 #

想用Python的solve_ivp实现两个摆锤摆动、碰撞弹跳的模拟，需要通过event选项检测摆锤位置交点（碰撞），中断求解后修改速度，再从碰撞点继续模拟。现有代码仅能找到交点，无法实现中断续算。

现有代码片段 #

变量定义

t_span = [0, 150]
y0_pend = [0, 3]
t = np.linspace(t_span[0], t_span[1], 10000)
p = [2, 3]
b = 0
n = np.linspace(1, 1000, 1000)

函数定义

def pend1_f(t, y):
    k = 1
    d = 0.1

    dx = k * y[1]
    dv = -y[0] - d * dx

    return [dx, dv]

def pend2_f(x, p, b):
    return [b + np.cos(p * x), -p * np.sin(p*x)]

现有ODE求解

zt = solve_ivp(pend_f, t_span, y0_pend, dense_output = True)
z = zt.sol(t)
pendulum_position = z[0]
pendulum_velocity = z[1]

解决方案 #

1. 整合双摆状态向量 #

将两个摆锤的位置、速度合并为一个状态向量y = [x1, v1, x2, v2]，让solve_ivp能同时求解两个摆的运动并检测碰撞。

2. 定义碰撞检测的event函数 #

event函数返回值穿过0时，求解器会触发事件并终止。这里用两个摆的位置差作为检测值：

def collision_event(t, y):
    # y[0] = 摆1位置，y[2] = 摆2位置
    return y[0] - y[2]

collision_event.terminal = True  # 触发事件后终止求解
collision_event.direction = 0   # 检测所有穿越0的情况（正负互转都算）

3. 合并双摆的ODE方程 #

把两个摆的运动方程整合为一个函数，供solve_ivp调用（注意原pend2_f的参数x推测为时间t，若不是请根据物理模型调整）：

def double_pendulum(t, y, p, b, d=0.1, k=1):
    x1, v1, x2, v2 = y
    # 摆1的导数
    dx1 = k * v1
    dv1 = -x1 - d * dx1
    # 摆2的导数
    dx2 = b + np.cos(p[0] * t)
    dv2 = -p[1] * np.sin(p[1] * t)
    return [dx1, dv1, dx2, dv2]

4. 循环求解：碰撞检测→修改速度→续算 #

通过循环反复调用solve_ivp，直到达到总模拟时间：

import numpy as np
from scipy.integrate import solve_ivp

# 初始化参数
t_total = 150
t_current = 0
# 初始状态：摆1[x1, v1]，摆2[x2, v2]
y_current = [0, 3, b + np.cos(p[0]*0), -p[1]*np.sin(p[1]*0)]
# 存储模拟历史数据
t_history = [t_current]
y_history = [y_current.copy()]

# 弹性碰撞恢复系数（e=1为完全弹性碰撞）
e = 1.0

while t_current < t_total:
    # 求解到下一次碰撞或总时间
    sol = solve_ivp(
        double_pendulum,
        [t_current, t_total],
        y_current,
        args=(p, b),
        events=collision_event,
        dense_output=True
    )
    
    # 记录本次求解结果
    t_history.extend(sol.t[1:])
    y_history.extend(sol.y.T[1:])
    
    # 更新当前时间与状态
    t_current = sol.t[-1]
    y_current = sol.y[:, -1].copy()
    
    # 若触发碰撞事件
    if sol.status == 1:
        # 获取碰撞前速度
        v1_before = y_current[1]
        v2_before = y_current[3]
        
        # 质量相同时的弹性碰撞速度公式，质量不同请自行调整
        v1_after = ((1 - e) * v1_before + (1 + e) * v2_before) / 2
        v2_after = ((1 + e) * v1_before + (1 - e) * v2_before) / 2
        
        # 更新状态中的速度
        y_current[1] = v1_after
        y_current[3] = v2_after
        
        # 记录碰撞后的初始状态
        t_history.append(t_current)
        y_history.append(y_current.copy())

# 转换为numpy数组便于后续分析/绘图
t_history = np.array(t_history)
y_history = np.array(y_history)
x1_history = y_history[:, 0]
v1_history = y_history[:, 1]
x2_history = y_history[:, 2]
v2_history = y_history[:, 3]

关键注意点 #

• event属性设置：terminal=True确保碰撞时求解器停止，direction=0避免漏检双向碰撞。
• 碰撞公式调整：如果两个摆锤质量不同，需要替换为对应质量的弹性碰撞速度公式。
• 状态整合必要性：必须将两个摆的状态放在同一向量中，才能让solve_ivp同时跟踪并检测碰撞。

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