Can I use a searching algorithm between arrays?

Question

I have produced a plot of two arrays of the same size, the x axis being a time array with varying time steps and the y axis being a pre-calculated array of values. This is the plot below:

Up until this point, I have been searching for the time for when delta = 0 (apart from when time=0) . To do this, I have been manually going into the array to find at what step it reaches this value (i.e. delta = 0 at step 3,000), then going into the time-value array and looking for the step with the same value (i.e. searching for the 3,000th step and retrieving the time).

Is there a way I can automate this and have it as an output rather than manually searching each time?

The base code is below:

import numpy as np
import matplotlib.pyplot as plt

from scipy import integrate

masses = [1, 1, 1]
r1, v1 = [0, 0], [-2*0.513938054919243, -2*0.304736003875733]
r2, v2 = [-1, 0], [0.513938054919243, 0.304736003875733]
r3, v3 = [1, 0], [0.513938054919243, 0.304736003875733]
u0 = np.concatenate([r1, v1, r2, v2, r3, v3])

def odesys(t, u):
    def force(a): return a / sum(a ** 2) ** 1.5
    r1, v1, r2, v2, r3, v3 = u.reshape([-1, 2])
    m1, m2, m3 = masses
    f12, f13, f23 = force(r1 - r2), force(r1 - r3), force(r2 - r3)
    a1, a2, a3 = -m2 * f12 - m3 * f13, m1 * f12 - m3 * f23, m1 * f13 + m2 * f23
    return np.concatenate([v1, a1, v2, a2, v3, a3])

# collect data
t_values = []
u_values = []
par_1_pos = []
d_values = []


# Time start, step, and finish point
t0, tf, t_step = 0, 18, 0.0001
nsteps = int((tf - t0) / t_step)

solution = integrate.RK45(odesys, t0, u0, tf, max_step=t_step)

# The loop for running the Runge-Kutta method over some time period.
u_values.append(solution.y)
t_values.append(t0)
par_1_pos.append(((solution.y[0] - u0[0])**2 + (solution.y[1] - u0[1])**2)**0.5)

d_values.append(((solution.y[0] - u0[0])**2 + (solution.y[1] - u0[1])**2)**0.5 + 
                ((solution.y[4] - u0[4])**2 + (solution.y[5] - u0[5])**2)**0.5 + 
                ((solution.y[8] - u0[8])**2 + (solution.y[9] - u0[9])**2)**0.5 +
                ((solution.y[2] - u0[2])**2 + (solution.y[3] - u0[3])**2)**0.5 +
                ((solution.y[6] - u0[6])**2 + (solution.y[7] - u0[7])**2)**0.5 +
                ((solution.y[10] - u0[10])**2 + (solution.y[11] - u0[11])**2)**0.5)

for step in range(nsteps):
    solution.step()
    u_values.append(solution.y)
    t_values.append(solution.t)
    
    par_1_pos.append(((solution.y[0] - u0[0])**2 + (solution.y[1] - u0[1])**2)**0.5)
    
    d_values.append(((solution.y[0] - u0[0])**2 + (solution.y[1] - u0[1])**2)**0.5 + 
                ((solution.y[4] - u0[4])**2 + (solution.y[5] - u0[5])**2)**0.5 + 
                ((solution.y[8] - u0[8])**2 + (solution.y[9] - u0[9])**2)**0.5 +
                ((solution.y[2] - u0[2])**2 + (solution.y[3] - u0[3])**2)**0.5 +
                ((solution.y[6] - u0[6])**2 + (solution.y[7] - u0[7])**2)**0.5 +
                ((solution.y[10] - u0[10])**2 + (solution.y[11] - u0[11])**2)**0.5)
    
    # break loop after modelling is finished
    if solution.status == 'finished':
        break

# Plotting of the individual particles
u = np.asarray(u_values).T

