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

Question

我正在求解下面给出的 2DOF spring-质量阻尼器系统：

这些是 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 或任何此类函数做同样的事情吗？

谢谢。

Answer 1

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,) 传递（它独立于 t 和 y）。或者只是函数的全局变量。

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

Passing matrices as input in scipy.integrate.solve_ivp python

python

matlab

scipy

differential-equations

ode45