Учебник: матрица монодромии и сдвиг орбиты на аналитическом седловом X-цикле#

Этот ноутбук демонстрирует основной рабочий процесс монодромии pyna на аналитически построенном седловом цикле (X-точке карты Пуанкаре), вложенном в 3D тороидальную геометрию.

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

  1. Как задать аналитическое поле с известным X-циклом m=3/n=1 в цилиндрических координатах (R, Z, φ)

  2. Как evolve_DPm_along_cycle интегрирует вариационное уравнение dDX_pol/dφ = A · DX_pol, чтобы получить матрицу монодромии DPm

  3. Почему эволюции DPm по уравнению коммутатора нужно DPm(φ₀) = DX_pol(φ_end), а не единичная матрица, и что это значит геометрически

  4. Как вычислить линейный сдвиг орбиты под действием поля возмущения с помощью orbit_shift_under_perturbation

  5. Как сравнить аналитическую формулу с результатом API pyna

Опорная геометрия:
Модель токамака с большим радиусом R₀=1 м, B_φ^ax=2.5 Тл, эллиптическим X-циклом с полуосями (R_ell=0.3, Z_ell=0.5) и вращательным преобразованием ι = n/m = 1/3.

1. Построение аналитического поля#

Строится поле, чей точный X-цикл m=3 известен аналитически.

1.1 Цикл и его собственные направления#

Цикл лежит на эллипсе

\[X_c(\varphi) = \begin{pmatrix} R_{\rm ax} + R_{\rm ell}\cos(\iota\varphi + \varphi_0) \ Z_{\rm ax} + Z_{\rm ell}\sin(\iota\varphi + \varphi_0) \end{pmatrix}\]

с \(\iota = 1/3\), поэтому орбита замыкается после \(m=3\) тороидальных оборотов.

Карта Пуанкаре в X-точке имеет собственные значения \(\lambda_u = e^{+1/5}\) (неустойчивое) и \(\lambda_s = e^{-1/5}\) (устойчивое). Собственные направления вращаются вдоль орбиты по заданному углу \(\theta(\varphi)\).

[1]:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# ── Cycle parameters ────────────────────────────────────────────────────────
Rax, Zax = 1.0, 0.0
Rell, Zell = 0.3, 0.5
phi0_cycle = 0.0
BPhiax = 2.5
m = 3       # toroidal period (orbit closes after m turns)
n = 1
iota = n / m  # rotational transform

# Собственные значения of the m-turn Poincaré map
lam_u = np.e ** (1/5)   # unstable multiplier
lam_s = np.e ** (-1/5)  # stable  multiplier
print(f"λ_u = {lam_u:.6f},  λ_s = {lam_s:.6f},  произведение = {lam_u * lam_s:.10f}  (должно быть 1)")

# Rotating eigendirections along the orbit
def theta_u(phi):
    """Angle of unstable eigenvector at toroidal position phi."""
    return phi/3 + np.pi/2 + np.pi/9

def theta_s(phi):
    """Angle of stable eigenvector at toroidal position phi."""
    return phi/3 + np.pi/2 - np.pi/9

dtheta_dphi = 1.0/3  # same rotation rate as the orbit

# ── The cycle position ──────────────────────────────────────────────────────
def cycle_RZ(phi):
    """Return (R, Z) of the X-cycle at toroidal angle phi."""
    return np.array([
        Rell * np.cos(iota * phi + phi0_cycle) + Rax,
        Zell * np.sin(iota * phi + phi0_cycle) + Zax,
    ])

# ── Toroidal field (axisymmetric model: R·B_φ = const) ─────────────────────
def BPhi(phi, RZ):
    """Toroidal field B_φ(R) = B_φ^ax · R_ax / R  (free-force-balance model)."""
    return BPhiax * Rax / RZ[0]

print("Цикл при phi=0:", cycle_RZ(0.0))
print("Цикл при phi=2π:", cycle_RZ(2*np.pi))
print("Цикл при phi=6π (замыкается):", cycle_RZ(6*np.pi))

