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

Posted 2023-03-12

【中文标题】在 scipy.integrate.solve_ivp python 中传递矩阵作为输入

【英文标题】：Passing matrices as input in scipy.integrate.solve_ivp python

【发布时间】：2021-11-25 05:45:38

【问题描述】：

我正在解决下面给出的 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 或任何此类函数？

谢谢。

【问题讨论】：

这看起来不像我多年前用于 ODE 的 MATLAB 代码。 F(t,y) 可以写成二维数组和y 数组的矩阵乘积。

这是python代码。我正在尝试在 python 中做到这一点。

您可以按顺序重写 (y1, y2, y3, y4) 向量，首先将一阶导数分组，然后是二阶导数。然后将问题一分为二：1) 通过用“ys”y[0,1] = f(....)(2D 矢量) 替换二阶导数来重写最后一组矩阵方程 2) 创建二阶导数方程 (dy[2,3]/dt = y[0,1])，也是 2D向量。然后追加两者。如果我能找到一些时间，我会尝试对其进行编码。但也许你明白了。

我的意思是，问题可能是您在 y 向量的定义中交替使用一阶和二阶导数，因此您需要手动编写所有内容。

我指的是最后一张据称是 MATLAB 的图像。

【参考方案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）。或者只是函数的全局变量。