¿Puedo utilizar un algoritmo de búsqueda entre matrices?
He producido una gráfica de dos matrices del mismo tamaño, siendo el eje x una matriz de tiempo con diferentes pasos de tiempo y el eje y siendo una matriz de valores precalculada. Esta es la trama a continuación:
Hasta este punto, he estado buscando el momento en que delta = 0 (aparte de cuando time=0). Para hacer esto, ingresé manualmente a la matriz para encontrar en qué paso alcanza este valor (es decir, delta = 0 en el paso 3000), luego ingresé a la matriz de valor de tiempo y busqué el paso con el mismo valor (es decir, buscando el paso número 3000 y recuperando el tiempo).
¿Hay alguna manera de automatizar esto y tenerlo como resultado en lugar de buscarlo manualmente cada vez?
El código base está a continuación:
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()
Por lo general, los datos de trazado tendrán un sistema de coordenadas como entrada, es decir, una lista de tuplas (eje x, eje y). Parece que los datos de su trazado tienen pasos en el eje x. ¿Está seguro de que no hay ninguna función en el paso => ¿tiempo? ¿También puedes compartir algún código que utilices para crear este gráfico? Cualquier solución que automatice este proceso dependerá de cómo esté haciendo lo que esté haciendo
-Glubus19/03/2021 a las 15:58
"buscando el paso número 3000": está redactado como si representara mucho trabajo. Pero si es un array, entonces con arr[3000] tienes el valor deseado de una sola vez, ¿no?
-trincot19/03/2021 a las 17:33
@Glubus El código que uso tiene un paso de tiempo variable, por lo que los valores de x cambian constantemente (esto se debe al método de integración que estoy usando. El código es bastante sustancial, pero lo subiré arriba.
- Ross Alan Slater20/03/2021 a las 16:13
@trincot El problema es la repetibilidad. Cambiaré constantemente los valores de delta que a veces pueden no ser iguales a cero. El punto es obtener el paso en el que llega a cero, en todo caso.
- Ross Alan Slater20/03/2021 a las 16:13
Por lo que puedo ver en tu código, parece que d_values y t_values son siempre igualese longitud, lo que significa que si el enésimo valor en d_values es 0, eso significa que necesita el enésimo valor en t_values. ¿Es esto correcto? Si es así, entonces escribiría una función que busque 0 en d_values, manteniendo el índice (= posición en d_values) y usándolo para obtener el valor correspondiente en t_values.
-Glubus22/03/2021 a las 15:15
------------------------------------
Haz tu vida menos complicada, utiliza la rutina de bucle de tiempo solve_ivp proporcionada.
Quiere calcular un punto de distancia mínima al punto inicial u0. Este es un punto donde la derivada de la distancia va de negativa a positiva. La derivada del cuadrado dela distancia es el producto escalar del vector tangente y el vector diferencia. Cerca del mínimo no hay mucha diferencia si el vector tangente es del punto o del punto inicial.
Defina así un evento del proceso de "conteo" tipo
Tu0 = odesys(t0,u0)
def dist_plane(u): return Tu0.dot(u-u0)
event0 = lambda t,u: dist_plane(u)
event0.direction = 1
y llama al solucionador con su paso a paso estándar RK45 y todos los parámetros
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]
Las últimas líneas para que puedas hacer todos los trazados como antes.
Sin embargo, para el mínimo solo tienes que evaluar los campos de eventos de la estructura devuelta. Imprimir las medidas de distancia da como resultado una tabla
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} |")
Esta lista más corta debería ser más fácil de buscar para el mínimo global.
2
Esto funciona perfectamente con mi código actual, ¡gracias! El único problema que veo actualmente es el cambio en la malla de tiempo, ya que antes cambiaría el valor de t_step para aumentar la trayectoria general de los tres cuerpos segúnseguridad, pero hacer esto ahora parece no lograr nada, ¿cuál sería el enfoque correcto con los cambios actuales?
- Ross Alan Slater26 de marzo de 2021 a las 6:36
1
Si utiliza un solucionador con tamaño de paso adaptable, la precisión se controla mediante las tolerancias de error absolutas y relativas, el <error local>/<tamaño de paso> se controla para que sea aproximadamente max(atol, rtol*norm(y)).
- Lutz Lehmann26 de marzo de 2021 a las 6:42
------------------------------------
Si las matrices tienen la misma longitud, puede mantener el índice de su búsqueda. De lo contrario, puede utilizar el porcentaje de la matriz como punto de partida para la búsqueda de la segunda matriz. (la posición 1000 en una matriz de 3000 es 1/3, por lo que aproximadamente la posición 1166 cuando se asigna a una matriz de 3500)
Si desea algún valor (o solo habrá uno), puede dividir los datos en dos con una búsqueda binaria.
Si quieres el primero y puede haber más de un valor, tendrás que hacer una búsqueda lineal.
------------------------------------
si tuLos valores de r x son pasos de tiempo
luego, primero conviértalos a valores de tiempo resumiéndolos (pero para eso primero necesita que la matriz esté ordenada por tiempo, lo cual sería una locura si no fuera así en este caso).
si tus valores de x están ordenados
luego utilice la búsqueda binaria. si no, utilice la búsqueda lineal.
si el valor buscado no está exactamente presente en sus datos
pero tiene puntos debajo y arriba, simplemente puede interpolarlos desde los vecinos más cercanos (interpolación lineal, cúbica o 2D superior), pero para eso es mejor si la matriz buscada está ordenada, de lo contrario tendrá dificultades para obtener vecinos más cercanos válidos.