如何正确更新我的集成函数数组？

Question

我一直在尝试使用 Runge-Kutta45 积分方法更新 space 中的一组粒子位置和速度，以获得某个时间步长的新状态。

最初，我创建了一个包含所有这些元素的数组并将它们组合 (y)：

r_vec_1 = np.array([0, 0])
v_vec_1 = np.array([-np.sqrt(2), -np.sqrt(2)])

r_vec_2 = np.array([-1, 0])
v_vec_2 = np.array([np.sqrt(2) / 2, np.sqrt(2) / 2])

r_vec_3 = np.array([1, 0])
v_vec_3 = np.array([np.sqrt(2) / 2, np.sqrt(2) / 2])

y_0 = np.concatenate((r_vec_1, v_vec_1, r_vec_2, v_vec_2, r_vec_3, v_vec_3))
y = y_0

现在，我使用这个数组作为我的初始条件并创建了一个函数，该函数为我提供了一个名为 F(y) 的新函数，它是我的函数 y 的导数，用一组一阶 ODE 表示：

def fun(t,y):
    np.array([y[2], y[3], x1_double_dot(y, mass_vector), y1_double_dot(y, mass_vector),
              y[6], y[7], x2_double_dot(y, mass_vector), y2_double_dot(y, mass_vector),
              y[10], y[11], x3_double_dot(y, mass_vector), y3_double_dot(y, mass_vector)])

获得新函数文件后，我使用了初始时间和最终时间以及 scipy.integrate.RK45 子例程所需的时间步长，得到以下代码：

#This will keep track of the y vector over time
history = [y_0]

#Time start, step, and finish point
t_0 = 0
t = 0
t_step = 0.01
t_final = 200
nsteps = int((t_final - t_0)/t_step)

#The loop for running the Runge-Kutta method over some time period.
for step in range(nsteps):
    y_new = sp.integrate.RK45(fun, t_0, y_0, t_final)
    history.append([y_new.y])
    y = y_new.y
    t += t_step

我要对RK45函数进行积分，然后将得到的新的y函数放回RK45函数中重新连续积分，直到时间结束

目前，该方法为我的历史数组中的所有值输出初始 y 函数 (y_0)，但我想要一个包含该时间段内所有 y 函数的更新列表。

任何建议将不胜感激，祝你有美好的一天！

Answer 1

class RK45 是一个步进器，它不会从自身集成。在循环中，您只是为这个 class 的新实例重复调用构造函数。由于您没有调用 .step() 方法，因此步进器保留在初始点的数据中。

stepper = sp.integrate.RK45(fun, t_0, y_0, t_final)
for step in range(nsteps):
    stepper.step()
    history.append([stepper.t,*stepper.y])

您可能会对步进器 classes 实现自适应步长控制感兴趣，因此您还需要记录时间。

或者，如果您想使用 t_eval 选项模拟 solve_ivp 的基本功能，请这样做

t0,tf,dt = 0,10,0.1
stepper = RK45(lambda t,y: [y[1],-np.sin(y[0])],t0,[0.5,0],tf)
history = [stepper.y]
times = np.arange(t0,tf,dt)
for t in times[1:]:
    if t > stepper.t:
        while t>stepper.t: stepper.step()
        sol = stepper.dense_output()
    history.append(sol(t))
plt.plot(times, history); plt.grid()

如何正确更新我的集成函数数组？

How do I update my integrated function array correctly?

python

arrays

numpy

scipy

runge-kutta