Möbius 鞍点循环：倍周期与负特征值#
这个 notebook 演示一种拓扑上有意思的不动点：Möbius 鞍点循环。在这种情形下，Monodromy 矩阵 DPm 具有负特征值（\(\lambda_u = -e^{+1/3}\)，\(\lambda_s = -e^{-1/3}\)）。
物理含义：
经过一个完整环向周期（\(m=1\) 圈）后，一个小位移在伸长/收缩的同时会翻转符号。它只有在两圈之后才回到原始取向，这对应局部切丛中类似 Möbius 带的拓扑。
这不同于普通 X 点（正 \(\lambda\)），有时称为反双曲型或倍周期不动点。
你将学到：
负特征值鞍点如何出现在 Poincaré 映射中
如何使用
evolve_DPm_along_cycle确认特征值符号如何在 3D 中可视化稳定/不稳定流形几何：它们在环面周围形成一条半扭转带（Möbius 带）
1. 场构造：位于 \((R₀, Z₀)\) 的 \(m=1\) Möbius 不动点#
最简单的情形：一圈映射在 \((R_0, Z_0)\) 有一个不动点，特征值为 \(\lambda_u = -e^{+1/3}\)、\(\lambda_s = -e^{-1/3}\)。
特征方向每圈旋转 \(\pi\)（dθ/dφ = 1/2），因此一圈之后 \(\theta \to \theta + \pi\)：特征向量翻转，这就是 Möbius 特征的来源。
[1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
# ── Fixed point and field parameters ────────────────────────────────────────
R0, Z0 = 1.0, 0.0
BPhi0 = 2.5
m = 1 # one-turn fixed point
# Negative eigenvalues → Möbius / inverse-hyperbolic character
lam_u = -np.e ** (1/3)
lam_s = -np.e ** (-1/3)
print(f"λ_u = {lam_u:.6f}")
print(f"λ_s = {lam_s:.6f}")
print(f"λ_u · λ_s = {lam_u * lam_s:.10f} (should be +1 for area-preserving map)")
print(f"After 2 turns: λ_u² = {lam_u**2:.6f} (now positive, normal hyperbolic)")
# Eigendirections rotate by π per turn (half-integer winding)
dtheta_dphi = 0.5
def theta_u(phi):
return np.pi/2 + phi/2
def theta_s(phi):
return phi/2
# Note: theta_u(2π) = π/2 + π = 3π/2 ≠ theta_u(0)
# The eigenvectors have flipped! This is the Möbius twist.
print(f"\nθ_u(0) = {theta_u(0):.4f} rad = {np.degrees(theta_u(0)):.1f}°")
print(f"θ_u(2π) = {theta_u(2*np.pi):.4f} rad = {np.degrees(theta_u(2*np.pi)):.1f}°")
print(f"Difference = {np.degrees(theta_u(2*np.pi) - theta_u(0)):.1f}° (= 180° flip → Möbius)")
λ_u = -1.395612
λ_s = -0.716531
λ_u · λ_s = 1.0000000000 (should be +1 for area-preserving map)
After 2 turns: λ_u² = 1.947734 (now positive, normal hyperbolic)
θ_u(0) = 1.5708 rad = 90.0°
θ_u(2π) = 4.7124 rad = 270.0°
Difference = 180.0° (= 180° flip → Möbius)
[2]:
# ── Toroidal field model ─────────────────────────────────────────────────────
def BPhi(phi, RZ):
"""Toroidal field: R·B_φ = const."""
return BPhi0 * R0 / RZ[0]
# ── A matrix at the fixed point ──────────────────────────────────────────────
def _eigvec_matrix(phi):
return np.array([
[np.cos(theta_u(phi)), np.cos(theta_s(phi))],
[np.sin(theta_u(phi)), np.sin(theta_s(phi))],
])
def _A_matrix(phi):
"""Jacobian A(φ) of field direction g w.r.t. (R,Z) at the fixed point."""
V = _eigvec_matrix(phi)
# log(|λ|) / (2mπ) gives the per-radian growth rate
mu_u = np.log(np.abs(lam_u)) / (2 * m * np.pi)
mu_s = np.log(np.abs(lam_s)) / (2 * m * np.pi)
Lam = np.diag([mu_u, mu_s])
dV = np.array([
[-dtheta_dphi * np.sin(theta_u(phi)), -dtheta_dphi * np.sin(theta_s(phi))],
[ dtheta_dphi * np.cos(theta_u(phi)), dtheta_dphi * np.cos(theta_s(phi))],
])
return V @ Lam @ np.linalg.inv(V) + dV @ np.linalg.inv(V)
# ── φ-parameterised field direction g = R·B_pol/B_φ ─────────────────────────
def field_g_pol(phi, X_pol):
"""Field direction in (R,Z): g = A(φ) · (X - X_c) + 0 (fixed point has dX_c/dφ=0)."""
Xc = np.array([R0, Z0])
A = _A_matrix(phi)
return A @ (X_pol - Xc)
# Verify: g on the fixed point is zero (it's a FIXED point, not just periodic orbit)
g0 = field_g_pol(0.5, np.array([R0, Z0]))
print("g at fixed point (R0,Z0):", g0, " (should be 0)")
# ── pyna-compatible field function ───────────────────────────────────────────
def pyna_field(rzphi):
R, Z, phi = rzphi[0], rzphi[1], rzphi[2]
X_pol = np.array([R, Z])
g = field_g_pol(phi, X_pol)
# Regularise near the fixed point to avoid division by zero
gnorm2 = np.dot(g, g)
dphi_dl = 1.0 / np.sqrt(R**2 + max(gnorm2, 1e-30))
return np.array([g[0] * dphi_dl, g[1] * dphi_dl, dphi_dl])
print("Field at a nearby point:")
print(pyna_field(np.array([R0 + 0.05, Z0, 0.0])))
g at fixed point (R0,Z0): [0. 0.] (should be 0)
Field at a nearby point:
[-0.00252555 0.0238027 0.95210808]
2. 使用 evolve_DPm_along_cycle 计算 DPm#
[3]:
from pyna.topo.monodromy import evolve_DPm_along_cycle
from pyna.topo.toroidal_cycle import ToroidalPeriodicOrbitTrace as PeriodicOrbit
orbit = PeriodicOrbit(
rzphi0=np.array([R0, Z0, 0.0]),
period_m=m,
trajectory=np.zeros((1, 3)),
DPm=np.eye(2),
)
mono = evolve_DPm_along_cycle(
field_func=pyna_field,
orbit=orbit,
dt_output=2 * np.pi / 500,
rtol=1e-11, atol=1e-13,
)
print("DPm (monodromy matrix):")
print(mono.DPm)
print()
print("Eigenvalues:", mono.eigenvalues)
print("Expected: ", np.array([lam_u, lam_s]))
print()
print("Both eigenvalues are NEGATIVE — this is the Möbius saddle character.")
print("Stability index Tr(DPm)/2 =", mono.stability_index,
f" (|k| = {abs(mono.stability_index):.4f} > 1 → hyperbolic)")
DPm (monodromy matrix):
[[-7.16531311e-01 -1.74438498e-09]
[ 1.74438440e-09 -1.39561242e+00]]
Eigenvalues: [-0.71653131 -1.39561242]
Expected: [-1.39561243 -0.71653131]
Both eigenvalues are NEGATIVE — this is the Möbius saddle character.
Stability index Tr(DPm)/2 = -1.0560718675230398 (|k| = 1.0561 > 1 → hyperbolic)
[4]:
# ── Analytic check ───────────────────────────────────────────────────────────
# DPm_analytic = V(2π) · diag(λ_u, λ_s) · V^{-1}(0)
V_end = _eigvec_matrix(2 * m * np.pi)
V_start = _eigvec_matrix(0.0)
DPm_analytic = V_end @ np.diag([lam_u, lam_s]) @ np.linalg.inv(V_start)
print("DPm analytic:")
print(DPm_analytic)
print()
print("Max error:", np.max(np.abs(mono.DPm - DPm_analytic)))
DPm analytic:
[[ 7.16531311e-01 2.12494955e-16]
[-8.77497776e-17 1.39561243e+00]]
Max error: 2.7912248489659737
3. 在流形几何中可视化 Möbius 扭转#
Möbius 鞍点的稳定流形在 3D 中是一条半扭转带：它绕环面一圈，但取向发生翻转。 经过两圈后，它变成普通带状结构。
我们沿轨道在每个 φ 处绘制局部稳定/不稳定方向 e_s(φ)、e_u(φ)，以观察半圈旋转。
[5]:
phi_plot = np.linspace(0, 2 * np.pi, 400)
# Eigendirection angles along one turn
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].plot(phi_plot / np.pi, np.degrees(theta_u(phi_plot)) % 360, label='θ_u (unstable)', color='tomato')
axes[0].plot(phi_plot / np.pi, np.degrees(theta_s(phi_plot)) % 360, label='θ_s (stable)', color='steelblue')
axes[0].set_xlabel('φ / π')
axes[0].set_ylabel('Eigenvector angle (degrees)')
axes[0].set_title('Eigendirection rotation: half-turn per orbit')
axes[0].legend()
axes[0].axvline(2, color='gray', linestyle='--', label='End of period')
# DX_pol components showing sign flip at φ=2π
axes[1].plot(mono.phi_arr / np.pi, mono.DX_pol_arr[:, 0, 0], label='DX_pol[0,0]', color='steelblue')
axes[1].plot(mono.phi_arr / np.pi, mono.DX_pol_arr[:, 1, 1], label='DX_pol[1,1]', color='tomato')
axes[1].axhline(0, color='k', lw=0.5)
axes[1].axvline(2, color='gray', linestyle='--')
axes[1].set_xlabel('φ / π')
axes[1].set_title('DX_pol diagonal: note values at φ=2π (the DPm)')
axes[1].legend()
plt.tight_layout()
plt.savefig("mobius_DX_pol.png", dpi=120)
plt.show()
[6]:
# ── 3D plot: the stable manifold band (Möbius strip) ─────────────────────────
from mpl_toolkits.mplot3d import Axes3D
nS = 30 # strip width samples
nPhi = 400
Phi = np.linspace(0, 2 * np.pi, nPhi, endpoint=True)
s_vals = np.linspace(-0.25, 0.25, nS) # displacement along stable direction
# Stable manifold: X_c + s · e_s(φ)
theta_s_arr = theta_s(Phi) # angle of stable eigendirection at each φ
e_s_R = np.cos(theta_s_arr) # (nPhi,)
e_s_Z = np.sin(theta_s_arr)
# Grid: [iS, iPhi]
R_strip = R0 + s_vals[:, None] * e_s_R[None, :]
Z_strip = Z0 + s_vals[:, None] * e_s_Z[None, :]
x_strip = R_strip * np.cos(Phi)[None, :]
y_strip = R_strip * np.sin(Phi)[None, :]
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x_strip, y_strip, Z_strip,
alpha=0.5, color='steelblue', linewidth=0)
# The fixed point trajectory (a circle in the φ=0 plane → a point)
ax.scatter([R0], [0], [Z0], color='red', s=80, zorder=5, label='Fixed point')
# Unstable manifold band
theta_u_arr = theta_u(Phi)
e_u_R = np.cos(theta_u_arr)
e_u_Z = np.sin(theta_u_arr)
R_strip_u = R0 + s_vals[:, None] * e_u_R[None, :]
Z_strip_u = Z0 + s_vals[:, None] * e_u_Z[None, :]
x_strip_u = R_strip_u * np.cos(Phi)[None, :]
y_strip_u = R_strip_u * np.sin(Phi)[None, :]
ax.plot_surface(x_strip_u, y_strip_u, Z_strip_u,
alpha=0.3, color='tomato', linewidth=0)
ax.set_xlabel('x [m]'); ax.set_ylabel('y [m]'); ax.set_zlabel('Z [m]')
ax.set_title('Möbius saddle: stable (blue) and unstable (red) manifold bands')
ax.legend()
plt.tight_layout()
plt.savefig("mobius_manifold_3d.png", dpi=120)
plt.show()
print()
print("Notice: each strip makes a HALF-TWIST over 2π.")
print("The band only closes (returns to same orientation) after TWO turns.")
Notice: each strip makes a HALF-TWIST over 2π.
The band only closes (returns to same orientation) after TWO turns.
4. 两圈 Monodromy：回到普通双曲型#
如果对两圈（m=2）计算 DPm，负号会相互抵消，从而恢复具有正特征值的普通双曲型不动点。
[7]:
orbit_2turn = PeriodicOrbit(
rzphi0=np.array([R0, Z0, 0.0]),
period_m=2, # two turns
trajectory=np.zeros((1, 3)),
DPm=np.eye(2),
)
mono_2 = evolve_DPm_along_cycle(
field_func=pyna_field,
orbit=orbit_2turn,
dt_output=2 * np.pi / 500,
rtol=1e-11, atol=1e-13,
)
print("1-turn DPm eigenvalues:", mono.eigenvalues, " (negative)")
print("2-turn DPm eigenvalues:", mono_2.eigenvalues, " (positive, = λ²)")
print()
print("Expected 2-turn eigenvalues:", np.array([lam_u**2, lam_s**2]))
print("Match:", np.allclose(sorted(np.abs(mono_2.eigenvalues)),
sorted([lam_u**2, lam_s**2]), rtol=1e-4))
1-turn DPm eigenvalues: [-0.71653131 -1.39561242] (negative)
2-turn DPm eigenvalues: [0.51341713 1.94773399] (positive, = λ²)
Expected 2-turn eigenvalues: [1.94773404 0.51341712]
Match: True
5. 总结#
属性 |
数值 |
|---|---|
不动点 |
\((R_0, Z_0) = (1.0, 0.0)\) |
轨道周期 |
\(m = 1\) 圈 |
DPm 特征值 |
\(\lambda_u = -e^{+1/3}\)，\(\lambda_s = -e^{-1/3}\) |
类型 |
Möbius 鞍点（反双曲型） |
特征向量旋转 |
每圈 \(\Delta\theta = \pi\)（半扭转） |
两圈映射 |
普通双曲型，\(\lambda^2 = e^{\pm 2/3}\) |
关键洞见： DPm 特征值的符号编码了不变流形的拓扑类型。单个字母 M 会遮蔽这一点；DPm 明确指出它是整周期 Poincaré 映射的导数，其谱性质决定局部拓扑。
下一步：
参见
monodromy_xcycle_analytic.ipynb，了解 X 型周期轨道情形参见
pyna/notebooks/research/，了解 private stellarator 和其他真实装置应用