Numerical integration for simulation of rigid body using odeint

Question

I am trying to numerically integrate the following equation of motion and solve for omega vector(3x1) :

So the above equation is I want to numerically integrate given initial Omega_naught vector, inertia matrix(=identity matrix) and moment vector(=zero vector).

At the moment I am trying to use odeint from scipy but it throws me a ValueError: Initial condition y0 must be one-dimensional.

Here is my approach

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)

I have never used odeint with matrix equations so I am not sure if I am using the right integration method for the equation. How can I resolve this using odeint or should I use a different integration method?

I am also open to MATLAB hints or answers.

Notations:

1. A - 3x1 vector
2. [I] - 3x3 matrix
3. tilde_A - skew symmetric matrix

Answer 1

The error means that there should only be one dimension in the initial point, that is, it should be a flat numpy array. The same then also goes inside the ODE function, the interface is flat arrays, you have to establish and destroy the structure of column vectors manually. But this seems to invite strange follow-up errors, so go the other way, make everything shapeless. The rule is that matrix multiplication with a shapeless array returns a shapeless array. Do not mix, that invites "broadcasting" to a matrix.

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

With these changes the code works for me.