How do I update my integrated function array correctly?

Question

I have been trying to use the Runge-Kutta45 integration method to update a set of positions and velocities of particles in space to get the new state at some time step.

Initially, I created an array with all these elements and combined them (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

Now, I used this array as my initial conditions and created a function that gave me a new function called F(y) which is the derivative of my function y represented in a set of 1st order ODEs:

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)])

Once I had obtained the new function file, I used an initial and final time as well as a times step which is needed in the scipy.integrate.RK45 subroutine, resulting in the following code:

#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

I want the RK45 function to integrate, then put the new y function obtained back into the RK45 function to be integrated again continuously until the time period has expired.

Currently, the method outputs the initial y function (y_0) for all values in my history array but I want an updated list of all y functions over the time period.

Any advice would be greatly appreciated, have an awesome day!

Answer 1

The class RK45 is a stepper, it does not integrate from itself. In the loop you were just calling repeatedly the constructor for a new instance of this class. As you did not call the .step() method, the stepper remains in its data at the initial point.

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

You might be interested that the stepper classes implement an adaptive step size control, so you also need to record the times.

Or if you want to simulate the base functionality of solve_ivp with the t_eval option, do it like this

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()