时间戳

首页 > 解决方案 > 在 scipy.integrate.solve_ivp python 中传递矩阵作为输入

python - 在 scipy.integrate.solve_ivp python 中传递矩阵作为输入

问题描述

我正在解决下面给出的 2DOF 弹簧质量阻尼器系统： 在此处输入图像描述

这些是 2 个控制方程 在此处输入图像描述

我已经通过以下方式解决了它：

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
m1 = 3
m2 = 5

k1 = 7
k2 = 9

c1 = 1
c2 = 2

f1 = 40
f2 = 4

start_time = 0
end_time = 60

initial_position_m_1 = 6
initial_velocity_m_1 = 0

initial_position_m_2 = 9
initial_velocity_m_2 = 4

delta_t = 0.1
def F(t, y):
    arr = np.array([
        y[1],
        (1/m1)*(f1*np.cos(3*t) - ((c1 + c1)*y[1] + (k1 + k2)*y[0]) + c2*y[3] + k2*y[2]),
        y[3],
        (1/m2)*(f2*np.sin(t**2) - c2*y[3] - k2*y[2] + c2*y[1] + k2*y[0])
    ])
    return arr

time_interval = np.array([start_time, end_time])
initial_conditions = np.array([initial_position_m_1, initial_velocity_m_1, initial_position_m_2, initial_velocity_m_2])

#######     solving the system of equations    ####
sol = solve_ivp(F, time_interval, initial_conditions, max_step = delta_t)

T = sol.t
Y = sol.y

现在，这是通过将 2 个控制方程转换为 4 个方程来完成的，如下所示： 在此处输入图像描述

问题是我必须分别编写每个方程（作为函数 F）

Matlab 有一种方法可以仅使用 Ode45 函数通过矩阵求解，即您不必在 Matlab 的函数 F 中单独编写所有方程。您可以在其中输入质量、刚度和阻尼系数作为矩阵。像这样： 在此处输入图像描述

我正在尝试解决一个涉及 30x30 矩阵的问题，如果我以上述方式解决问题，我将不得不为函数 F 编写 60 个单独的方程，而在 Matlab 中，我可以将先前计算的 30x30 矩阵直接传递给函数。有什么方法可以对 python 中的solve_ivp 或任何此类函数做同样的事情吗？

谢谢你。

标签： pythonmatlabscipydifferential-equationsode45

解决方案

arr = np.array([
        y[1],
        (1/m1)*(f1*np.cos(3*t) - ((c1 + c1)*y[1] + (k1 + k2)*y[0]) + c2*y[3] + k2*y[2]),
        y[3],
        (1/m2)*(f2*np.sin(t**2) - c2*y[3] - k2*y[2] + c2*y[1] + k2*y[0])
    ])

可以写成（大致）：

f = np.array([0, f1*np.cos(3*t),0,f2*np.sin(t**2)])
M = np.array([
        [0, 1, 0, 0],
        [(k1+K2), (c1+c1), k2, c2],
        [0,0,0,1],
        [k2, c2, ....]])
arr = f[:,None] + M.dot(y)

该M数组可以通过args=(M,)（它独立于tand y）传递。或者只是作为函数的全局变量。

推荐阅读