# Plot for The trajectory of the three bodies over the time period
plt.plot(u[0], u[1], '-o', lw=1, ms=3, label="body 1")
plt.plot(u[4], u[5], '-x', lw=1, ms=3, label="body 2")
plt.plot(u[8], u[9], '-s', lw=1, ms=3, label="body 3")

plt.title('Trajectories of the three bodies')
plt.xlabel('X Position')
plt.ylabel('Y Position')

plt.legend()
plt.grid()
plt.show()
plt.close()

# Plot for d(delta_t) values
plt.plot(t_values, d_values)

plt.title('Delta number for the three bodies')
plt.xlabel('Time (s)')
plt.ylabel('Delta')

plt.grid()
plt.show()
plt.close()

# Plot of distance between P1 and IC
plt.plot(t_values, par_1_pos)

plt.title('Plot of distance between P1 and IC')
plt.xlabel('Time (s)')
plt.ylabel('Distance from origin')

plt.grid()
plt.show()
plt.close()

Answer 1

Make your life less complicated, use the provided time-loop routine solve_ivp.

You want to compute a point of minimal distance to the initial point u0. This is a point where the derivative of the distance goes from negative to positive. The derivative of the square of the distance is the dot product of tangent vector and difference vector. Close to the minimum it makes not much difference if the tangent vector is of the point or of the initial point.

Thus define an event of the "counting" type

Tu0 = odesys(t0,u0)
def dist_plane(u): return Tu0.dot(u-u0)
event0 = lambda t,u: dist_plane(u)
event0.direction = 1

and call the solver with its standard stepper RK45 and all the parameters

solution = solve_ivp(odesys, [t0,tf], u0, events=event0, atol=1e-10, rtol=1e-11)

u_values = solution.y.T
t_values = solution.t

def norm(u): return sum(u**2)**0.5

def dist_1(u): return norm(u[:2]-u0[:2])
def dist_all(u): return sum(norm(uu) for uu in (u-u0).reshape([-1,2]))

par_1_pos = [ dist_1(uu) for uu in u_values]
d_values = [dist_all(uu) for uu in u_values]
d_plane = [dist_plane(uu) for uu in u_values]

The last lines so you can do all the plots as before.

For the minimal however you just have to evaluate the event fields of the returned structure. Printing the distance measures results in a table

for tt,uu in zip(solution.t_events[0], solution.y_events[0]):
    print(rf"|{tt:12.8f} | {dist_1(uu):12.8f} | {dist_all(uu):12.8f} | {dist_plane(uu):12.8f} |")

t	dist pos1	dist all	dist deriv
0.00000000	0.00000000	0.00000000	0.00000000
2.81292221	0.58161236	3.54275380	0.00000000
4.17037860	0.35583855	5.77531098	0.00000000
5.71976151	0.63111430	3.98764796	-0.00000000
8.66440460	0.00000019	3.73331800	0.00000000
11.60904921	0.63111445	3.98764485	-0.00000000
13.15843018	0.35583804	5.77530951	0.00000000
14.51588605	0.58161284	3.54275265	-0.00000000
17.32881023	0.00000078	0.00000328	0.00000000

This shorter list should be easier to search for the global minimum.

Answer 2

1. if your x values are time steps

   then convert them to time values first by summing them up (but for that you need the array is sorted by time first which would be insane if not already in this case).

2. if your x values are sorted

   then use binary search. if not use linear search.

3. if your searched value is not exactly present in your data

   but you have points below and above you can simply interpolate it from closest neighbors (linear,cubic or higher 2D interpolation) but for that is best if your searched array is sorted otherwise you will have hard time getting valid closest neighbors.

Answer 3

If the arrays are the same length, you can keep the index of your search.. Otherwise, you may be able to use the percentage through of the array as your starting point for the search of the second array. (position 1000 in an array of 3000 is 1/3rd, so about position 1166 when mapped into an array of 3500)

If you want any value (or there will only be one), you can bisect the data with a binary search.

If you want the first and there may be more than one value, you'll have to do a linear search.