数学推导与算法

本页补全以下推导：A. Kuramoto 自洽方程与 $K_c=2/(\pi g(0))$；B. 共识方程的显式解与 $\mu_2$ 的作用；C. 主稳定函数如何导出特征值比判据；D. 洛伦兹分布下的 Ott–Antonsen 约化；E. 簇同步的商网络；以及三个数值算法。

A. Kuramoto 自洽方程与临界耦合

假设：全连接、$N\to\infty$；$g(\omega)$ 关于 0 对称、单峰（在 $\omega=0$ 处取最大值）；只寻找 $r$ 与 $\psi$ 不随时间变化的稳态解，不讨论其稳定性。由对称性可取平均频率为 0，并选坐标使 $\psi=0$。

第 1 步：两类振荡器。 平均场方程变为 $\dot\theta_i=\omega_i-Kr\sin\theta_i$。

• 若 $|\omega_i|\le Kr$，存在不动点 $\sin\theta_i=\omega_i/(Kr)$；其中满足 $\cos\theta_i\ge0$ 的那个稳定，即 $\theta_i=\arcsin(\omega_i/Kr)\in[-\pi/2,\pi/2]$。这些是锁定振荡器。
• 若 $|\omega_i|>Kr$，$\dot\theta_i$ 不变号，振荡器绕圆漂移，在 $\theta$ 处停留的时间与 $1/|\dot\theta|$ 成正比，平稳密度 $\rho(\theta\,|\,\omega)\propto1/|\omega-Kr\sin\theta|$。由于 $\rho(\theta+\pi\,|\,-\omega)=\rho(\theta\,|\,\omega)$ 且 $g$ 对称，漂移振荡器对 $\langle\cos\theta\rangle$ 和 $\langle\sin\theta\rangle$ 的贡献两两抵消。

第 2 步：自洽条件。 序参量 $r=\langle e^{i\theta}\rangle$ 只剩锁定振荡器的贡献；虚部因 $g$ 对称而为零，实部为 $$r=\int_{-Kr}^{Kr}\cos\theta(\omega)\,g(\omega)\,d\omega .$$ 换元 $\omega=Kr\sin\theta$，$d\omega=Kr\cos\theta\,d\theta$： $$r=Kr\int_{-\pi/2}^{\pi/2}\cos^2\theta\;g(Kr\sin\theta)\,d\theta .$$ 这就是自洽方程。$r=0$ 永远是解（不相干态）。

第 3 步：分岔点。 若 $r>0$，两边除以 $r$： $$1=K\int_{-\pi/2}^{\pi/2}\cos^2\theta\;g(Kr\sin\theta)\,d\theta .$$ 令 $r\to0^+$，$g(Kr\sin\theta)\to g(0)$，且 $\int_{-\pi/2}^{\pi/2}\cos^2\theta\,d\theta=\pi/2$，得 $1=K g(0)\pi/2$，即 $$K_c=\frac{2}{\pi g(0)} .$$

第 4 步：临界点附近的 $r$。 把 $g$ 在 0 处展开到二阶，$g(Kr\sin\theta)\approx g(0)+\tfrac12 g''(0)K^2r^2\sin^2\theta$，并用 $\int_{-\pi/2}^{\pi/2}\cos^2\theta\sin^2\theta\,d\theta=\pi/8$： $$1=\frac{K}{K_c}+\frac{\pi}{16}g''(0)K^3r^2 \;\Longrightarrow\; r^2\approx\frac{16\,(K-K_c)}{\pi K_c^4\,[-g''(0)]}.$$ 单峰分布 $g''(0)<0$，所以只有 $K>K_c$ 时才有实解，且 $r\propto\sqrt{K-K_c}$——这就是Kuramoto 模型与序参量中所说的连续相变及其平方根标度。对洛伦兹分布，自洽积分可以精确计算，得到 $r=\sqrt{1-2\gamma/K}$（另一种推导见 D 节）。

B. 共识方程的显式解与代数连通度

设网络无向、连通，边权 $w_{ij}=w_{ji}\ge0$，$L=D-A$。

第 1 步：$L$ 半正定。 对任意 $\mathbf x$， $$\mathbf x^\top L\mathbf x=\sum_{(i,j)\in E}w_{ij}(x_i-x_j)^2\ge0,$$ 等号成立当且仅当每条边两端相等；连通时即 $\mathbf x\propto\mathbf 1$。所以 0 是单重特征值，$\mu_2>0$。

