I am trying to numerically integrate the following equation where the omega(3x1) vector and theta are (4X1) vector and W(w) is 4x4 skew-symmetric matrix
def skew_omega(w):
"""skew symmetric omega (4x4)"""
print(w)
return np.array([[0, -w[0], -w[1], -w[2]],
[w[0], 0, w[2], -w[1]],
[w[1], -w[2], 0, w[0]],
[w[2], w[1], -w[0], 0]])
def theta_model(euler_param, w, t):
"""Ode for theta"""
print(w)
theta = 0.5*np.matmul(skew_omega(w),euler_param)
return theta
#initial conditions
t = np.array([0, 1])
theta0 = np.array([1, 0, 0, 0])
omega0 = np.array([0.1, 10, 0])
#solve ode
theta = odeint(theta_model, theta0, t, args = (omega0,))
0.0
0.0
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-72-1324ec16154f> in <module>
----> 1 theta = odeint(theta_model, theta0, t, args = (om,))
W:\Anaconda\lib\site-packages\scipy\integrate\odepack.py in odeint(func, y0, t, args, Dfun, col_deriv, full_output, ml, mu, rtol, atol, tcrit, h0, hmax, hmin, ixpr, mxstep, mxhnil, mxordn, mxords, printmessg, tfirst)
243 full_output, rtol, atol, tcrit, h0, hmax, hmin,
244 ixpr, mxstep, mxhnil, mxordn, mxords,
--> 245 int(bool(tfirst)))
246 if output[-1] < 0:
247 warning_msg = _msgs[output[-1]] + " Run with full_output = 1 to get quantitative information."
<ipython-input-55-c581c2c300b8> in theta_model(euler_param, w, t)
1 def theta_model(euler_param, w, t):
2 print(w)
----> 3 theta = 0.5*np.matmul(skew_omega(w),euler_param)
4
5 return theta
<ipython-input-49-bb06a7f0e258> in skew_omega(w)
1 def skew_omega(w):
2 print(w)
----> 3 return np.array([[0, -w[0], -w[1], -w[2]],
4 [w[0], 0, w[2], -w[1]],
5 [w[1], -w[2], 0, w[0]],
TypeError: 'float' object is not subscriptable
The printed zeros((the print statements in the functions) indicate that the odeint solver is not taking in the inital value- omega0. Instead it gives a 0.0 float as a place holder and hence the TypeError is thrown.
Why the ode solver not taking omega0 as input?
From odeint
documentation
Parameters
----------
func : callable(y, t, ...) or callable(t, y, ...)
Computes the derivative of y at t.
Your model function
def theta_model(euler_param, w, t):
Looks like you have to exchange t
and w
in the argument list.