使用 odeint 模拟刚体的数值积分答案

【问题标题】：Numerical integration for simulation of rigid body using odeint使用 odeint 模拟刚体的数值积分
【发布时间】：2020-12-01 05:44:51
【问题描述】：

我正在尝试对以下运动方程进行数值积分并求解欧米茄矢量（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 矩阵
波浪号_A - 斜对称矩阵

【问题讨论】：

标签： python matrix scipy ode odeint

【解决方案1】：

错误的意思是初始点应该只有一维，也就是应该是一个扁平的numpy数组。 ODE函数内部也是如此，接口是平面数组，您必须手动建立和销毁列向量的结构。但这似乎会引发奇怪的后续错误，所以反过来，让一切变得无形。规则是与无形数组的矩阵乘法返回无形数组。不要混合，这会邀请“广播”到矩阵。

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])

通过这些更改，代码适用于我。

【讨论】：