I want to integrate a particle path in 2D using the integrate.ode module
. Things that are a bit different in my case are that, I only want to integrate up to a certain position, determined by the maximum allowed x coordinate of the particle: x_max
.
The main issue I have is that the particle may first move very slowly and then gather more speed later on. Hence I don't want to waste effort with small time steps in this region. The algorithm should be able to adjust such that smaller time steps are used when the particle velocity becomes high.
In the end result if I plot the particle trajectory in the "phase-space" I should have a smooth line.
I have some rough pseudo-code below for this purpose:
backend = "dopri5"
x_max = 1
solver = ode(f)
solver.set_integrator(backend)
solver.set_initial_value(y0, t0)
t, y = [t0], [y0]
k = 1.2
while solver.successful() and solver.y[0] < x_max
solver.integrate(solver.t+dt)
t.append(solver.t)
y.append(solver.y)
v_current = numpy.linalg.norm(y[-1])
v_previous = numpy.linalg.norm(y[-2])
if numpy.abs( v_current-v_previous ) > k * v_previous:
dt = 0.8*dt
del y[-1]
else:
dt = dt*1.2
Trouble is this algorithm may not be that robust, as choosing the values k, 1.2, 0.8
is somewhat arbitrary and may cause some stability issues with the algorithm.
EDIT: I also want to be able to plot points that are equally spaced in time on the trajectory, to give an indication of how the speed of the particle changes.
Can anyone suggest a better way of doing this?