将矩阵作为scipy.integrate.solve_ivp python中的输入传递
将矩阵作为scipy.integrate.solve_ivp python中的输入传递
提问于 2021-10-04 20:31:37
我正在求解一个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或任何类似的函数来做同样的事情？
谢谢。
回答 1
Stack Overflow用户
回答已采纳
发布于 2021-10-04 21:04:00
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)。或者仅仅是函数的全局变量。
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/69442104
复制相关文章
相似问题
腾讯云开发者