第 2 步：谱分解求解。 $L$ 对称，有正交归一特征向量 $\mathbf v_1=\mathbf 1/\sqrt N,\mathbf v_2,\dots,\mathbf v_N$。线性方程 $\dot{\mathbf x}=-L\mathbf x$ 的解为 $$\mathbf x(t)=e^{-Lt}\mathbf x(0)=\sum_{k=1}^{N}e^{-\mu_kt}\,(\mathbf v_k^\top\mathbf x(0))\,\mathbf v_k =\bar x\,\mathbf 1+\sum_{k=2}^{N}e^{-\mu_kt}c_k\mathbf v_k ,$$ 其中 $c_k=\mathbf v_k^\top\mathbf x(0)$，$\bar x=\frac1N\sum_ix_i(0)$。平均值守恒，因为 $\frac{d}{dt}\mathbf 1^\top\mathbf x=-\mathbf 1^\top L\mathbf x=0$。

第 3 步：收敛速度。 偏差 $\boldsymbol\delta(t)=\mathbf x(t)-\bar x\mathbf 1$ 满足 $$\|\boldsymbol\delta(t)\|^2=\sum_{k\ge2}e^{-2\mu_kt}c_k^2\le e^{-2\mu_2t}\|\boldsymbol\delta(0)\|^2 ,$$ 所以 $\mu_2$ 给出最慢的衰减速率；当 $\boldsymbol\delta(0)$ 在 $\mathbf v_2$ 上有分量时，长时间行为就由 $e^{-\mu_2t}$ 主导，$1/\mu_2$ 是达成共识的特征时间。由 Courant–Fischer 定理，$\mu_2=\min_{\mathbf x\perp\mathbf 1,\,\mathbf x\ne0}\dfrac{\mathbf x^\top L\mathbf x}{\mathbf x^\top\mathbf x}$：使这个比值小的 $\mathbf x$ 在两个“几乎不相连”的部分上取相反符号，所以存在“瓶颈”的网络 $\mu_2$ 小——这正是 Fiedler 向量可用于二分（见社团结构与社区发现）的原因。对无权网络，Fiedler（1973）还给出上界 $\mu_2\le\frac{N}{N-1}k_{\min}$（$k_{\min}$ 为最小度），即一个度很低的节点就足以让共识变慢。

第 4 步：Fiedler 上界的一行证明。 设 $v$ 为度最小的节点，取 $\mathbf x=\mathbf e_v-\frac1N\mathbf 1$（与 $\mathbf 1$ 正交）。$\mathbf x^\top L\mathbf x=\mathbf e_v^\top L\mathbf e_v=k_{\min}$（$L\mathbf 1=0$），$\mathbf x^\top\mathbf x=1-\frac1N$，代入 Rayleigh 商即得 $\mu_2\le\frac{N}{N-1}k_{\min}$。

C. 主稳定函数与特征值比

假设：$N$ 个相同的 $n$ 维系统，耦合函数 $\mathbf H$ 相同，$L$ 行和为零（保证同步流形 $\mathbf x_1=\cdots=\mathbf x_N=\mathbf s(t)$ 是不变集）且可对角化（无向网络自动满足）；只讨论局部（线性）稳定性。

第 1 步：线性化。 令 $\mathbf x_i=\mathbf s+\boldsymbol\xi_i$，$\dot{\mathbf s}=\mathbf F(\mathbf s)$，保留一阶项： $$\dot{\boldsymbol\xi}_i=D\mathbf F(\mathbf s)\,\boldsymbol\xi_i-\sigma\sum_jL_{ij}\,D\mathbf H(\mathbf s)\,\boldsymbol\xi_j .$$ 堆叠成 $\dot{\boldsymbol\xi}=\big[I_N\otimes D\mathbf F-\sigma\,L\otimes D\mathbf H\big]\boldsymbol\xi$。

