理解scipy.integrate.RK45要求的问题

rn0zuynd  于 2022-11-10  发布在  其他
关注(0)|答案(1)|浏览(271)

我试图用scipy.integrate.RK45()来求解一阶微分方程组。我已经编写出了我希望绘制的模型函数(位移与时间),但RK45()要求该函数有两个参数,即'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

def model(t, y):

    m = 2
    c = 10
    k = 1500

    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

time_step = 0.01
t_upper = 0.5
t_lower = 0

initial_conditions = [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)是什么?我如何将它们有效地合并到我的函数中?

niknxzdl

niknxzdl1#

你已经使用了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中的矩阵乘积。*执行元素式乘积。
编辑:我还没有发现RK45(fun=model(t, y),调用。这样做,你传递的是t,y处的函数模型的值。你需要给予函数本身:RK45(fun=model, ...

相关问题