import numpy as np
import scipy.linalg
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
##############################################################################
1. Parámetros y configuración
##############################################################################
hbar = 1.0
t_hop = 1.0
G = 1.0
M = 1.0
lambda_g = 0.2
dt = 0.02
steps = 50
factor para la evolución cuántica: exp(-iH dt / hbar)
factor_evol = -1j * dt / hbar
##############################################################################
2. Definición del icosaedro base (sin deformación)
##############################################################################
phi = (1 + np.sqrt(5)) / 2.0
base_vertices = np.array([
[-1, phi, 0],
[ 1, phi, 0],
[-1, -phi, 0],
[ 1, -phi, 0],
[ 0, -1, phi],
[ 0, 1, phi],
[ 0, -1, -phi],
[ 0, 1, -phi],
[ phi, 0, -1],
[ phi, 0, 1],
[-phi, 0, -1],
[-phi, 0, 1]
], dtype=float)
Normalizamos para que un vértice tenga radio 1
r0 = np.linalg.norm(base_vertices[0])
base_vertices /= r0
Lista de aristas del icosaedro (12 vértices, 30 aristas)
edges = [
(0,1), (0,5), (0,7), (0,10), (0,11),
(1,5), (1,7), (1,8), (1,9),
(2,3), (2,6), (2,7), (2,9), (2,11),
(3,6), (3,7), (3,8), (3,9),
(4,5), (4,6), (4,9), (4,10), (4,11),
(5,9), (5,10),
(6,7), (6,9), (6,10),
(7,8),
(8,9), (8,10),
(10,11)
]
def build_adjacency_matrix():
"""
Construye la matriz de adyacencia (12x12) con -t_hop en las posiciones
correspondientes a las aristas.
"""
A = np.zeros((12, 12), dtype=complex)
for (i, j) in edges:
A[i, j] = -t_hop
A[j, i] = -t_hop
return A
##############################################################################
3. Potencial logarítmico estático (sin deformar vértices)
##############################################################################
def build_potential_log(vertices):
"""
Calcula el potencial gravitatorio discreto V_G(i) = - sum(G*M / dist(i,j))
y luego aplica la corrección logarítmica V_log(i) = lambda_g * log(1 + V_G(i)).
Como la geometría no cambia, este potencial es constante en el tiempo.
"""
N = len(vertices)
V_G = np.zeros(N, dtype=float)
for i in range(N):
for j in range(N):
if i != j:
dist_ij = np.linalg.norm(vertices[i] - vertices[j])
V_G[i] += - (G * M / dist_ij)
# Evitamos valores no válidos en la log
min_val = np.min(V_G)
shift = 0.0
if min_val <= -1.0:
shift = -(min_val + 0.9999)
V_log = lambda_g * np.log(1.0 + V_G + shift)
return V_log
##############################################################################
4. Construcción del Hamiltoniano estático
##############################################################################
adjacency = build_adjacency_matrix()
El icosaedro permanece siempre igual (sin "update_vertices")
Por tanto, calculamos el potencial logarítmico solo una vez:
V_log_static = build_potential_log(base_vertices)
Hamiltoniano total fijo: H = adjacency + diag(V_log_static)
H_static = adjacency + np.diag(V_log_static)
Operador de evolución para un paso de tiempo dt
U_static = scipy.linalg.expm(factor_evol * H_static)
##############################################################################
5. Estado cuántico inicial y bucle de evolución
##############################################################################
psi = np.zeros(12, dtype=complex)
psi[0] = 1.0 # partícula localizada en el vértice 0 al inicio
psi_t_history = [psi.copy()]
for _ in range(steps):
psi = U_static @ psi
# Renormalizamos para evitar errores numéricos
norm_psi = np.linalg.norm(psi)
if norm_psi > 1e-12:
psi /= norm_psi
psi_t_history.append(psi.copy())
##############################################################################
6. Visualización de resultados
##############################################################################
prob_t = np.array([np.abs(state)**2 for state in psi_t_history])
Evolución de la probabilidad
plt.figure(figsize=(10,6))
for i in range(12):
plt.plot(prob_t[:, i], label=f'Vértice {i}')
plt.xlabel('Paso de tiempo')
plt.ylabel('Probabilidad')
plt.title('Evolución de la probabilidad en un icosaedro estático (potencial fijo)')
plt.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left')
plt.tight_layout()
plt.show()
Visualización 3D de la configuración final (los vértices no cambian)
final_prob = prob_t[-1]
fig = plt.figure(figsize=(7,7))
ax = fig.add_subplot(111, projection='3d')
colors = plt.cm.viridis(final_prob / final_prob.max())
for i in range(12):
x, y, z = base_vertices[i]
ax.scatter(x, y, z, color=colors[i], s=50 + 400*final_prob[i])
ax.set_title('Configuración final en el icosaedro estático')
plt.show()
Antony Padilla Morales