CODEWITHSUNDEEP
pythonscipysimulationode
relent95
relent95

Reputation: 4731

How can I use scipy.integrate.solve_ivp() to solve an ODE in an interactive simulation?

I have implemented a simple simulation on a uniform circular motion using the low level API scipy.integrate.RK45() like the following.

import numpy as np
import scipy.integrate
import matplotlib.pyplot as plt

r = np.array([1, 0], 'float')
v = np.array([0, 1], 'float')
dt = 0.1

def motion_eq(t, y):
    r, v = y[0:2], y[2:4]
    return np.hstack([v, -r])
motion_solver = scipy.integrate.RK45(motion_eq, 0, np.hstack([r, v]),
    t_bound = np.inf, first_step = dt, max_step = dt)

particle, *_ = plt.plot(*r.T, 'o')
plt.gca().set_aspect(1)
plt.xlim([-2, 2])
plt.ylim([-2, 2])
def update():
    motion_solver.step()
    r = motion_solver.y[0:2]
    particle.set_data(*r.T)
    plt.draw()
timer = plt.gcf().canvas.new_timer(interval = 50)
timer.add_callback(update)
timer.start()

plt.show()

I tried at first the high level API scipy.integrate.solve_ivp(), but it seems it doesn't provide an interface to create an instance containing the state of the system and get states of the system iteratively.(I am calling this an interactive simulation, because you can pause, change the system state and resume, although it's not implemented in the sample code.)

Is this possible with the solve_ivp()? If it is not, am I doing right with the RK45, especially on specifying the t_bound, first_step and max_step options? I can find plenty of resources on solving for a given time interval on Internet, but I couldn't find one on solving like this.

Upvotes: 1

Views: 2094

Answers (1)

relent95
relent95

Reputation: 4731

Although I'm not an expert on the numerical analysis, I'm answering my own question, because I have come to a conclusion after analyzing the solve_ivp() and RK45 implementation and experimenting with APIs.

Firstly, the high level API solve_ivp() does not provide an interface to create an instance containing the state of the system and get states of the system iteratively. So, I should do like this even though it's inefficient.

dt = 0.1
t = 0
r = np.array([1, 0], 'float')
v = np.array([0, 1], 'float')

def motion_eq(t, y):
    r, v = y[0:2], y[2:4]
    return np.hstack([v, -r])

...
def update():
    global t, r, v
    sol = scipy.integrate.solve_ivp(motion_eq, (t, t + dt), np.hstack([r, v]),
        t_eval = (t + dt,), method = 'RK45')
    t = sol.t[0]
    r, v = sol.y[0:2, 0], sol.y[2:4, 0]
    ...
...

Secondly, the low level API RK45 provides the interface to create an instance containing the state of the system(the RK45() initializer) and get states of the system iteratively(the step() method), but it does not provide an interface to control the time step at each iteration. But if the desired time step is constant and small enough, you can control the time step in some degree via the first_step and max_step arguments of the RK45() initializer, like the example in the question.

This can be seen in this and this. At the end of the _step_impl(), the self.h_abs gets larger then the previous value before the call. But at the start of the _step_impl() on the next call, it is limited to the self.max_step. But as Lutz Lehmann said, if the estimated error(the error_norm) is large, the self.h_abs at the next call will be smaller than the self.max_step. So in this case, it is required to loop and interpolate like the following.(The interpolation part is attributed to Lutz Lehmann.)

...
dt = 0.1
t = 0
...
motion_solver = scipy.integrate.RK45(motion_eq, 0, np.hstack([r, v]),
    t_bound = np.inf, first_step = dt, max_step = dt)
...
def update():
    global t
    t += dt
    while motion_solver.t < t:
        motion_solver.step()
    sol = motion_solver.dense_output()
    y = sol(t) # interpolation
    r, v = y[0:2], y[2:4]
    ...
...

Upvotes: 1

Related Questions