Учебник: седловой цикл Мёбиуса - удвоение периода и отрицательные собственные значения#

Этот ноутбук демонстрирует топологически интересную неподвижную точку: седловой цикл Мёбиуса, где матрица монодромии DPm имеет отрицательные собственные значения (\(\lambda_u = -e^{+1/3}\), \(\lambda_s = -e^{-1/3}\)).

Физический смысл:
После одного полного тороидального периода (\(m=1\) оборот) малое смещение меняет знак при растяжении/сжатии. Оно возвращается к исходной ориентации только после двух оборотов - это топология, похожая на ленту Мёбиуса в локальном касательном расслоении.

Это отличается от обычной X-точки (положительное \(\lambda\)) и иногда называется обратно-гиперболической или удвоенной по периоду неподвижной точкой.

Что вы узнаете:

  1. Как седловые точки с отрицательными собственными значениями возникают в картах Пуанкаре

  2. Как использовать evolve_DPm_along_cycle, чтобы подтвердить знак собственных значений

  3. Как визуализировать геометрию устойчивого/неустойчивого многообразия в 3D: они образуют полуперекрученную ленту (ленту Мёбиуса) вокруг тора

1. Построение поля: неподвижная точка Мёбиуса m=1 в (R₀, Z₀)#

Простейший случай: одна неподвижная точка однооборотной карты в \((R_0, Z_0)\) с собственными значениями \(\lambda_u = -e^{+1/3}\), \(\lambda_s = -e^{-1/3}\).

Собственные направления поворачиваются на \(\pi\) за каждый оборот (dθ/dφ = 1/2), поэтому после одного оборота \(\theta \to \theta + \pi\) собственные векторы меняют знак, что и дает характер Мёбиуса.

[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}  (должно быть +1 для сохраняющей площадь карты)")
print(f"После 2 оборотов: λ_u² = {lam_u**2:.6f}  (теперь положительно, обычная гиперболика)")

# 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"Разность = {np.degrees(theta_u(2*np.pi) - theta_u(0)):.1f}°  (= 180° переворот -> Мебиус)")

λ_u = -1.395612
λ_s = -0.716531
λ_u · λ_s = 1.0000000000  (должно быть +1 для сохраняющей площадь карты)
После 2 оборотов: λ_u² = 1.947734  (теперь положительно, обычная гиперболика)

θ_u(0) = 1.5708 рад = 90.0°
θ_u(2π) = 4.7124 рад = 270.0°
Разность = 180.0°  (= переворот на 180° -> Мебиус)

[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 в неподвижной точке (R0,Z0):", g0, "  (должно быть 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("Поле в соседней точке:")
print(pyna_field(np.array([R0 + 0.05, Z0, 0.0])))

g в неподвижной точке (R0,Z0): [0. 0.]   (должно быть 0)
Поле в соседней точке:
[-0.00252555  0.0238027   0.95210808]

2. Вычисление DPm с evolve_DPm_along_cycle#

[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_вывод=2 * np.pi / 500,
    rtol=1e-11, atol=1e-13,
)

print("DPm (матрица монодромии):")
print(mono.DPm)
print()
print("Собственные значения:", mono.eigenvalues)
print("Ожидается:   ", np.array([lam_u, lam_s]))
print()
print("Оба собственных значения ОТРИЦАТЕЛЬНЫ - это характер седла Мебиуса.")
print("Индекс устойчивости Tr(DPm)/2 =", mono.stability_index,
      f" (|k| = {abs(mono.stability_index):.4f} > 1 → hyperbolic)")

DPm (матрица монодромии):
[[-7.16531311e-01 -1.74438498e-09]
 [ 1.74438440e-09 -1.39561242e+00]]

Собственные значения: [-0.71653131 -1.39561242]
Ожидается:           [-1.39561243 -0.71653131]

Оба собственных значения ОТРИЦАТЕЛЬНЫ - это характер седла Мебиуса.
Индекс устойчивости Tr(DPm)/2 = -1.0560718675230398  (|k| = 1.0561 > 1 -> гиперболическая)

[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:")
print(DPm_analytic)
print()
print("Макс. ошибка:", np.max(np.abs(mono.DPm - DPm_analytic)))

Аналитическая DPm:
[[ 7.16531311e-01  2.12494955e-16]
 [-8.77497776e-17  1.39561243e+00]]

Макс. ошибка: 2.7912248489659737

3. Визуализация перекручивания Мёбиуса в геометрии многообразий#

Устойчивое многообразие седла Мёбиуса является полуперекрученной лентой в 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('Угол собственного вектора (градусы)')
axes[0].set_title('Поворот собственного направления: пол-оборота за орбиту')
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: значения при φ=2π (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 [м]')
ax.set_title('Седло Мебиуса: полосы устойчивого (синего) и неустойчивого (красного) многообразий')
ax.legend()
plt.tight_layout()
plt.savefig("mobius_manifold_3d.png", dpi=120)
plt.show()

print()
print("Обратите внимание: каждая полоса делает ПОЛУПЕРЕКРУТ за 2π.")
print("Лента замыкается (возвращается к той же ориентации) только после ДВУХ оборотов.")


Обратите внимание: каждая полоса делает ПОЛУПЕРЕКРУТ за 2π.
Лента замыкается (возвращается к той же ориентации) только после ДВУХ оборотов.

4. Двухоборотная монодромия: возврат к обычной гиперболике#

Если вычислить DPm за 2 оборота (m=2), отрицательные знаки сокращаются, и получается обычная гиперболическая неподвижная точка с положительными собственными значениями.

[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_вывод=2 * np.pi / 500,
    rtol=1e-11, atol=1e-13,
)

print("Собственные значения DPm за 1 оборот:", mono.eigenvalues,   "  (отрицательные)")
print("Собственные значения DPm за 2 оборота:", mono_2.eigenvalues, "  (положительные, = λ²)")
print()
print("Ожидаемые собственные значения за 2 оборота:", np.array([lam_u**2, lam_s**2]))
print("Совпадение:", np.allclose(sorted(np.abs(mono_2.eigenvalues)),
                             sorted([lam_u**2, lam_s**2]), rtol=1e-4))

Собственные значения DPm за 1 оборот: [-0.71653131 -1.39561242]   (отрицательные)
Собственные значения DPm за 2 оборота: [0.51341713 1.94773399]   (положительные, = λ²)

Ожидаемые собственные значения за 2 оборота: [1.94773404 0.51341712]
Совпадение: True

5. Итоги#

Свойство

Значение

Неподвижная точка

\((R_0, Z_0) = (1.0, 0.0)\)

Период орбиты

\(m = 1\) оборот

Собственные значения DPm

\(\lambda_u = -e^{+1/3}\), \(\lambda_s = -e^{-1/3}\)

Характер

седло Мёбиуса (обратно-гиперболическое)

Поворот собственного вектора

\(\Delta\theta = \pi\) за оборот (полуповорот)

Двухоборотная карта

обычная гиперболическая с \(\lambda^2 = e^{\pm 2/3}\)
Ключевая идея: знак собственных значений DPm кодирует топологический тип инвариантных многообразий. Одна буква M скрывает это;
DPm явно показывает, что это производная полнооборотной карты Пуанкаре, чьи спектральные свойства определяют локальную топологию.

Дальше:

  • См. monodromy_xcycle_analytic.ipynb для случая аналитического седлового цикла

  • См. notebooks/research/ для шаблонов публичных аналитических и экспериментальных приложений