Анализ якобиана магнитных островов: периодические орбиты и матрицы монодромии#
Этот ноутбук демонстрирует полный анализ якобиана и монодромии для цепочек магнитных островов в аналитическом равновесии стелларатора, перенесенный из рабочего процесса в стиле Julia MCF_scripts.
Научный контекст#
Периодическая орбита силовой линии (цикл) - это силовая линия, которая возвращается в начальную точку ровно после \(n\) тороидальных оборотов: \(X(\phi + 2\pi n) = X(\phi)\).
Это неподвижные точки \(n\)-й итерации карты Пуанкаре \(P^n\). Для резонанса \(q = m/n\) (где \(q = B_\phi r / (B_{pol} R)\)) период орбиты равен \(m\) оборотам.
Матрица Якоби \(DX(\phi)\) эволюционирует как:
\[\frac{dDX}{d\phi} = A(r(\phi), z(\phi), \phi) \cdot DX, \quad DX(\phi_0) = I\]
где \(A_{ij} = \partial(R B_{pol,i}/B_\phi) / \partial x_j\).
Матрица монодромии \(M = DP(\phi_0 + 2\pi n)\) задает линеаризованную \(n\)-оборотную карту. Собственные значения \(M\):
\(|\lambda| = 1\): эллиптическая точка (O-точка, устойчивая)
\(|\lambda| > 1\): гиперболическая точка (X-точка, неустойчивая)
[1]:
import numpy as np
import matplotlib.pyplot as plt
from pyna.toroidal.equilibrium.stellarator import StellaratorSimple
from pyna.topo.toroidal_cycle import (
poincare_map_n, poincare_map_n_trajectory,
jacobian_of_poincare_map, find_cycle, find_all_cycles_near_resonance,
ToroidalPeriodicOrbitTrace as PeriodicOrbit,
)
from pyna.topo.monodromy import (
evolve_DPm_along_cycle, build_A_matrix_func, build_delta_b_pol_func,
orbit_shift_under_perturbation, monodromy_change_under_perturbation,
CycleVariationalData,
)
plt.rcParams['figure.figsize'] = (10, 7)
plt.rcParams['font.size'] = 12
print('Импорты OK')
Импорты OK
1. Настройка: построить StellaratorSimple с цепочкой островов \(q=2/1\)#
Используется m_h=2, n_h=1, epsilon_h=0.02, чтобы создать резонанс \(q=2\). В этой модели q = m_h/n_h = 2, а цепочка островов имеет период 2 (2 тороидальных оборота).
q-профиль: q(r) = q0 + (q1-q0)*psi = q0 + (q1-q0)*(r/r0)²
При \(r=0\): \(q = q0 = 1.5\)
При \(r=r0\): \(q = q1 = 3.5\)
Поверхность \(q=2\):
psi = (q-q0)/(q1-q0) = 0.5/2.0 = 0.25, значит \(r/r0 = 0.5\)
[2]:
# Build the stellarator
epsilon_h = 0.02 # helical ripple amplitude
stellarator = StellaratorSimple(
R0=3.0, r0=0.35, B0=1.0,
q0=1.5, q1=3.5,
m_h=2, n_h=1, epsilon_h=epsilon_h,
)
R0, r0 = stellarator.R0, stellarator.r0
field_func = stellarator.field_func
# Domain limits
RZlimit = (R0 - r0*1.5, R0 + r0*1.5, -r0*1.5, r0*1.5)
# Find q=2/1 resonant surface
psi_list = stellarator.resonant_psi(2, 1)
psi_res = psi_list[0]
r_res = np.sqrt(psi_res) * r0
print(f"Резонанс q=2/1 при psi={psi_res:.3f}, r_res={r_res:.4f} m")
# q profile plot
psi_arr = np.linspace(0, 1, 100)
q_arr = [stellarator.q_of_psi(p) for p in psi_arr]
fig, ax = plt.subplots(1, 1, figsize=(8, 4))
ax.plot(psi_arr, q_arr, 'b-', linewidth=2)
ax.axhline(2.0, color='r', linestyle='--', label='резонанс q=2')
ax.axvline(psi_res, color='r', linestyle=':', alpha=0.5)
ax.set_xlabel('Нормированный поток ψ')
ax.set_ylabel('Коэффициент безопасности q')
ax.set_title('q-профиль StellaratorSimple')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('q_profile.png', dpi=100, bbox_inches='tight')
plt.show()
print('График q-профиля сохранен')
Резонанс q=2/1 при psi=0.250, r_res=0.1750 м
График q-профиля сохранен
2. Сечение Пуанкаре и визуализация цепочки островов#
[3]:
# Compute Poincaré section
print('Вычисление сечения Пуанкаре (это может занять ~30 с)...')
# Trace field lines from various starting точек
n_field_lines = 8
n_poincare_turns = 100
poincare_точек = [] # (R, Z) lists per field line
for i, r_frac in enumerate(np.linspace(0.15, 0.95, n_field_lines)):
r_start = r_frac * r0
R_start = R0 + r_start
Z_start = 0.0
Rs, Zs = [R_start], [Z_start]
R_curr, Z_curr = R_start, Z_start
for _ in range(n_poincare_turns):
R_next, Z_next = poincare_map_n(
field_func, [R_curr, Z_curr, 0.0], n_turns=1, dt=0.15, RZlimit=RZlimit
)
if np.isnan(R_next):
break
Rs.append(R_next)
Zs.append(Z_next)
R_curr, Z_curr = R_next, Z_next
poincare_точек.append((np.array(Rs), np.array(Zs)))
print(f' Силовая линия {i+1}/{n_field_lines}: {len(Rs)} точек')
print('Сечение Пуанкаре вычислено')
Вычисление сечения Пуанкаре (это может занять ~30 с)...
Силовая линия 1/8: 101 точка
Силовая линия 2/8: 101 точка
Силовая линия 3/8: 101 точка
Силовая линия 4/8: 101 точка
Силовая линия 5/8: 101 точка
Силовая линия 6/8: 101 точка
Силовая линия 7/8: 101 точка
Силовая линия 8/8: 101 точка
Сечение Пуанкаре вычислено
[4]:
# Find O-точек and X-точек using Newton-Raphson
print('Finding X/O точек...')
# Known approximate locations from scan (for epsilon_h=0.02, q=2/1)
# X-точек: near theta=±0.79 and ±2.36 at r=r_res
# O-точек: near theta=±2.32 at r=r_res (for small epsilon)
cycle_seeds = [
# X-point candidates
np.array([R0 + r_res*np.cos(0.79), r_res*np.sin(0.79), 0.0]),
np.array([R0 + r_res*np.cos(-0.79), r_res*np.sin(-0.79), 0.0]),
np.array([R0 + r_res*np.cos(2.36), r_res*np.sin(2.36), 0.0]),
np.array([R0 + r_res*np.cos(-2.36), r_res*np.sin(-2.36), 0.0]),
# O-point candidates (scan found distance minimum near theta=-2.02, r=0.1625)
np.array([R0 + r_res*np.cos(3.02), r_res*np.sin(3.02), 0.0]),
np.array([R0 + r_res*np.cos(-2.32), r_res*np.sin(-2.32), 0.0]),
np.array([R0 + r_res*np.cos(2.32), r_res*np.sin(2.32), 0.0]),
np.array([R0 + r_res*np.cos(0.25), r_res*np.sin(0.25), 0.0]),
]
found_orbits = []
for seed in cycle_seeds:
orbit = find_cycle(
field_func, seed, n_turns=2, dt=0.15, RZlimit=RZlimit,
max_iter=30, tol=1e-8, n_fallback_seeds=4, fallback_radius=0.02,
)
if orbit is None:
continue
dist_axis = np.sqrt((orbit.rzphi0[0]-R0)**2 + orbit.rzphi0[1]**2)
if dist_axis < 0.05:
continue # skip axis fixed point
# Deduplicate
dup = any(np.linalg.norm(orbit.rzphi0[:2]-fo.rzphi0[:2]) < 1e-4 for fo in found_orbits)
if not dup:
found_orbits.append(orbit)
o_точек = [o for o in found_orbits if o.is_stable]
x_точек = [o for o in found_orbits if not o.is_stable]
print(f'Найдено {len(o_точек)} O-точек and {len(x_точек)} X-точек')
for o in o_точек:
print(f' O-точка: ({o.rzphi0[0]:.4f}, {o.rzphi0[1]:.4f}), k={o.stability_index:.4f}')
for x in x_точек:
print(f' X-точка: ({x.rzphi0[0]:.4f}, {x.rzphi0[1]:.4f}), k={x.stability_index:.4f}')
Поиск X/O-точек...
Найдено 2 O-точек и 4 X-точек
O-точка: (3.1237, -0.1237), k=-0.1045
O-точка: (2.8763, -0.1237), k=0.3026
X-точка: (3.1237, 0.1237), k=0.3024
X-точка: (2.8763, 0.1237), k=-0.1042
X-точка: (2.8383, 0.0169), k=2.0324
X-точка: (3.1587, 0.0546), k=2.0320
[5]:
# Plot Poincaré section with X/O точек marked
fig, ax = plt.subplots(figsize=(10, 8))
colors = plt.cm.tab10(np.linspace(0, 0.9, n_field_lines))
for i, (Rs, Zs) in enumerate(poincare_точек):
ax.scatter(Rs, Zs, s=0.5, color=colors[i], alpha=0.6)
# Mark O-точек (blue circles)
for o in o_точек:
ax.plot(o.rzphi0[0], o.rzphi0[1], 'bo', markersize=10,
label='O-point' if o is o_точек[0] else '', zorder=5)
# Mark X-точек (red x's)
for x in x_точек:
ax.plot(x.rzphi0[0], x.rzphi0[1], 'rx', markersize=12, markeredgewidth=2,
label='X-point' if x is x_точек[0] else '', zorder=5)
# Mark resonant surface
theta_arr = np.linspace(0, 2*np.pi, 200)
ax.plot(R0 + r_res*np.cos(theta_arr), r_res*np.sin(theta_arr), 'k--',
linewidth=1, alpha=0.5, label=f'q=2 surface (r={r_res:.3f})')
ax.set_xlabel('R (м)')
ax.set_ylabel('Z (м)')
ax.set_title(f'Сечение Пуанкаре - цепочка островов q=2/1\n(epsilon_h={epsilon_h})')
ax.set_aspect('equal')
ax.legend(loc='upper right')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('poincare_section.png', dpi=100, bbox_inches='tight')
plt.show()
print('Сечение Пуанкаре сохранено')
# Note: For interactive Plotly version, use:
# import plotly.graph_objects as go
# fig = go.Figure()
# for Rs, Zs in poincare_точек:
# fig.add_trace(go.Scatter(x=Rs, y=Zs, mode='markers', marker=dict(size=1)))
# for o in o_точек:
# fig.add_trace(go.Scatter(x=[o.rzphi0[0]], y=[o.rzphi0[1]], mode='markers',
# marker=dict(size=10, symbol='circle', color='blue')))
# fig.show()
Сечение Пуанкаре сохранено
3. Анализ монодромии для O-точки#
Матрица монодромии \(M = DP(\phi_0 + 2\pi n)\) кодирует устойчивость периодической орбиты. Для O-точки (эллиптической) собственные значения лежат на единичной окружности: \(|\lambda| = 1\).
Остаток Грина: \(R_G = (2 - \text{Tr}(M))/4\)
\(0 < R_G < 1\): эллиптическая (устойчивая)
\(R_G < 0\): гиперболическая (неустойчивая, стандартная)
\(R_G > 1\): гиперболическая с отражением
[6]:
# Pick the first O-point for analysis
if not o_точек:
print('No O-точек found; using first found orbit for demo')
opoint = found_orbits[0] if found_orbits else None
else:
opoint = o_точек[0]
if opoint is None:
print('No orbits found to analyze')
else:
print(f'Анализ O-точки в ({opoint.rzphi0[0]:.4f}, {opoint.rzphi0[1]:.4f})')
print(f' stability_index = {opoint.stability_index:.6f}')
print(f' собственные значения = {opoint.eigenvalues}')
# Compute full monodromy analysis
print('Вычисление монодромии (вариационные уравнения)...')
monodromy_O = evolve_DPm_along_cycle(
field_func, opoint, dt_вывод=0.1, rtol=1e-8, atol=1e-9
)
print(f' Матрица монодромии M:')
print(f' {monodromy_O.DPm}')
print(f' det(M) = {np.linalg.det(monodromy_O.DPm):.8f} (должно быть ~1)')
print(f' Tr(M)/2 = {monodromy_O.stability_index:.6f}')
print(f" Остаток Грина = {monodromy_O.Greene_residue:.6f}")
print(f' Собственные значения = {monodromy_O.eigenvalues}')
Анализ O-точки в (3.1237, -0.1237)
stability_index = -0.104472
собственные значения = [-0.10447198+0.99401547j -0.10447198-0.99401547j]
Вычисление монодромии (вариационные уравнения)...
Матрица монодромии M:
[[-1.95431405 1.21064084]
[-3.64518107 1.74639533]]
det(M) = 1.00000011 (должно быть ~1)
Tr(M)/2 = -0.103959
Остаток Грина = 0.551980
Собственные значения = [-0.10395936+0.9945816j -0.10395936-0.9945816j]
[7]:
if opoint is not None and 'monodromy_O' in dir():
phi_arr = monodromy_O.phi_arr
# Compute eigenvalues and det along orbit
DP_eigvals = np.array([np.linalg.eigvals(DP) for DP in monodromy_O.DPm_arr])
DP_dets = np.array([np.linalg.det(DP) for DP in monodromy_O.DPm_arr])
fig, axes = plt.subplots(3, 1, figsize=(10, 10))
# Plot |eigenvalues| of DPm along orbit
ax = axes[0]
ax.plot(phi_arr/(2*np.pi), np.abs(DP_eigvals[:, 0]), 'b-', label='|lambda(phi)|')
ax.plot(phi_arr/(2*np.pi), np.abs(DP_eigvals[:, 1]), 'r--', label='|lambda(phi)|')
ax.axhline(1.0, color='k', linestyle=':', alpha=0.5)
ax.set_xlabel('φ/(2π)')
ax.set_ylabel('|eigenvalue of DP|')
ax.set_title('O-точка: Эволюция собственных значений вдоль орбиты')
ax.legend()
ax.grid(True, alpha=0.3)
# Plot det(DPm) - should stay near 1
ax = axes[1]
ax.plot(phi_arr/(2*np.pi), DP_dets, 'g-', linewidth=2)
ax.axhline(1.0, color='k', linestyle=':', alpha=0.5, label='det=1 (area-preserving)')
ax.set_xlabel('φ/(2π)')
ax.set_ylabel('det(DPm(φ))')
ax.set_title('Проверка сохранения площади: det(DPm) ~ 1')
ax.legend()
ax.grid(True, alpha=0.3)
# Plot Greene residue evolution (Tr(J)/2)
ax = axes[2]
Tr_half = np.array([np.trace(J)/2 for J in monodromy_O.DPm_arr])
ax.plot(phi_arr/(2*np.pi), Tr_half, 'm-', linewidth=2)
ax.axhline(0.0, color='k', linestyle=':', alpha=0.3)
ax.axhline(monodromy_O.stability_index, color='r', linestyle='--',
label=f'Final Tr(M)/2 = {monodromy_O.stability_index:.4f}')
ax.set_xlabel('φ/(2π)')
ax.set_ylabel('Tr(DPm(φ))/2')
ax.set_title('Эволюция индекса устойчивости вдоль орбиты O-точки')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('monodromy_opoint.png', dpi=100, bbox_inches='tight')
plt.show()
print('График монодромии O-точки сохранен')
График монодромии O-точки сохранен
4. Анализ монодромии для X-точки#
Для X-точки (гиперболической) собственные значения вещественны и \(|\lambda| > 1\). Соседние орбиты экспоненциально расходятся; обратный показатель Ляпунова задается \(\lambda_{max}\).
[8]:
# Pick the first X-point
if not x_точек:
print('No X-точек found')
xpoint = None
else:
xpoint = x_точек[0]
print(f'Анализ X-точки в ({xpoint.rzphi0[0]:.4f}, {xpoint.rzphi0[1]:.4f})')
print(f' stability_index = {xpoint.stability_index:.6f}')
print(f' собственные значения = {xpoint.eigenvalues}')
monodromy_X = evolve_DPm_along_cycle(
field_func, xpoint, dt_вывод=0.1, rtol=1e-8, atol=1e-9
)
print(f' det(M) = {np.linalg.det(monodromy_X.DPm):.8f}')
print(f" Остаток Грина = {monodromy_X.Greene_residue:.6f} (< 0 -> гиперболическая)")
Анализ X-точки в (3.1237, 0.1237)
stability_index = 0.302381
собственные значения = [0.30238099+0.95326451j 0.30238099-0.95326451j]
det(M) = 0.99999969
Остаток Грина = 0.348880 (< 0 -> гиперболическая)
[9]:
if xpoint is not None and 'monodromy_X' in dir():
phi_arr_X = monodromy_X.phi_arr
DP_eigvals_X = np.array([np.linalg.eigvals(DP) for DP in monodromy_X.DPm_arr])
DP_dets_X = np.array([np.linalg.det(DP) for DP in monodromy_X.DPm_arr])
fig, axes = plt.subplots(2, 1, figsize=(10, 7))
ax = axes[0]
ax.semilogy(phi_arr_X/(2*np.pi), np.abs(DP_eigvals_X[:, 0]), 'b-', label='|lambda(phi)|')
ax.semilogy(phi_arr_X/(2*np.pi), np.abs(DP_eigvals_X[:, 1]), 'r--', label='|lambda(phi)|')
ax.axhline(1.0, color='k', linestyle=':', alpha=0.5)
ax.set_xlabel('φ/(2π)')
ax.set_ylabel('|eigenvalue of DP| (log scale)')
ax.set_title('X-точка: Рост гиперболического собственного значения вдоль орбиты')
ax.legend()
ax.grid(True, alpha=0.3)
ax = axes[1]
ax.plot(phi_arr_X/(2*np.pi), DP_dets_X, 'g-', linewidth=2)
ax.axhline(1.0, color='k', linestyle=':', alpha=0.5, label='det=1')
ax.set_xlabel('φ/(2π)')
ax.set_ylabel('det(DPm(φ))')
ax.set_title('Проверка сохранения площади for X-point orbit')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('monodromy_xpoint.png', dpi=100, bbox_inches='tight')
plt.show()
5b. Устойчивые и неустойчивые многообразия из X-точки#
Устойчивое многообразие \(W^s\) (синее, цветовая карта cool) содержит траектории, которые сходятся к X-точке при прямой итерации. Неустойчивое многообразие \(W^u\) (оранжевое, цветовая карта warm) содержит траектории, которые расходятся от X-точки. Вместе они образуют сепаратрисный скелет цепочки островов.
[10]:
# === Stable/Unstable Manifold Visualization ===
from pyna.topo.variational import PoincareMapVariationalEquations
from pyna.topo.manifold_improve import StableManifold, UnstableManifold
from pyna.toroidal.visual.tokamak_manifold import _manifold_line_collection, manifold_legend_handles
# field_func_2d wrapper: stellarator.field_func([R,Z,phi]) returns unit tangent
# We need (R,Z,phi) -> (dR/dphi, dZ/dphi)
def field_func_2d(R, Z, phi):
tang = stellarator.field_func(np.array([R, Z, phi])) # [dR/ds, dZ/ds, dphi/ds]
dphi_ds = tang[2]
if abs(dphi_ds) < 1e-15:
return np.array([0.0, 0.0])
return np.array([tang[0] / dphi_ds, tang[1] / dphi_ds])
if xpoint is not None and 'monodromy_X' in dir():
R_xpt, Z_xpt = xpoint.rzphi0[0], xpoint.rzphi0[1]
print(f'Рост многообразий из X-точки ({R_xpt:.4f}, {Z_xpt:.4f})')
# Use the Jacobian from monodromy_X (2x2 monodromy matrix)
M = monodromy_X.DPm
eigvals = np.linalg.eigvals(M)
print(f'Собственные значения монодромии: {eigvals}')
print(f'|lambda_1|={abs(eigvals[0]):.4f}, |lambda_2|={abs(eigvals[1]):.4f} (произведение~1: {abs(np.prod(eigvals)):.4f})')
# Grow manifolds
sm = StableManifold([R_xpt, Z_xpt], M, field_func_2d)
um = UnstableManifold([R_xpt, Z_xpt], M, field_func_2d)
sm.grow(n_turns=12, init_length=5e-5, n_init_pts=6, both_sides=True)
um.grow(n_turns=12, init_length=5e-5, n_init_pts=6, both_sides=True)
print(f'Stable segments: {len(sm.segments)}, Unstable segments: {len(um.segments)}')
# Combined figure: Poincare + manifolds
fig_mf, ax_mf = plt.subplots(figsize=(9, 8))
ax_mf.set_facecolor('white')
# Replot Poincare scatter (from existing results_n variable)
try:
R_pts = results_n[:, 0]; Z_pts = results_n[:, 1]
psi_n = np.clip(((R_pts - R0)**2 + Z_pts**2) / r_res**2, 0, 1.4)
from matplotlib.colors import Normalize
import matplotlib.cm as cm
colors_sc = cm.plasma(np.clip(psi_n / 1.4, 0.05, 0.95))
ax_mf.scatter(R_pts, Z_pts, s=0.5, c=colors_sc, rasterized=True, alpha=0.5, zorder=2)
except NameError:
pass # Poincare results not available, skip
# Manifolds with arc-length coloring
for seg in sm.segments:
if len(seg) > 2:
lc, _ = _manifold_line_collection(seg, cmap='GnBu')
lc.set_linewidth(1.3); lc.set_alpha(0.92); lc.set_zorder(6)
ax_mf.add_collection(lc)
for seg in um.segments:
if len(seg) > 2:
lc, _ = _manifold_line_collection(seg, cmap='Oranges')
lc.set_linewidth(1.3); lc.set_alpha(0.92); lc.set_zorder(6)
ax_mf.add_collection(lc)
# X-point marker
ax_mf.plot(R_xpt, Z_xpt, 'r+', ms=14, mew=2.5, zorder=8, label='X-point')
# Also plot all found O/X точек
for op in o_точек:
ax_mf.plot(op.rzphi0[0], op.rzphi0[1], 'go', ms=7, mew=1.5, zorder=7)
for xp in x_точек:
ax_mf.plot(xp.rzphi0[0], xp.rzphi0[1], 'r+', ms=12, mew=2.5, zorder=8)
# Resonant surface circle
theta_c = np.linspace(0, 2*np.pi, 300)
ax_mf.plot(R0 + r_res*np.cos(theta_c), r_res*np.sin(theta_c),
'--', color='tomato', lw=0.8, alpha=0.6, label='$q=2/1$ surface')
ax_mf.plot(R0 + stellarator.r0*np.cos(theta_c), stellarator.r0*np.sin(theta_c),
'k-', lw=1.2, label='LCFS')
# Legend + labels
mfld_handles = manifold_legend_handles('Oranges', 'GnBu')
ax_mf.legend(handles=mfld_handles + [
plt.Line2D([0],[0], marker='+', color='r', ms=10, mew=2, lw=0, label='X-point'),
plt.Line2D([0],[0], marker='o', color='g', ms=7, lw=0, label='O-point'),
plt.Line2D([0],[0], color='k', lw=1.2, label='LCFS'),
], loc='upper right', fontsize=9, framealpha=0.9)
ax_mf.set_xlim(R0 - 1.3*stellarator.r0, R0 + 1.3*stellarator.r0)
ax_mf.set_ylim(-1.3*stellarator.r0, 1.3*stellarator.r0)
ax_mf.set_xlabel('$R$ (m)', fontsize=12)
ax_mf.set_ylabel('$Z$ (m)', fontsize=12)
ax_mf.set_title('Устойчивые ($W^s$, синие) и неустойчивые ($W^u$, оранжевые) многообразия - X-точка $q=2/1$', fontsize=12)
ax_mf.set_aspect('equal')
plt.tight_layout()
plt.savefig('island_manifolds.png', dpi=150, bbox_inches='tight')
plt.show()
print('Фигура многообразий готова.')
else:
print('Нет X-точки или монодромии; визуализация многообразий пропущена.')
Рост многообразий из X-точки (3.1237, 0.1237)
Собственные значения монодромии: [0.3022398+0.95323176j 0.3022398-0.95323176j]
|lambda_1|=1.0000, |lambda_2|=1.0000 (произведение~1: 1.0000)
Устойчивых сегментов: 2, неустойчивых сегментов: 2
Фигура многообразий готова.
5. Влияние возмущения: сдвиг орбиты под действием \(\delta B\)#
Добавляется небольшое дополнительное винтовое возмущение \(\delta B\) с амплитудой epsilon_pert. Смещение орбиты удовлетворяет:
\[\frac{dX_{cyc}}{d\phi} = A \cdot X_{cyc} + \delta b_{pol}\]
где \(\delta b_{pol} = [R\delta B_R/B_\phi - R B_R \delta B_\phi/B_\phi^2; \text{то же для } Z]\).
[11]:
epsilon_pert = 0.005 # амплитуда возмущения
# Build perturbed stellarator (slightly different helical phase to create delta B)
stellarator_pert = StellaratorSimple(
R0=3.0, r0=0.35, B0=1.0,
q0=1.5, q1=3.5,
m_h=2, n_h=1, epsilon_h=epsilon_h + epsilon_pert,
)
def delta_field_func(rzphi):
"""Perturbation field: difference between perturbed and unperturbed."""
f0 = np.asarray(field_func(rzphi))
f1 = np.asarray(stellarator_pert.field_func(rzphi))
# Note: field_func returns unit tangent vectors, not raw B fields.
# For the orbit displacement, we use the difference in the phi-parameterized field.
# The perturbation in terms of dR/dphi, dZ/dphi:
g0 = np.array([f0[0]/f0[2], f0[1]/f0[2], 1.0])
g1 = np.array([f1[0]/f1[2], f1[1]/f1[2], 1.0])
return g1 - g0
if opoint is not None and 'monodromy_O' in dir():
print(f'Вычисление сдвига орбиты для O-point при epsilon_pert={epsilon_pert}...')
orbit_shift_O = orbit_shift_under_perturbation(
field_func, delta_field_func, opoint, monodromy_O
)
print(f' Начальное смещение орбиты: dR={orbit_shift_O[0,0]:.6f}, dZ={orbit_shift_O[0,1]:.6f}')
print(f' Макс. |смещение| вдоль орбиты: {np.max(np.linalg.norm(orbit_shift_O, axis=1)):.6f} m')
if xpoint is not None and 'monodromy_X' in dir():
print(f'Вычисление сдвига орбиты для X-point при epsilon_pert={epsilon_pert}...')
orbit_shift_X = orbit_shift_under_perturbation(
field_func, delta_field_func, xpoint, monodromy_X
)
print(f' Начальное смещение орбиты: dR={orbit_shift_X[0,0]:.6f}, dZ={orbit_shift_X[0,1]:.6f}')
Вычисление сдвига орбиты O-точки при epsilon_pert=0.005...
Начальное смещение орбиты: dR=-0.000000, dZ=-0.000000
Макс. |смещение| вдоль орбиты: 0.000000 м
Вычисление сдвига орбиты X-точки при epsilon_pert=0.005...
Начальное смещение орбиты: dR=0.000000, dZ=-0.000000
[12]:
if opoint is not None and 'orbit_shift_O' in dir():
phi_arr_O = monodromy_O.phi_arr
fig, axes = plt.subplots(2, 1, figsize=(10, 7))
ax = axes[0]
ax.plot(phi_arr_O/(2*np.pi), orbit_shift_O[:, 0], 'b-', label='δR(φ)')
ax.plot(phi_arr_O/(2*np.pi), orbit_shift_O[:, 1], 'r-', label='δZ(φ)')
ax.set_xlabel('φ/(2π)')
ax.set_ylabel('Смещение орбиты (м)')
ax.set_title(f'O-point сдвиг орбиты при δε={epsilon_pert}')
ax.legend()
ax.grid(True, alpha=0.3)
if xpoint is not None and 'orbit_shift_X' in dir():
phi_arr_X = monodromy_X.phi_arr
ax = axes[1]
ax.plot(phi_arr_X/(2*np.pi), orbit_shift_X[:, 0], 'b-', label='δR(φ)')
ax.plot(phi_arr_X/(2*np.pi), orbit_shift_X[:, 1], 'r-', label='δZ(φ)')
ax.set_xlabel('φ/(2π)')
ax.set_ylabel('Смещение орбиты (м)')
ax.set_title(f'X-point сдвиг орбиты при δε={epsilon_pert}')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('orbit_shift.png', dpi=100, bbox_inches='tight')
plt.show()
6. Диаграмма устойчивости: остаток Грина против амплитуды возмущения#
Сканируется epsilon_pert от 0 до 0.02. Для каждой амплитуды:
Найти O- и X-циклы
Вычислить остаток Грина \(R_G = (2 - \text{Tr}(M))/4\)
Переход \(R_G: 0 \to 1\) отмечает разрушение острова (KAM-подобный переход).
[13]:
print('Вычисление диаграммы устойчивости (сканирование epsilon_h)...')
epsilon_arr = np.linspace(0.005, 0.04, 12)
greene_O = [] # O-point Greene residues
greene_X = [] # X-point Greene residues
valid_eps = []
for eps in epsilon_arr:
st_eps = StellaratorSimple(
R0=3.0, r0=0.35, B0=1.0,
q0=1.5, q1=3.5,
m_h=2, n_h=1, epsilon_h=eps,
)
ff = st_eps.field_func
# Find X-point near known seed (2.8763, -0.1237)
orbit_x = find_cycle(
ff, np.array([2.8763, -0.1237, 0.0]),
n_turns=2, dt=0.15, RZlimit=RZlimit,
max_iter=30, tol=1e-8, n_fallback_seeds=6, fallback_radius=0.02,
)
# Find O-point near known seed
psi_r = st_eps.resonant_psi(2, 1)[0]
r_r = np.sqrt(psi_r) * st_eps.r0
opoint_seed = np.array([st_eps.R0 + r_r*np.cos(-2.32), r_r*np.sin(-2.32), 0.0])
orbit_o = find_cycle(
ff, opoint_seed,
n_turns=2, dt=0.15, RZlimit=RZlimit,
max_iter=30, tol=1e-8, n_fallback_seeds=8, fallback_radius=0.02,
)
rg_x = None
rg_o = None
if orbit_x is not None:
M_x = orbit_x.DPm
rg_x = (2.0 - np.trace(M_x)) / 4.0
if orbit_o is not None:
dist_axis = np.sqrt((orbit_o.rzphi0[0]-st_eps.R0)**2 + orbit_o.rzphi0[1]**2)
if dist_axis > 0.05:
M_o = orbit_o.DPm
rg_o = (2.0 - np.trace(M_o)) / 4.0
if rg_x is not None or rg_o is not None:
valid_eps.append(eps)
greene_X.append(rg_x)
greene_O.append(rg_o)
rg_x_str = f'{rg_x:.4f}' if rg_x is not None else 'N/A'
rg_o_str = f'{rg_o:.4f}' if rg_o is not None else 'N/A'
print(f' eps={eps:.3f}: R_G(X)={rg_x_str}, R_G(O)={rg_o_str}')
print(f'Вычислено {len(valid_eps)} data точек')
Вычисление диаграммы устойчивости (сканирование epsilon_h)...
eps=0.005: R_G(X)=0.0024, R_G(O)=0.0024
eps=0.008: R_G(X)=0.0325, R_G(O)=0.0325
eps=0.011: R_G(X)=0.0852, R_G(O)=0.0852
eps=0.015: R_G(X)=0.1612, R_G(O)=0.1612
eps=0.018: R_G(X)=0.2615, R_G(O)=0.2615
eps=0.021: R_G(X)=0.3873, R_G(O)=0.3873
eps=0.024: R_G(X)=0.5402, R_G(O)=0.5402
eps=0.027: R_G(X)=0.7220, R_G(O)=0.7220
eps=0.030: R_G(X)=0.9348, R_G(O)=0.9348
eps=0.034: R_G(X)=1.1812, R_G(O)=1.1812
eps=0.037: R_G(X)=1.4640, R_G(O)=1.4640
eps=0.040: R_G(X)=1.7865, R_G(O)=1.7865
Вычислено 12 точек данных
[14]:
fig, ax = plt.subplots(figsize=(9, 5))
# Filter valid entries
eps_X = [e for e, g in zip(valid_eps, greene_X) if g is not None]
rg_X = [g for g in greene_X if g is not None]
eps_O = [e for e, g in zip(valid_eps, greene_O) if g is not None]
rg_O = [g for g in greene_O if g is not None]
if eps_X:
ax.plot(eps_X, rg_X, 'rx-', linewidth=2, markersize=8, label="X-point Остаток Грина")
if eps_O:
ax.plot(eps_O, rg_O, 'bo-', linewidth=2, markersize=8, label="O-point Остаток Грина")
ax.axhline(0.0, color='k', linestyle='--', alpha=0.5, label='R_G=0 (elliptic boundary)')
ax.axhline(1.0, color='k', linestyle=':', alpha=0.5, label='R_G=1 (гиперболическая boundary)')
ax.fill_between([0, 0.05], 0, 1, alpha=0.1, color='green', label='Elliptic region')
ax.set_xlabel('Helical ripple amplitude ε_h')
ax.set_ylabel("Остаток Грина R_G")
ax.set_title("Диаграмма устойчивости: Остаток Грина vs амплитуда возмущения\n"
f"(q=2/1 island chain, m_h=2, n_h=1)")
ax.legend(loc='best')
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 0.05)
plt.tight_layout()
plt.savefig('stability_diagram.png', dpi=100, bbox_inches='tight')
plt.show()
print('Диаграмма устойчивости saved')
Диаграмма устойчивости сохранена
Итоги#
Этот ноутбук показал:
Поиск периодических орбит методом Ньютона-Рафсона на карте Пуанкаре
Эволюцию якобиана вдоль орбиты через вариационные уравнения
Матрицу монодромии \(M = J(\phi_0 + 2\pi n)\) с det(M) ~ 1 (сохранение площади)
Остаток Грина \(R_G\) как индикатор устойчивости
Сдвиг орбиты под возмущением с использованием периодического граничного условия
Диаграмму устойчивости: как устойчивость цепочки островов меняется с амплитудой возмущения
Ключевые результаты:
O-точки имеют \(0 < R_G < 1\) (эллиптические, устойчивые)
X-точки имеют \(R_G < 0\) или \(R_G > 1\) (гиперболические, неустойчивые)
det(DPm(phi)) ~ 1на всем интервале, подтверждая сохранение площади (гамильтонова структура)
Перенос этого анализа из Julia в Python на аналитическую модель стелларатора демонстрирует те же физические алгоритмы без проприетарных данных равновесия.