使用 odeint 进行刚体模拟的数值积分

Question

我正在尝试对以下运动方程进行数值积分并求解 omega 矢量 (3x1)：

所以上面的等式是我想对给定的初始Omega_naught向量、惯性矩阵（=单位矩阵）和力矩向量（=零向量）进行数值积分。

目前我正在尝试使用 scipy 中的 odeint 但它抛出一个 ValueError: Initial condition y0 must be one-dimensional .

这是我的方法

I = np.array([[10, 0, 0], [0, 15, 0], [0, 0, 20]]) 
C0 = np.eye(3)
M = np.zeros(3).reshape(3, 1)
I_inv = np.linalg.inv(I)

def skew(x):
    return np.array([[0, -x[2], x[1]],
                     [x[2], 0, -x[0]],
                     [-x[1], x[0], 0]])
def model(w, t):
    H = np.matmul(I, w) #angular momentum
    A = M - np.matmul(skew(w), H)
    
    dwdt = np.matmul(I_inv, A)
    
    return dwdt

#Initial condition
w0 = np.array([0.01, 10, 0]).reshape(3, 1)

#time points
t = np.linspace(0, 20)

#solve ode
w = odeint(model, w0, t)

我从来没有在矩阵方程中使用过 odeint，所以我不确定我是否对方程使用了正确的积分方法。我如何使用 odeint 解决这个问题，还是应该使用不同的集成方法？

我也愿意接受 MATLAB 提示或答案。

符号：

A - 3x1 向量
[I] - 3x3 矩阵
tilde_A - 偏斜对称矩阵

Answer 1

错误的意思是初始点应该只有一维，即应该是一个扁平的numpy数组。 ODE函数内部也一样，接口是平面数组，你必须手动建立和破坏列向量的结构。但这似乎会引起奇怪的后续错误，所以走另一条路，让一切都变得无形。规则是矩阵乘以无形数组 returns 无形数组。不混，那请“广播”到一个矩阵。

M = np.zeros(3)

def model(w, t):
    H = np.matmul(I, w) #angular momentum
    A = M - np.matmul(skew(w), H)
    
    dwdt = np.matmul(I_inv, A)
    
    return dwdt

#Initial condition
w0 = np.array([0.01, 10, 0])

经过这些更改，代码对我有用。

使用 odeint 进行刚体模拟的数值积分

Numerical integration for simulation of rigid body using odeint

python

matrix

scipy

ode

odeint