第 2 步：解耦。 令 $L=V\Lambda V^\top$，作变换 $\boldsymbol\eta=(V^\top\otimes I_n)\boldsymbol\xi$。利用 Kronecker 积的混合乘积性质 $(V^\top\otimes I_n)(L\otimes D\mathbf H)(V\otimes I_n)=\Lambda\otimes D\mathbf H$，方程分解为 $N$ 个互不耦合的 $n$ 维方程： $$\dot{\boldsymbol\eta}_k=\big[D\mathbf F(\mathbf s)-\sigma\mu_k\,D\mathbf H(\mathbf s)\big]\boldsymbol\eta_k,\qquad k=1,\dots,N .$$ $k=1$（$\mu_1=0$）对应沿同步流形本身的扰动，不影响同步。

第 3 步：主稳定函数。 对实参数 $\alpha$ 定义“主方程”$\dot{\boldsymbol\zeta}=[D\mathbf F(\mathbf s)-\alpha D\mathbf H(\mathbf s)]\boldsymbol\zeta$，其最大李雅普诺夫指数记为 $\Lambda(\alpha)$——它只依赖节点动力学与耦合方式，与网络无关。同步局部稳定当且仅当 $\Lambda(\sigma\mu_k)<0$ 对所有 $k=2,\dots,N$ 成立。若 $\Lambda(\alpha)<0$ 的区域是有限区间 $(\alpha_1,\alpha_2)$，需要 $\sigma\mu_2>\alpha_1$ 且 $\sigma\mu_N<\alpha_2$，存在这样的 $\sigma$ 当且仅当 $\mu_N/\mu_2<\alpha_2/\alpha_1$。若稳定区域是无界的 $(\alpha_1,\infty)$，则只需 $\sigma\mu_2>\alpha_1$，此时只有 $\mu_2$ 起作用。

第 4 步：推广。 有向或非对称耦合时，$L$ 的特征值 $\mu_k$ 可能为复数；若 $L$ 可对角化（$L=V\Lambda V^{-1}$），同样的变换给出 $\dot{\boldsymbol\eta}_k=[D\mathbf F-\sigma\mu_kD\mathbf H]\boldsymbol\eta_k$，只是 $\alpha$ 取复数，主稳定函数成为复平面上的函数 $\Lambda(\alpha_r+i\alpha_i)$，稳定条件变为所有 $\sigma\mu_k$ 落在复平面上的稳定区内。

D. 洛伦兹分布下的 Ott–Antonsen 约化

Ott 与 Antonsen（2008，Chaos）发现，对 $N\to\infty$ 的全连接 Kuramoto 模型，振荡器相位密度 $f(\theta,\omega,t)$ 存在一个不变流形，其上傅里叶系数满足 $f_n=\alpha(\omega,t)^n$（几何级数）。在这一流形上，连续性方程化为关于 $\alpha$ 的单个方程： $$\dot\alpha=-i\omega\alpha+\frac K2\big(z^*-z\,\alpha^2\big),\qquad z^*=\int g(\omega)\,\alpha(\omega,t)\,d\omega ,$$ 其中 $z=re^{i\psi}$ 为序参量。对洛伦兹分布，积分可以用留数定理计算：把 $\alpha$ 解析延拓到 $\omega$ 的下半平面，只需在极点 $\omega=-i\gamma$ 处取值，得到 $z^*=\alpha(-i\gamma,t)$。代入上式（取 $\omega=-i\gamma$）并分离模长，得到关于 $r$ 的一维常微分方程： $$\dot r=-\gamma r+\frac K2\,r\,(1-r^2).$$ 它的不动点为 $r=0$ 与 $r^*=\sqrt{1-2\gamma/K}$（$K>2\gamma$），后者稳定。这与 A 节的自洽结果完全一致，但额外给出了动力学：$r(t)$ 如何随时间趋近稳态，以及 $r=0$ 在 $K<2\gamma$ 时稳定、在 $K>2\gamma$ 时失稳。Ott–Antonsen 约化是分析爆发式同步、奇异态（两组振荡器各自一个 $\alpha$）以及带外部驱动的 Kuramoto 模型的标准工具；它对带噪声的系统不严格成立。

E. 簇同步的商网络

