对称性决定的簇同步（星图中的“远程同步”）

# 星图：中心节点与 6 个叶节点都是 Stuart–Landau 振子 dz/dt = (1 + i w - |z|^2) z + K sum_j A_ij (z_j - z)。
# 中心固有频率 w_h=3，叶节点 w=1，K=0.3。叶节点之间没有直接连接，但它们在图的自同构下
# 属于同一轨道（可以任意互换），因此“所有叶节点状态相同”是动力学的不变子空间。
# 先用 networkx 列出自同构轨道，再模拟，比较叶节点之间、叶节点与中心之间的状态差。
# 解读：叶节点彼此完全同步（状态差≈0），中心与它们频率锁定，但振幅和相位都不同，
#       不属于这个簇——互不相连的叶节点因对称位置而同步，这就是簇同步/远程同步。
import numpy as np
import networkx as nx
from networkx.algorithms.isomorphism import GraphMatcher

G = nx.star_graph(6)                                  # 节点 0 为中心
orbits = {}
for iso in GraphMatcher(G, G).isomorphisms_iter():
    for u, v in iso.items():
        orbits.setdefault(u, set()).add(v)
print("自同构轨道:", sorted({tuple(sorted(o)) for o in orbits.values()}))

A = nx.to_numpy_array(G, nodelist=range(7)); deg = A.sum(1)
w = np.r_[3.0, np.ones(6)]
rng = np.random.default_rng(2)
z = rng.standard_normal(7) + 1j * rng.standard_normal(7)
K, dt, steps = 0.3, 0.005, 60000
phases = []
for n in range(steps):
    z = z + dt * ((1 + 1j * w - np.abs(z) ** 2) * z + K * (A @ z - deg * z))
    if n >= steps // 2:
        phases.append(np.angle(z))
ph = np.unwrap(np.array(phases), axis=0)
freq = (ph[-1] - ph[0]) / ((len(ph) - 1) * dt)
print("叶节点之间最大状态差 max|z_i - z_j| = %.1e" % np.abs(z[1:, None] - z[None, 1:]).max())
print("振幅：中心 %.2f，叶节点 %.2f；中心与叶节点相位差 %.2f rad" % (abs(z[0]), abs(z[1]), np.angle(z[0] / z[1])))
print("平均角频率：中心 %.3f，叶节点" % freq[0], np.round(freq[1:], 3).tolist())

import matplotlib.pyplot as plt
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4), gridspec_kw={"width_ratios": [1, 1.6]})
pos = nx.spring_layout(G, seed=0)
nx.draw_networkx_edges(G, pos, ax=ax1, edge_color="0.6", alpha=0.4)
nx.draw_networkx_nodes(G, pos, ax=ax1, node_size=120, node_color=["tab:red"] + ["tab:blue"] * 6)
ax1.set_title("Orbits: hub (red) vs leaves (blue)"); ax1.axis("off")
tt = np.arange(len(ph))[-1500:] * dt
for i in range(7):
    ax2.plot(tt, np.cos(ph[-1500:, i]), c="tab:red" if i == 0 else "tab:blue",
             lw=2.5 if i == 0 else 1.2, ls="-" if i == 0 else "--", label=("hub" if i == 0 else ("leaves" if i == 1 else None)))
ax2.set_xlabel("time (last window)"); ax2.set_ylabel("cos(phase)"); ax2.legend(loc="upper right")
ax2.set_title("Unconnected leaves overlap exactly: cluster synchrony")
plt.tight_layout(); plt.show()