使用SymPy Mechanics模块时遇‘无法解包非可迭代Particle对象’错误

问题现象 #

使用SymPy Mechanics模块构建三重摆模型时，通过循环向particles列表添加Particle对象，调用KM.kanes_equations(forces, particles)时触发错误：cannot unpack non-iterable Particle object，但forces列表采用相同构建逻辑却无问题。

错误原因 #

该错误核心是**kanes_equations方法对传入的质点参数格式/顺序不符合要求**：

• 若使用较旧版本SymPy，kanes_equations的第二个参数不接受Particle对象，而是需要(质点位置点, 质量)的元组序列；
• 若参数顺序写反（将质点列表作为第一个参数传入），方法会尝试将Particle对象解包为力的元组(Point, Force)，导致解包失败。

修复方案 #

针对两种情况，分别提供修复方式：

方案1：适配旧版SymPy（不使用Particle对象） #

修改循环中质点的存储逻辑，将Particle对象替换为(位置点, 质量)元组：

# 替换循环内的质点创建代码
# Pai = mechanics.Particle('Pa' + str(i), Pi, m[i])
# particles.append(Pai)
particles.append((Pi, m[i]))

方案2：确保参数顺序正确（新版SymPy） #

若使用新版SymPy（支持传入Particle对象），检查kanes_equations的参数顺序，确保第一个参数是力的列表，第二个是质点/刚体列表（若之前顺序写反则交换）：

# 确保参数顺序：(力列表, 质点/刚体列表)
fr, fr_star = KM.kanes_equations(forces, particles)

注：新版SymPy官方文档中，kanes_equations的标准调用格式为kanes_equations(loads, bodies)，其中loads是力的元组序列，bodies是Particle或RigidBody对象序列。

完整修正代码（方案1适配版） #

def integrate_pendulum(n, times,
                       initial_positions=135,
                       initial_velocities=0,
                       lengths=None, masses=1):
    """Integrate a multi-pendulum with `n` sections"""
    #-------------------------------------------------
    # Step 1: construct the pendulum model
    
    # Generalized coordinates and velocities
    # (in this case, angular positions & velocities of each mass) 
    q = mechanics.dynamicsymbols('q:{0}'.format(n))
    u = mechanics.dynamicsymbols('u:{0}'.format(n))

    # mass and length
    m = symbols('m:{0}'.format(n))
    l = symbols('l:{0}'.format(n))

    # gravity and time symbols
    g, t = symbols('g,t')
    
    #--------------------------------------------------
    # Step 2: build the model using Kane's Method

    # Create pivot point reference frame
    A = mechanics.ReferenceFrame('A')
    P = mechanics.Point('P')
    P.set_vel(A, 0)

    # lists to hold particles, forces, and kinetic ODEs
    # for each pendulum in the chain
    particles = []
    forces = []
    kinetic_odes = []

    for i in range(n):
        # Create a reference frame following the i^th mass
        Ai = A.orientnew('A' + str(i), 'Axis', [q[i], A.z])
        Ai.set_ang_vel(A, u[i] * A.z)

        # Create a point in this reference frame
        Pi = P.locatenew('P' + str(i), l[i] * Ai.x)
        Pi.v2pt_theory(P, A, Ai)

        # 存储(位置点, 质量)元组而非Particle对象
        particles.append((Pi, m[i]))

        # Set forces & compute kinematic ODE
        forces.append((Pi, m[i] * g * A.x))
        kinetic_odes.append(q[i].diff(t) - u[i])

        P = Pi

    # Generate equations of motion
    KM = mechanics.KanesMethod(A, q_ind=q, u_ind=u,
                               kd_eqs=kinetic_odes)
    fr, fr_star = KM.kanes_equations(forces, particles)
    
    #-----------------------------------------------------
    # Step 3: numerically evaluate equations and integrate

    # initial positions and velocities – assumed to be given in degrees
    y0 = np.deg2rad(np.concatenate([np.broadcast_to(initial_positions, n),
                                    np.broadcast_to(initial_velocities, n)]))
        
    # lengths and masses
    if lengths is None:
        lengths = np.ones(n) / n
    lengths = np.broadcast_to(lengths, n)
    masses = np.broadcast_to(masses, n)

    # Fixed parameters: gravitational constant, lengths, and masses
    parameters = [g] + list(l) + list(m)
    parameter_vals = [9.81] + list(lengths) + list(masses)

    # define symbols for unknown parameters
    unknowns = [Dummy() for i in q + u]
    unknown_dict = dict(zip(q + u, unknowns))
    kds = KM.kindiffdict()

    # substitute unknown symbols for qdot terms
    mm_sym = KM.mass_matrix_full.subs(kds).subs(unknown_dict)
    fo_sym = KM.forcing_full.subs(kds).subs(unknown_dict)

    # create functions for numerical calculation 
    mm_func = lambdify(unknowns + parameters, mm_sym)
    fo_func = lambdify(unknowns + parameters, fo_sym)

    # function which computes the derivatives of parameters
    def gradient(y, t, args):
        vals = np.concatenate((y, args))
        sol = np.linalg.solve(mm_func(*vals), fo_func(*vals))
        return np.array(sol).T[0]

    # ODE integration
    return odeint(gradient, y0, times, args=(parameter_vals,))

验证步骤 #

运行以下代码测试修复后的模型：

import numpy as np
import matplotlib.pyplot as plt
from sympy import symbols, Dummy
from sympy.physics import mechanics
from scipy.integrate import odeint

def get_xy_coords(p, lengths=None):
    """Convert pendulum simulation results to x/y coordinates"""
    n = p.shape[1] // 2
    if lengths is None:
        lengths = np.ones(n) / n
    x = np.zeros((p.shape[0], n))
    y = np.zeros((p.shape[0], n))
    for i in range(n):
        x[:, i] = np.sum(lengths[:i+1] * np.sin(p[:, :i+1]), axis=1)
        y[:, i] = -np.sum(lengths[:i+1] * np.cos(p[:, :i+1]), axis=1)
    return x, y

t = np.linspace(0, 10, 1000)
p = integrate_pendulum(n=3, times=t)
x, y = get_xy_coords(p)
plt.plot(x, y);
plt.show()

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