设节点划分 $\{\mathcal C_1,\dots,\mathcal C_q\}$ 是等价划分：对任意 $a,b$，$\mathcal C_a$ 中每个节点连向 $\mathcal C_b$ 的边权之和都等于同一个数 $Q_{ab}$。令 $\mathbf x_i=\mathbf y_a$（$i\in\mathcal C_a$），代入 $\dot{\mathbf x}_i=\mathbf F(\mathbf x_i)+\sigma\sum_jA_{ij}\mathbf H(\mathbf x_j)$（这里用邻接形式，拉普拉斯形式同理）： $$\dot{\mathbf y}_a=\mathbf F(\mathbf y_a)+\sigma\sum_{b=1}^{q}Q_{ab}\,\mathbf H(\mathbf y_b).$$ 右端对 $\mathcal C_a$ 内所有节点相同，所以“簇内相同”的子空间是不变的；簇同步的动力学由 $q$ 个节点的商网络 $Q$ 描述。自同构轨道一定构成等价划分，反之不然。簇同步是否稳定，还需要对横向扰动做类似 C 节的线性化，Pecora 等（2014）用对称群的不可约表示把横向扰动分块，从而降低计算量。

算法

算法 1  网络 Kuramoto 模拟与 K 扫描（含前向/后向扫描）
输入：A、ω、K 网格、步长 h、预热步数 T0、统计步数 T1
1. θ ← 均匀随机初值
2. for K in K 网格（前向：递增；后向：递减，从强同步态开始）:
3.     （沿用上一 K 的终态 θ，以观察滞后；若只要无滞后的曲线，可每次重新随机初值）
4.     for t = 1 … T0+T1:
5.         f(θ)_i = ω_i + K Σ_j A_ij sin(θ_j − θ_i)   # 稀疏：cosθ_i (A sinθ)_i − sinθ_i (A cosθ)_i
6.         用 RK4（或小步长欧拉）以步长 h 推进 θ
7.         if t > T0: 记录 r(t) = |Σ_j e^{iθ_j}| / N
8.     输出 K 与时间平均 ⟨r⟩（以及标准差，作为亚稳态度量）

复杂度：每次求右端为 $O(m)$（全连接时利用 $\sum_j\sin(\theta_j-\theta_i)=N r\sin(\psi-\theta_i)$ 可降为 $O(N)$），每个 $K$ 为 $O\big((T_0+T_1)\,m\big)$。判断同步时与序参量的 $1/\sqrt N$ 基线比较；判断一级相变时比较前向与后向曲线（见爆发式同步）。

算法 2  数值计算主稳定函数 Λ(α)
输入：F、DF、DH、α 网格、步长 h、积分时长、重归一化间隔 τ
1. 积分 ds/dt = F(s) 越过暂态，得到同步轨道上的初值 s0
2. for α in α 网格:
3.     s ← s0；ζ ← 随机单位向量；S ← 0
4.     同时积分 ds/dt = F(s) 与 dζ/dt = [DF(s) − α DH(s)] ζ
5.     每隔 τ：S ← S + ln‖ζ‖；ζ ← ζ/‖ζ‖
6.     Λ(α) ← S / 总积分时间
7. 在 Λ 变号处用二分法细化 α1、α2；对给定网络计算 μ2、μN 并检验 μN/μ2 < α2/α1

复杂度：MSF 只需对每个 $\alpha$ 积分一次 $n$ 维系统（$n$ 通常很小），与网络规模无关；网络部分只需求 $L$ 的两个极端特征值，稀疏网络用 Lanczos 方法即可。李雅普诺夫指数的估计误差随积分时长约按 $1/\sqrt T$ 减小，$\Lambda$ 接近 0 处的边界 $\alpha_1,\alpha_2$ 最不确定。

算法 3  从多通道信号估计相位同步网络
输入：N 通道信号、频段 [f1, f2]、指标（PLV / PPC / wPLI …）、替代数据次数 S
1. 零相位带通滤波（如双向 Butterworth），去掉边缘若干周期
2. 希尔伯特变换得到 X_i(t) = A_i(t) e^{iφ_i(t)}
3. for 每对 (i, j):
4.     S_ij(t) = X_i(t) X_j(t)*；按定义计算指标 W_ij
5.     for s = 1 … S: 对 X_j 做随机循环平移（或相位随机化），重算指标，得零分布
6.     p_ij = 零分布中 ≥ W_ij 的比例
7. 多重比较校正（FDR 等）后得到网络；报告频段、数据长度、指标与参考方式

复杂度：$O(N^2T\,S)$。替代数据只能检验“两信号独立”这一零假设，不能排除共同源或容积传导（见从实测信号估计同步的算例）。