λ_u = 1.221403,  λ_s = 0.818731,  произведение = 1.0000000000  (должно быть 1)
Цикл при phi=0: [1.3 0. ]
Цикл при phi=2π: [0.85      0.4330127]
Цикл при phi=6π (замыкается): [ 1.3000000e+00 -1.2246468e-16]

[2]:

# ── Poloidal field on the cycle ──────────────────────────────────────────────
#
# The poloidal field on the exact cycle is just B_pol = (dR_c/dl, dZ_c/dl)
# re-parameterised by phi.  We compute the φ-derivatives:
#   dR_c/dφ = -ι R_ell sin(ι φ + φ₀)
#   dZ_c/dφ = +ι Z_ell cos(ι φ + φ₀)

def BR0_BZ0_on_cycle(phi):
    """Return (B_R, B_Z) of the unperturbed field on the X-cycle."""
    Rc = cycle_RZ(phi)[0]
    Bp = BPhi(phi, cycle_RZ(phi))
    dRc = -iota * Rell * np.sin(iota * phi + phi0_cycle)
    dZc =  iota * Zell * np.cos(iota * phi + phi0_cycle)
    return np.array([dRc / Rc * Bp, dZc / Rc * Bp])


# ── Full field B(R, Z, φ) with exact X-cycle ────────────────────────────────
#
# We construct B_pol(X, φ) so that:
#  1. B_pol(X_c(φ), φ) = B_pol^0(φ)  (matches cycle tangent)
#  2. The Jacobian A(φ) of g = R·B_pol/B_φ along X_c has the prescribed
#     eigenstructure (λ_u, λ_s with rotating eigenvectors θ_u, θ_s).
#
# Construction: let V(φ) = [e_u | e_s] be the eigenvector matrix and
# Λ = diag(log|λ_u|, log|λ_s|) / (2mπ)  (per-φ growth rates).
# Then  A(φ) = V Λ V^{-1} + dV/dφ · V^{-1}  (accounts for rotating frame).
# The field is B_pol(X, φ) = B_pol^0(φ) + R·A(φ)·(X - X_c(φ)) / R_c

