问题理解scipy.integrate.RK45要求
Issue understanding scipy.integrate.RK45 requirements
我正在尝试用 scipy.integrate.RK45()
求解一阶微分方程组。我已经编写了我希望绘制的模型函数(位移与时间),但是 RK45()
要求此函数采用 2 个参数,即 't' 和 'y',其中 't' 是一个标量, 'y' 在我的例子中是一个数组。这在下面有更好的描述:
https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.RK45.html
我的脚本如下:
import numpy as np
from scipy.integrate import RK45, RK23
# Model the calling function [X_dot]:
#------------------------------------------------------------------
def model(t, y):
# Define the mass, stiffness and damping [m, c, k]:
m = 2
c = 10
k = 1500
# Define the piecewise forcing function [F]:
if (t >= 0 and t < 0.1):
F = 200 * t
if (t >= 0.1 and t < 0.25):
F = 20
else:
F = 0
E_matrix = [[0, 1], [(-k / m), (-c / m)]]
Q_matrix = [0, F / m]
return E_matrix*X + Q_matrix
# Define the initial conditions and integration boundaries:
#------------------------------------------------------------------
time_step = 0.01
t_upper = 0.5
t_lower = 0
initial_conditions = [0, 0] # [displacement(t=0) = 0, velocity(t=0) = 0]
points_to_plot = RK45(fun=model(t, y), t0=t_lower, y0=initial_conditions, t_bound=t_upper, vectorized=True)
下面是我要解决的系统图片:
我发现很少有这样的例子,因为大多数解决方案都使用 odeint()
。
这两个参数 (t, y) 是什么?我如何将它们有效地合并到我的函数中?非常感谢任何帮助。
您已经使用过 t
。现在,将 def model(t, y):
更改为 def model(t, X):
,您也将使用 X
。请注意,t 和 y 是位置参数,您可以在函数中随意调用它们。您还有另一个问题,那就是您乘以 Python 列表!在Python中,相对于Matlab,你需要指定你做一个数组:
改变
E_matrix = [[0, 1], [(-k / m), (-c / m)]]
Q_matrix = [0, F / m]
到
E_matrix = np.array([[0, 1], [(-k / m), (-c / m)]])
Q_matrix = np.array([0, F / m])
和
return E_matrix*X + Q_matrix
到
return E_matrix @ X + Q_matrix
因为@
是NumPy中的矩阵乘积。 *
执行 element-wise 产品。
编辑:我没有发现调用 RK45(fun=model(t, y),
。通过这样做,您将在 t,y
处传递函数模型的值。你需要给函数本身:RK45(fun=model, ...
我正在尝试用 scipy.integrate.RK45()
求解一阶微分方程组。我已经编写了我希望绘制的模型函数(位移与时间),但是 RK45()
要求此函数采用 2 个参数,即 't' 和 'y',其中 't' 是一个标量, 'y' 在我的例子中是一个数组。这在下面有更好的描述:
https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.RK45.html
我的脚本如下:
import numpy as np
from scipy.integrate import RK45, RK23
# Model the calling function [X_dot]:
#------------------------------------------------------------------
def model(t, y):
# Define the mass, stiffness and damping [m, c, k]:
m = 2
c = 10
k = 1500
# Define the piecewise forcing function [F]:
if (t >= 0 and t < 0.1):
F = 200 * t
if (t >= 0.1 and t < 0.25):
F = 20
else:
F = 0
E_matrix = [[0, 1], [(-k / m), (-c / m)]]
Q_matrix = [0, F / m]
return E_matrix*X + Q_matrix
# Define the initial conditions and integration boundaries:
#------------------------------------------------------------------
time_step = 0.01
t_upper = 0.5
t_lower = 0
initial_conditions = [0, 0] # [displacement(t=0) = 0, velocity(t=0) = 0]
points_to_plot = RK45(fun=model(t, y), t0=t_lower, y0=initial_conditions, t_bound=t_upper, vectorized=True)
下面是我要解决的系统图片:
我发现很少有这样的例子,因为大多数解决方案都使用 odeint()
。
这两个参数 (t, y) 是什么?我如何将它们有效地合并到我的函数中?非常感谢任何帮助。
您已经使用过 t
。现在,将 def model(t, y):
更改为 def model(t, X):
,您也将使用 X
。请注意,t 和 y 是位置参数,您可以在函数中随意调用它们。您还有另一个问题,那就是您乘以 Python 列表!在Python中,相对于Matlab,你需要指定你做一个数组:
改变
E_matrix = [[0, 1], [(-k / m), (-c / m)]]
Q_matrix = [0, F / m]
到
E_matrix = np.array([[0, 1], [(-k / m), (-c / m)]])
Q_matrix = np.array([0, F / m])
和
return E_matrix*X + Q_matrix
到
return E_matrix @ X + Q_matrix
因为@
是NumPy中的矩阵乘积。 *
执行 element-wise 产品。
编辑:我没有发现调用 RK45(fun=model(t, y),
。通过这样做,您将在 t,y
处传递函数模型的值。你需要给函数本身:RK45(fun=model, ...