我试图用Python中的solve_ivp函数来计算姿态运动方程的常微分方程，但问题是其中一个参数，角速度ω，发生了变化，我想把这个考虑进去，我之前已经从另一个常微分方程中计算了ω，现在我想把得到的结果作为另一个常微分方程的输入。
我是这么做的：

def fun2(time, euler):
    omegax = 5.2928e-10; omegay = -2.5347e-11; omegaz = 2.6609e-6
    dot1 = (omegax*np.sin(euler[2]) + omegay*np.cos(euler[2]))/np.sin(euler[1])
    dot2 = omegax*np.cos(euler[2]) - omegay*np.sin(euler[2])
    dot3 = omegaz - (omegax*np.sin(euler[2]) + omegay*np.cos(euler[2]))/np.tan(euler[1])

    return np.array([dot1, dot2, dot3])

angles = integrate.solve_ivp(fun2, tspan, euler0, t_eval = t, method = 'RK45', dense_output = True, rtol=1e-13, atol=1e-22)

这里我使用了一个常量omega来运行代码，但是我希望它可以改变。我在另一个ODE中得到的omega（总是使用solve_ivp）是一个矩阵的形式，其中所有的omegax，y和z分别在第1，2和3列，有1000000行。我尝试过在fun2中求解前面的ODE，如下所示：

def fun2(time, euler):
    x_t = integrate.solve_ivp(fun, tspan, omega0, t_eval=t, method='RK45', dense_output=True, rtol=1e-13, atol=1e-22)
    omegax = x_t.y[0];   omegay = x_t.y[1];   omegaz = x_t.y[2]
    dot1 = (omegax*np.sin(euler[2]) + omegay*np.cos(euler[2]))/np.sin(euler[1])
    dot2 = omegax*np.cos(euler[2]) - omegay*np.sin(euler[2])
    dot3 = omegaz - (omegax*np.sin(euler[2]) + omegay*np.cos(euler[2]))/np.tan(euler[1])

    return np.array([dot1, dot2, dot3])

angles = integrate.solve_ivp(fun2, tspan, euler0, t_eval = t, method = 'RK45', dense_output = True, rtol=1e-13, atol=1e-22)

不幸的是，它没有工作，我得到这个错误消息："操作数无法与形状一起广播（3，1000000）（3，）"，现在我卡住了。有人能帮我吗？