def _eigvec_matrix(phi):
    """2x2 matrix of eigenvectors [e_u | e_s] at 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):
    """A(φ) = Jacobian of g = R·B_pol/B_φ w.r.t. (R,Z) along cycle."""
    V = _eigvec_matrix(phi)
    # per-radian growth rates
    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])
    # rotating-frame correction: dV/dφ · V^{-1}
    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)

def field_g_pol(phi, X_pol):
    """φ-parameterised field direction g = R·B_pol/B_φ  (units: m)."""
    Xc = cycle_RZ(phi)
    g0 = Xc  # Xc gives the cycle tangent when dividing by R·B_phi/R_B_phi^ax
    # Actually g = dX_c/dphi + A·(X - X_c)
    dXc = np.array([
        -iota * Rell * np.sin(iota * phi + phi0_cycle),
         iota * Zell * np.cos(iota * phi + phi0_cycle),
    ])
    A = _A_matrix(phi)
    return dXc + A @ (X_pol - Xc)

# Verify: g on the cycle equals dX_c/dphi
phi_test = 0.7
Xc_test = cycle_RZ(phi_test)
g_on_cycle = field_g_pol(phi_test, Xc_test)
dXc_expected = np.array([
    -iota * Rell * np.sin(iota * phi_test + phi0_cycle),
     iota * Zell * np.cos(iota * phi_test + phi0_cycle),
])
print("g на цикле  :", g_on_cycle)
print("dX_c/dφ     :", dXc_expected)
print("Совпадение:", np.allclose(g_on_cycle, dXc_expected))

g на цикле  : [-0.02312218  0.16215018]
dX_c/dφ     : [-0.02312218  0.16215018]
Совпадение: True

2. Обертка как функция поля pyna#

pyna ожидает функции поля с сигнатурой field_func(rzphi) -> (dR/dl, dZ/dl, dφ/dl), то есть параметризованные длиной дуги. Мы преобразуем из φ-параметризации с помощью \(d\varphi/dl = B_\varphi / (R |B|)\).

[3]:

def pyna_field(rzphi):
    """pyna-compatible field function: rzphi=(R,Z,φ) ?(dR/dl, dZ/dl, dφ/dl)."""
    R, Z, phi = rzphi[0], rzphi[1], rzphi[2]
    X_pol = np.array([R, Z])
    g = field_g_pol(phi, X_pol)          # (dR/dφ, dZ/dφ)
    Bp = BPhi(phi, X_pol)                # B_φ at (R, Z)
    # dφ/dl = B_φ / (R · |B|),  |B|² = B_R² + B_Z² + B_φ²
    # B_R = g[0] · B_φ / R,  B_Z = g[1] · B_φ / R
    # So  |B|² = B_φ² · (g²/R² + 1)  ? dphi/dl = 1/sqrt(R² + |g|²)
    dphi_dl = 1.0 / np.sqrt(R**2 + np.dot(g, g))
    dR_dl = g[0] * dphi_dl
    dZ_dl = g[1] * dphi_dl
    return np.array([dR_dl, dZ_dl, dphi_dl])

# Quick sanity: integrate the cycle once and check it closes
def rhs_phi(phi, XZ):
    f = pyna_field(np.array([XZ[0], XZ[1], phi]))
    dphi_dl = f[2]
    return np.array([f[0] / dphi_dl, f[1] / dphi_dl])  # dX/dphi

X0 = cycle_RZ(0.0)
sol = solve_ivp(rhs_phi, [0, 2 * m * np.pi], X0, rtol=1e-10, atol=1e-12, max_step=0.01)
X_end = sol.y[:, -1]
print(f"Начало цикла: {X0}")
print(f"Конец цикла  : {X_end}")
print(f"Ошибка замыкания: {np.linalg.norm(X_end - X0):.2e}  (должна быть < 1e-6)")

Начало цикла: [1.3 0. ]
Конец цикла : [1.30000000e+00 3.20923843e-17]
Ошибка замыкания: 4.45e-16  (должна быть < 1e-6)

3. Вычисление матрицы монодромии с pyna.topo.evolve_DPm_along_cycle#

Функция evolve_DPm_along_cycle интегрирует две связанные системы вдоль орбиты:

Фаза 1 - вариационное уравнение для DX_pol:

\[\frac{d\,\mathtt{DX\_pol}}{d\varphi} = A(\varphi)\cdot \mathtt{DX\_pol}, \quad \mathtt{DX\_pol}(\varphi_0) = I\]

После одного полного периода (\(\varphi_0 \to \varphi_0 + 2m\pi\)) получаем \(\mathtt{DPm} = \mathtt{DX\_pol}(\varphi_{\rm end})\) - матрицу монодромии.

Фаза 2 - уравнение коммутатора для эволюции DPm:

\[\frac{d\,\mathtt{DPm}}{d\varphi} = A\,\mathtt{DPm} - \mathtt{DPm}\,A, \quad \mathtt{DPm}(\varphi_0) = \mathtt{DX\_pol}(\varphi_{\rm end})\]

Обратите внимание: начальное условие не является единичной матрицей - это матрица монодромии, только что вычисленная в фазе 1. Это гарантирует, что DPm(φ) отслеживает, как выглядел бы якобиан полнооборотной карты, если начинать период в φ, а не в φ₀.

[4]:

from pyna.topo.monodromy import evolve_DPm_along_cycle
from pyna.topo.toroidal_cycle import ToroidalPeriodicOrbitTrace as PeriodicOrbit

# Build a PeriodicOrbit object that pyna needs
rzphi0 = np.array([cycle_RZ(0.0)[0], cycle_RZ(0.0)[1], 0.0])
orbit = PeriodicOrbit(
    rzphi0=rzphi0,
    period_m=m,
    trajectory=np.zeros((1, 3)),  # placeholder; evolve_DPm_along_cycle re-integrates
    DPm=np.eye(2),                # placeholder
)

# Run monodromy computation
mono = evolve_DPm_along_cycle(
    field_func=pyna_field,
    orbit=orbit,
    dt_вывод=2 * np.pi / 200,  # ~200 вывод points per turn
    rtol=1e-10, atol=1e-12,
)

print("DPm (матрица монодромии):")
print(mono.DPm)
print()
print("Собственные значения :", mono.eigenvalues)
print("Индекс устойчивости (Tr/2) :", mono.stability_index)
print("Остаток Грина         :", mono.Greene_residue)

DPm (матрица монодромии):
[[ 1.02006674 -0.07328026]
 [-0.55316615  1.02006674]]

Собственные значения : [0.81873081 1.22140267]
Индекс устойчивости (Tr/2) : 1.0200667408183737
Остаток Грина             : -0.010033370409186837

[5]:

# ── Analytic verification ────────────────────────────────────────────────────
#
# The analytic DPm at phi=0 should be V(2mπ) · diag(λ_u, λ_s) · V^{-1}(0)
# where V(φ) is the eigenvector matrix.

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

Аналитическая DPm:
[[ 1.02006676 -0.07328031]
 [-0.55316612  1.02006676]]

DPm из pyna:
[[ 1.02006674 -0.07328026]
 [-0.55316615  1.02006674]]

Макс. ошибка: 5.4677180783002655e-08

4. DX_pol вдоль орбиты и начальное условие DPm#

Здесь визуализируется DX_pol(φ) - матрица перехода состояния в каждой точке φ - и объясняется, почему уравнение DPm фазы 2 начинается с DX_pol(φ_end), а не с I.

Геометрический смысл: DPm(φ) - это монодромия “сдвинутой” орбиты, которая начинается в φ вместо φ₀. Поскольку полнооборотная карта меняется при сдвиге сечения Пуанкаре, DPm(φ) в общем случае не равно DX_pol(φ_end); оно эволюционирует по уравнению коммутатора, а его начальное значение должно равняться DX_pol(φ_end), чтобы быть согласованным с определением.

[6]:

phi_arr = mono.phi_arr
fig, axes = plt.subplots(2, 2, figsize=(12, 8), sharex=True)
fig.suptitle("Компоненты DX_pol(φ) вдоль X-цикла m=3", fontsize=14)
labels = [["DX_pol[0,0]", "DX_pol[0,1]"], ["DX_pol[1,0]", "DX_pol[1,1]"]]

for i in range(2):
    for j in range(2):
        axes[i, j].plot(phi_arr / np.pi, mono.DX_pol_arr[:, i, j], color='steelblue')
        axes[i, j].set_ylabel(labels[i][j])
        axes[i, j].axvline(0, color='k', lw=0.5)
        axes[i, j].axvline(2 * m, color='gray', lw=0.5, linestyle='--', label='φ_end')
        # Mark the initial value of DPm (= DX_pol at phi_end)
        axes[i, j].scatter([2 * m], [mono.DX_pol_arr[-1, i, j]],
                            color='red', zorder=5, label='DPm IC')

axes[0, 0].legend(fontsize=9)
for ax in axes[1]:
    ax.set_xlabel("φ / π")
plt.tight_layout()
plt.savefig("DX_pol_components.png", dpi=120)
plt.show()

print("DX_pol(φ_end) = DPm_init (НУ фазы 2):")
print(mono.DX_pol_arr[-1])
print()
print("DPm (результат фазы 1 = монодромия):")
print(mono.DPm)
print()
print("Они совпадают:", np.allclose(mono.DX_pol_arr[-1], mono.DPm))

DX_pol(φ_end) = DPm_init (НУ фазы 2):
[[ 1.02006674 -0.07328026]
 [-0.55316615  1.02006674]]

DPm (результат фазы 1 = монодромия):
[[ 1.02006674 -0.07328026]
 [-0.55316615  1.02006674]]

Они совпадают: True

5. Линейный сдвиг орбиты под действием поля возмущения#

Теперь применяется малое возмущение δB, и вычисляется линейный сдвиг X-цикла с помощью orbit_shift_under_perturbation.

Сдвиг δX(φ) удовлетворяет:

\[\frac{d\,\delta X}{d\varphi} = A(\varphi)\,\delta X + \delta b_{\rm pol}(\varphi)\]

с периодическим граничным условием δX(φ₀ + 2mπ) = δX(φ₀), что дает

\[\delta X(\varphi_0) = (\mathtt{DPm} - I)^{-1}\,(-\delta X_{\rm particular}(\varphi_{\rm end}))\]

Используется простое аналитическое возмущение: равномерный вертикальный сдвиг δB_Z = ε.

[7]:

from pyna.topo.monodromy import orbit_shift_under_perturbation

epsilon = 0.01  # perturbation amplitude

def delta_field(rzphi):
    """Perturbation: uniform δB_Z = ε (arc-length parameterised)."""
    R, Z, phi = rzphi[0], rzphi[1], rzphi[2]
    X_pol = np.array([R, Z])
    Bp = BPhi(phi, X_pol)
    # dφ/dl from unperturbed field
    g = field_g_pol(phi, X_pol)
    dphi_dl = 1.0 / np.sqrt(R**2 + np.dot(g, g))
    # δB_Z = ε ?δ(dZ/dl) = ε · |B| / |B| = ε / |B| · |B| = ε dphi_dl · R / Bp · Bp = ...
    # Simple: δdZ/dl = ε · dphi_dl / dphi_dl = ε (unit B)
    return np.array([0.0, epsilon * dphi_dl, 0.0])

delta_X = orbit_shift_under_perturbation(
    field_func=pyna_field,
    delta_field_func=delta_field,
    orbit=orbit,
    monodromy_analysis=mono,
)

print(f"Сдвиг орбиты при φ=0: δR={delta_X[0, 0]:.6f} m,  δZ={delta_X[0, 1]:.6f} m")
print(f"Проверка периодичности: |δX(end) - δX(start)| = "
      f"{np.linalg.norm(delta_X[-1] - delta_X[0]):.2e}  (должна быть < 1e-6)")

Сдвиг орбиты при φ=0: δR=-0.029622 м,  δZ=0.000000 м
Проверка периодичности: |δX(end) - δX(start)| = 1.36e-15  (должно быть < 1e-6)

[8]:

# ── Visualise the orbit shift in 3D ─────────────────────────────────────────
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')

Rc = mono.trajectory[:, 0]
Zc = mono.trajectory[:, 1]
phi_plot = mono.phi_arr

# Unperturbed cycle (xyz)
ax.plot(Rc * np.cos(phi_plot), Rc * np.sin(phi_plot), Zc,
        color='black', linewidth=2, label='Unperturbed X-cycle')

# Shifted cycle
R_shifted = Rc + delta_X[:, 0]
Z_shifted = Zc + delta_X[:, 1]
ax.plot(R_shifted * np.cos(phi_plot), R_shifted * np.sin(phi_plot), Z_shifted,
        color='tomato', linewidth=2, linestyle='--', label=f'Shifted cycle (ε={epsilon})')

ax.set_xlabel('x [m]')
ax.set_ylabel('y [m]')
ax.set_zlabel('Z [м]')
ax.legend()
ax.set_title('Линейный сдвиг орбиты под равномерным возмущением δB_Z')
plt.tight_layout()
plt.savefig("orbit_shift_3d.png", dpi=120)
plt.show()

6. Ключевые выводы#

Понятие

Символ

API pyna

Вариационная матрица вдоль орбиты

DX_pol(φ)

mono.DX_pol_arr, mono.DX_pol_at(φ)

Якобиан полнооборотной карты Пуанкаре

DPm

mono.DPm

Начальное условие фазы 2

DPm(φ₀) = DX_pol(φ_end)

автоматически в evolve_DPm_along_cycle

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

λ_u, λ_s

mono.eigenvalues

Индекс устойчивости

Tr(DPm)/2

mono.stability_index

Остаток Грина

(2 - Tr(DPm))/4

mono.Greene_residue

Линейный сдвиг орбиты

δX(φ)

orbit_shift_under_perturbation(...)
Почему ``DX_pol``, а не ``J``?
Суффикс _pol явно показывает, что матрица живет в полоидальном сечении \((R, Z)\), а не в декартовом пространстве. Имя переменной сразу несет физический смысл.
Почему ``DPm``, а не ``M``?
DPm = производная (D) карты Пуанкаре (P) за m оборотов. Одна буква M не раскрывает эту структуру.