Метрическая модель динамики твёрдого тела: параметрический резонанс, как проверяемое следствие

МЕТРИЧЕСКАЯ МОДЕЛЬ ДИНАМИКИ ТВЁРДОГО ТЕЛА: ПАРАМЕТРИЧЕСКИЙ РЕЗОНАНС КАК ПРОВЕРЯЕМОЕ СЛЕДСТВИЕ

А. Алексеенко
ORCID 0009-0004-4408-0097
Цикл «Рецепт хрупкого бессмертия»

Аннотация

Предложена минимальная метрическая модель динамики свободного твёрдого тела, в которой моменты инерции зависят от локального скалярного поля α: I_i(α) = I_i⁰(1 + ξ_i α). Уравнения движения выведены в формализме Эйлера–Пуанкаре из вариационного принципа с лагранжианом L = ½ Σ I_i(α) ω_i² + ½ α̇² − ½ m² α². Показано, что модель внутренне согласована (энергия сохраняется с точностью 3.8×10⁻¹⁴ в численной симуляции) и предсказывает параметрический резонанс при частоте модуляции m ≈ 2Ω_prec, где Ω_prec — частота прецессии. Численная симуляция подтверждает: при ξ = 0.1 и m = 2Ω_prec возникает резонанс с отношением частот α:ω₁ = 2.00, отсутствующий в классическом пределе ξ = 0. При физичных параметрах (α ~ 6×10⁻³, I₁ > 0) эффект умеренный и не требует нефизичных условий. Параметры ξ_i и m феноменологические; их значения требуют экспериментального определения. Работа не доказывает реальность поля α, но показывает, что метрическая рамка способна давать нетривиальные предсказания, согласующиеся с физикой.

Ключевые слова: динамика твёрдого тела, уравнения Эйлера–Пуанкаре, параметрический резонанс, уравнение Матьё, скалярное поле, метрическая теория.

1. Введение

Эффект Джанибекова — переворот свободно вращающегося тела вокруг промежуточной оси инерции — полностью объясняется классической механикой: теоремой о промежуточной оси и диссипацией (пример «Эксплореров», где гибкие антенны вносили внутреннее трение). Это объяснение полное и достаточное для механизма. Однако классика не даёт замкнутой аналитической формулы для периода переворота T(ω): динамика моделируется численно, но аналитического выражения нет из-за нелинейности и зависимости от начальных условий.

В цикле «Рецепт хрупкого бессмертия» предложена метрическая рамка, в которой постоянная тонкой структуры α рассматривается как локальное поле, а не как константа. В настоящей работе мы проверяем, может ли эта рамка дать проверяемое количественное следствие для вращающегося тела. Мы не заменяем классическое объяснение механизма, а добавляем гипотезу о параметре, которая формулируется как модель с дифференциальными уравнениями и проверяется численно.

2. Модель

2.1. Лагранжиан

Рассматривается свободное твёрдое тело с угловыми скоростями ω_i и динамическими моментами инерции

I_i(α) = I_i⁰(1 + ξ_i α), (1)

где α — безразмерное скалярное поле, ξ_i — безразмерные коэффициенты связи. Лагранжиан системы:

L = ½ Σ_{i=1}^{3} I_i(α) ω_i² + ½ α̇² − ½ m² α². (2)

Связь поля с вращением инкапсулирована в I_i(α); отдельный член взаимодействия не вводится, чтобы избежать двойного счёта.

2.2. Уравнения движения

Поскольку ω_i — квазикоординаты, применяется редукция Эйлера–Пуанкаре:

d/dt(∂L/∂ω_i) + (ω × ∂L/∂ω)_i = 0. (3)

Подстановка (1)–(2) даёт:

I_i(α) ω̇_i + ξ_i I_i⁰ α̇ ω_i + (ω × I(α) ω)_i = 0, (4)

где I(α) — диагональный тензор с компонентами (1). Вариация по α даёт уравнение поля:

α̈ + m² α = ½ Σ_i ξ_i I_i⁰ ω_i². (5)

Система (4)–(5) замкнута: ω_i(t) входят в источник для α, а α(t) входит в моменты инерции. Энергия

E = ½ Σ_i I_i(α) ω_i² + ½ α̇² + ½ m² α² (6)

сохраняется, что проверяется численно.

2.3. Условие бифуркации

Граница устойчивости промежуточной оси определяется знаком произведения

(I₂(α) − I₁(α))(I₃(α) − I₂(α)). (7)

При I₁ < I₂ < I₃ произведение положительно, и ось I₂ неустойчива. Поскольку I_i зависят от α(t), граница является динамической. Критическое значение поля, при котором происходит инверсия порядка I₁ и I₂:

α_crit = (I₂⁰ − I₁⁰) / (ξ(I₁⁰ + I₂⁰)). (8)

Состояние с инверсией достижимо при α > α_crit; оно не запрещено линейной связью, но требует достаточной амплитуды «дыхания» метрики.

2.4. Частота прецессии и предсказание резонанса

Частота малых колебаний оси волчка вокруг вектора момента:

Ω_prec² = (I₂ − I₃)(I₃ − I₁) ω₃² / (I₁ I₂). (9)

При периодической модуляции I₁(t) уравнение для возмущений сводится к уравнению Матьё:

δω̈₁ + Ω₀²(1 − ε cos(m t)) δω₁ = 0. (10)

Параметрический резонанс возникает при

m ≈ 2Ω_prec / n, n = 1, 2, … (11)

Наиболее сильный эффект — при n = 1, то есть m ≈ 2Ω_prec.

3. Численная реализация

Система (4)–(5) интегрировалась методом Radau (scipy.integrate.solve_ivp) с допусками rtol = 10⁻¹⁰, atol = 10⁻¹². Параметры: I₁⁰ = 1.0×10⁻⁶ кг·м², I₂⁰ = 1.2×10⁻⁶ кг·м², I₃⁰ = I₁⁰ + I₂⁰, ξ₁ = 0.1, ξ₂ = −0.1, ξ₃ = (ξ₁I₁⁰ + ξ₂I₂⁰)/I₃⁰. Начальные условия: ω₁ = 0.1, ω₂ = 0.05, ω₃ = 100 рад/с, α = 0, α̇ = 0.

Для сканирования по m использовался диапазон m ∈ [0.5Ω_prec, 4Ω_prec] с 30 значениями. Для каждого m вычислялась максимальная амплитуда ω₁ за время симуляции. Контрольный расчёт проводился при ξ₁ = ξ₂ = 0 (классический предел).

Для физичного режима использовались параметры: I₁⁰ = 1.0 кг·м², I₂⁰ = 1.2 кг·м², ω₃ = 10 рад/с, m = 20 рад/с, ξ₁ = 0.1, ξ₂ = −0.1, время симуляции 200 с.

4. Результаты

4.1. Сохранение энергии

При m = 200 рад/с отклонение энергии за 10 с составило 3.8×10⁻¹⁴ относительно начального значения. Это подтверждает корректность реализации уравнений (4)–(5): ошибка в знаке или пропуск члена İ_i ω_i привели бы к дрейфу на уровне 10⁻³–10⁻⁵.

4.2. Сканирование по частоте модуляции

На зависимости максимальной амплитуды ω₁ от m обнаружен пик при m ≈ 2Ω_prec = 200 рад/с для метрической модели (ξ = 0.1). В классическом пределе (ξ = 0) пик отсутствует: кривая плоская.

Рис. 1.Зависимость максимальной амплитуды ω₁ от частоты модуляции m. Синяя кривая — метрическая модель (ξ = 0.1), оранжевая — классический предел (ξ = 0). Красная пунктирная линия — предсказанная частота 2Ω_prec = 200 рад/с. Пик присутствует только у метрической модели.

 Это указывает на то, что резонанс является следствием метрической связи, а не артефактом интегратора.

4.3. Физичный режим

При I ~ 1 кг·м², ω₃ = 10 рад/с, m = 20 рад/с, ξ = 0.1 получено:

• max|α| = 6.46×10⁻³ — поле физично, инерция меняется на 0.06%;
• min I₁(α) = 0.9994 — момент инерции положителен на всём интервале;
• max|ω₁| = 0.117 рад/с — модуляция ~1% от ω₃;
• отношение числа нулевых переходов α и ω₁ = 2.00 — подпись параметрического резонанса (отклик на субгармонике модуляции).

Рис. 2.Динамика поля α(t) при физичных параметрах (I ~ 1 кг·м², ω₃ = 10 рад/с, m = 20 рад/с, ξ = 0.1). Амплитуда α ~ 6×10⁻³, поле колеблется около нуля; момент инерции I₁(α) остаётся положительным на всём интервале.

Рис. 3.Угловые скорости ω₁, ω₂, ω₃ для того же режима. ω₂ демонстрирует кувырки (эффект Джанибекова), ω₃ стабильна, ω₁ модулируется в такт с α (отношение частот 2:1).

5. Обсуждение

Результаты показывают, что метрическая модель (1)–(5) внутренне согласована (энергия сохраняется), предсказательна (резонанс на 2Ω_prec, отсутствующий в классике) и физична при разумных параметрах (α ~ 6×10⁻³, I₁ > 0). Это косвенное подтверждение того, что рамка способна давать нетривиальные следствия, согласующиеся с физикой.

Рис. 4. Сравнение амплитуды α в физичном (ξ = 0.1, m = 200) и демонстрационном (ξ = 1, m = 0.05) режимах. В физичном режиме α ~ 10⁻⁹ (плоская линия), в демонстрационном — α ~ 10⁻² (видимые колебания). Это иллюстрирует разницу масштабов: эффект реален, но при физичных ξ мал.

Однако работа не доказывает реальность поля α. Форма I_i(α) = I_i⁰(1 + ξ_i α) является анзацем; параметры ξ_i и m феноменологические и не выведены из первых принципов. Их значения должны определяться экспериментально — например, по положению и ширине резонансного пика.

Разделительная диагностика. Классический параметрический резонанс требует внешней вибрации корпуса (I_i(t) = I_i⁰(1 + ε cos νt)). Метрический резонанс самодостаточен: α колеблется из-за вращения, и это же вращение испытывает резонанс. Эксперимент: волчок в вакууме, без внешнего привода; если резонанс на 2Ω_prec возникает сам — это подтверждение модели.

6. Ограничения

1. Феноменологический статус параметров. ξ_i и m не выведены из первых принципов; их значения не известны.
2. Анзац для I_i(α). Линейная форма — простейшее предположение; возможны нелинейные связи.
3. Отсутствие эксперимента. Численная симуляция подтверждает внутреннюю согласованность и предсказание, но не реальность поля α.
4. Идеализация. Модель не учитывает диссипацию, конечную жёсткость тела, внешние поля.

7. Заключение

Предложена минимальная метрическая модель динамики твёрдого тела с динамическими моментами инерции. Уравнения выведены из вариационного принципа в формализме Эйлера–Пуанкаре. Численная симуляция показала: (1) энергия сохраняется с точностью 3.8×10⁻¹⁴; (2) при m ≈ 2Ω_prec возникает параметрический резонанс с отношением частот 2:1, отсутствующий в классическом пределе; (3) при физичных параметрах эффект умеренный и не требует нефизичных условий. Это косвенное подтверждение работоспособности метрической модели как физической конструкции. Прямое подтверждение (реальность α, значения ξ_i и m) требует эксперимента.

Мы уточнили не только Ω_prec, а связь между Ω_prec и α. В классике этой связи нет. В нашей модели Ω_prec² зависит от I_i(α), а I_i(α) зависит от α(t). Значит, Ω_prec «дышит» вместе с метрикой. И это «дыхание» даёт резонанс.

Так что краткая формула: Ω_prec(α) вместо Ω_prec = const. И это — уточнение параметра, которое проверяемо.

Список литературы

1. Алексеенко А. Эффект Джанибекова как метрический срыв: полная лекция с формулами // Цикл «Рецепт хрупкого бессмертия».
2. Алексеенко А. Вселенная, которая дышит: уравнения Фридмана как мост между ОТО и метрической теорией // Цикл «Рецепт хрупкого бессмертия».
3. Алексеенко А. Цвет как метрика: изотопы, лестница состояний, прочность связей и толщина // Цикл «Рецепт хрупкого бессмертия».
4. Алексеенко А. Квантование метрики граничными условиями: эффект Казимира, лестница состояний и рубцевание // Цикл «Рецепт хрупкого бессмертия».
5. Landau L.D., Lifshitz E.M. Mechanics. — Course of Theoretical Physics, Vol. 1.
6. Arnold V.I. Mathematical Methods of Classical Mechanics.
7. McLachlan N.W. Theory and Application of Mathieu Functions.
8. Brans C., Dicke R.H. Mach's Principle and a Relativistic Theory of Gravitation // Phys. Rev. 124, 925 (1961).

Карта достоверности (для внутреннего протокола)

Классическое объяснение эффекта Джанибекова — факт, полное для механизма.
Вывод уравнений (4)–(5) из лагранжиана (2) — структурный вывод; форма корректна, параметры феноменологические.
Сохранение энергии 3.8×10⁻¹⁴ — численный результат, подтверждает корректность реализации.
Пик при 2Ω_prec только у ξ=0.1 — численный результат, подтверждает предсказание.
Отношение частот 2:1 — численный результат, подпись параметрического резонанса.
Физичность при ξ=0.1, m=20 — численный результат, модель не требует нефизичных параметров.
Косвенное подтверждение работоспособности модели — следствие трёх результатов; не является доказательством реальности α.
Реальность α и значения ξ_i, m — гипотеза, требующая эксперимента.

Приложение 1: Код.

Вот модели кода для запуска.*

Первая модель: (Модель проверяет качество формул)

Модель качества кода

# -*- coding: utf-8 -*-
"""
Симуляция метрической модели вращающегося тела.
Требует: numpy, scipy, matplotlib.
Запуск: python metric_sim.py
"""

import sys
import numpy as np
from scipy.integrate import solve_ivp

# --- Matplotlib: пробуем интерактивный backend, если не выходит — сохраняем в файл ---
import matplotlib
try:
   matplotlib.use('TkAgg')  # стандартный backend Windows
except Exception:
   matplotlib.use('Agg')    # если нет GUI — рисуем в файлы
import matplotlib.pyplot as plt

# --- Параметры ---
I1_0 = 1.0e-6
I2_0 = 1.2e-6
I3_0 = I1_0 + I2_0
xi1 = 0.1
xi2 = -0.1
xi3 = (xi1 * I1_0 + xi2 * I2_0) / I3_0
m = 200.0  # рад/с; для видимого эффекта попробуй 2.0

def I1(alpha): return I1_0 * (1 + xi1 * alpha)
def I2(alpha): return I2_0 * (1 + xi2 * alpha)
def I3(alpha): return I3_0 * (1 + xi3 * alpha)

def equations(t, y):
   w1, w2, w3, alpha, dalpha = y
   I1a, I2a, I3a = I1(alpha), I2(alpha), I3(alpha)
   dw1 = ((I2a - I3a) * w2 * w3 - xi1 * I1_0 * dalpha * w1) / I1a
   dw2 = ((I3a - I1a) * w3 * w1 - xi2 * I2_0 * dalpha * w2) / I2a
   dw3 = ((I1a - I2a) * w1 * w2 - xi3 * I3_0 * dalpha * w3) / I3a
   ddalpha = 0.5 * (xi1 * I1_0 * w1**2 + xi2 * I2_0 * w2**2 + xi3 * I3_0 * w3**2) - m**2 * alpha
   return [dw1, dw2, dw3, dalpha, ddalpha]

def energy(y):
   w1, w2, w3, alpha, dalpha = y
   return 0.5 * (I1(alpha) * w1**2 + I2(alpha) * w2**2 + I3(alpha) * w3**2) \
         + 0.5 * dalpha**2 + 0.5 * m**2 * alpha**2

# --- Начальные условия ---
y0 = [0.1, 0.05, 100.0, 0.0, 0.0]
t_span = (0.0, 10.0)      # секунд; для быстрого теста поставь (0.0, 1.0)
t_eval = np.linspace(t_span[0], t_span[1], 10000)

print("Запуск solve_ivp...")
sol = solve_ivp(
   equations, t_span, y0, t_eval=t_eval,
   method='Radau', rtol=1e-10, atol=1e-12
)

if not sol.success:
   print("Интегратор не сошёлся:", sol.message)
   sys.exit(1)

# --- Проверка сохранения энергии (векторизовано) ---
w1, w2, w3, alpha, dalpha = sol.y
E = (0.5 * (I1(alpha) * w1**2 + I2(alpha) * w2**2 + I3(alpha) * w3**2)
   + 0.5 * dalpha**2 + 0.5 * m**2 * alpha**2)

drift = np.max(np.abs(E - E[0])) / E[0]
print(f"Отклонение энергии: {drift:.2e}")

# --- График 1: угловые скорости ---
fig1, ax1 = plt.subplots(figsize=(10, 6))
ax1.plot(sol.t, sol.y[0], label='ω₁')
ax1.plot(sol.t, sol.y[1], label='ω₂')
ax1.plot(sol.t, sol.y[2] / 100.0, label='ω₃/100')
ax1.set_xlabel('t, с'); ax1.set_ylabel('ω, рад/с')
ax1.set_title('Динамика угловых скоростей')
ax1.legend(); ax1.grid()
fig1.tight_layout()
fig1.savefig('omega_plot.png', dpi=150)
print("Сохранено: omega_plot.png")

# --- График 2: поле alpha ---
fig2, ax2 = plt.subplots(figsize=(8, 4))
ax2.plot(sol.t, sol.y[3], label='α')
ax2.set_xlabel('t, с'); ax2.set_ylabel('α')
ax2.set_title('Динамика метрического поля α')
ax2.legend(); ax2.grid()
fig2.tight_layout()
fig2.savefig('alpha_plot.png', dpi=150)
print("Сохранено: alpha_plot.png")

# --- Пытаемся показать окна (если backend интерактивный) ---
try:
   plt.show()
except Exception as e:
   print("Окна не открылись (это нормально для backend 'Agg'):", e)

print("Готово.")

Код Python 3

Вторая модель: (Наглядно показывающая резонанс)

Модель для наглядности резонанса

# -*- coding: utf-8 -*-
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# --- Параметры ---
I1_0 = 1.0e-6
I2_0 = 1.2e-6
I3_0 = I1_0 + I2_0

def run_sim(m, xi1, xi2, t_end=5.0, n_points=5000):
    xi3 = (xi1*I1_0 + xi2*I2_0) / I3_0

    def I1(a): return I1_0*(1 + xi1*a)
    def I2(a): return I2_0*(1 + xi2*a)
    def I3(a): return I3_0*(1 + xi3*a)

    def eq(t, y):
       w1, w2, w3, a, da = y
      I1a, I2a, I3a = I1(a), I2(a), I3(a)
      dw1 = ((I2a - I3a)*w2*w3 - xi1*I1_0*da*w1) / I1a
      dw2 = ((I3a - I1a)*w3*w1 - xi2*I2_0*da*w2) / I2a
      dw3 = ((I1a - I2a)*w1*w2 - xi3*I3_0*da*w3) / I3a
      dda = 0.5*(xi1*I1_0*w1**2 + xi2*I2_0*w2**2 + xi3*I3_0*w3**2) - m**2*a
      return [dw1, dw2, dw3, da, dda]

   y0 = [0.1, 0.05, 100.0, 0.0, 0.0]
   t_eval = np.linspace(0, t_end, n_points)
   sol = solve_ivp(eq, (0, t_end), y0, t_eval=t_eval,
                  method='Radau', rtol=1e-9, atol=1e-11)
   return sol

# --- Ω_prec при alpha=0 ---
def omega_prec(w3=100.0):
   return np.sqrt(abs((I2_0 - I3_0)*(I3_0 - I1_0)) / (I1_0*I2_0)) * w3

Op = omega_prec()
print(f"Ω_prec = {Op:.2f} рад/с,  2Ω_prec = {2*Op:.2f} рад/с")

# --- Шаг 1: один прогон с физичными параметрами, три панели ---
xi1, xi2 = 0.1, -0.1
m_test = 2*Op
sol = run_sim(m_test, xi1, xi2, t_end=5.0)

w1, w2, w3, alpha, _ = sol.y
print(f"max|alpha| = {np.max(np.abs(alpha)):.3e}  (должно быть < {1/xi1:.1f})")

fig, axes = plt.subplots(3, 1, figsize=(10, 8), sharex=True)
axes[0].plot(sol.t, w1); axes[0].set_ylabel('ω₁'); axes[0].grid()
axes[1].plot(sol.t, w2); axes[1].set_ylabel('ω₂'); axes[1].grid()
axes[2].plot(sol.t, w3); axes[2].set_ylabel('ω₃'); axes[2].set_xlabel('t, с'); axes[2].grid()
fig.suptitle(f'Угловые скорости при m = 2Ω_prec = {m_test:.1f}')
fig.tight_layout(); fig.savefig('omega_3panels.png', dpi=150)

fig2, ax2 = plt.subplots(figsize=(8, 4))
ax2.plot(sol.t, alpha); ax2.set_xlabel('t, с'); ax2.set_ylabel('α'); ax2.grid()
ax2.set_title(f'Поле α при m = {m_test:.1f}')
fig2.tight_layout(); fig2.savefig('alpha_single.png', dpi=150)

# --- Шаг 2: сканирование по m ---
print("Сканирование по m...")
m_values = np.linspace(0.5*Op, 4*Op, 30)
amp_xi = []
amp_classic = []

for m in m_values:
   # с метрической связью
   s = run_sim(m, 0.1, -0.1, t_end=5.0, n_points=3000)
   amp_xi.append(np.max(np.abs(s.y[0])))
   # классика (xi=0)
   s0 = run_sim(m, 0.0, 0.0, t_end=5.0, n_points=3000)
   amp_classic.append(np.max(np.abs(s0.y[0])))

amp_xi = np.array(amp_xi)
amp_classic = np.array(amp_classic)

fig3, ax3 = plt.subplots(figsize=(10, 5))
ax3.plot(m_values, amp_xi, 'o-', label='с метрической связью (ξ=0.1)')
ax3.plot(m_values, amp_classic, 's--', label='классика (ξ=0)')
ax3.axvline(2*Op, color='red', linestyle=':', label=f'2Ω_prec = {2*Op:.1f}')
ax3.set_xlabel('m, рад/с'); ax3.set_ylabel('max |ω₁|'); ax3.grid(); ax3.legend()
ax3.set_title('Амплитуда ω₁ vs частота «дыхания» m')
fig3.tight_layout(); fig3.savefig('resonance_scan.png', dpi=150)

print("Сохранено: omega_3panels.png, alpha_single.png, resonance_scan.png")
print("Готово.")

Код Python 3

Третья модель: (Реальная физическая)

Физическая модель

# -*- coding: utf-8 -*-
"""
Физичный демо-режим: alpha физична (|alpha| < 0.1), резонанс 2:1 виден.
Тело большое (I ~ 1 кг·м²), omega3 = 10 рад/с, m = 2*Omega_prec = 20 рад/с.
"""
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# --- Параметры: большое тело, физичные alpha ---
I1_0 = 1.0      # кг·м² (не гайка, а большое тело)
I2_0 = 1.2
I3_0 = I1_0 + I2_0

xi1 = 0.1
xi2 = -0.1
xi3 = (xi1*I1_0 + xi2*I2_0) / I3_0

omega3_0 = 10.0  # рад/с
m = 20.0          # рад/с = 2*Omega_prec (Omega_prec = omega3 = 10)

t_end = 200.0
n_points = 10000

def I1(a): return I1_0*(1 + xi1*a)
def I2(a): return I2_0*(1 + xi2*a)
def I3(a): return I3_0*(1 + xi3*a)

def eq(t, y):
    w1, w2, w3, a, da = y
    I1a, I2a, I3a = I1(a), I2(a), I3(a)
    dw1 = ((I2a - I3a)*w2*w3 - xi1*I1_0*da*w1) / I1a
    dw2 = ((I3a - I1a)*w3*w1 - xi2*I2_0*da*w2) / I2a
    dw3 = ((I1a - I2a)*w1*w2 - xi3*I3_0*da*w3) / I3a
    dda = 0.5*(xi1*I1_0*w1**2 + xi2*I2_0*w2**2 + xi3*I3_0*w3**2) - m**2*a
    return [dw1, dw2, dw3, da, dda]

y0 = [0.1, 0.05, omega3_0, 0.0, 0.0]
t_eval = np.linspace(0, t_end, n_points)

sol = solve_ivp(eq, (0, t_end), y0, t_eval=t_eval,
    method='Radau', rtol=1e-10, atol=1e-12)

if not sol.success:
    print("Интегратор не сошёлся:", sol.message)
    exit(1)

w1, w2, w3, alpha, dalpha = sol.y
print(f"max|alpha| = {np.max(np.abs(alpha)):.3e}  (должно быть < 0.1)")
print(f"max|w1|    = {np.max(np.abs(w1)):.3e}")
print(f"min I1(alpha) = {np.min(I1(alpha)):.3e}  (должно быть > 0)")

# --- График 1: три панели omega ---
fig, axes = plt.subplots(3, 1, figsize=(10, 8), sharex=True)
axes[0].plot(sol.t, w1); axes[0].set_ylabel('ω₁'); axes[0].grid()
axes[1].plot(sol.t, w2); axes[1].set_ylabel('ω₂'); axes[1].grid()
axes[2].plot(sol.t, w3); axes[2].set_ylabel('ω₃'); axes[2].set_xlabel('t, с'); axes[2].grid()
fig.suptitle(f'Физичный режим: m = {m} рад/с = 2Ω_prec, max|α| = {np.max(np.abs(alpha)):.2e}')
fig.tight_layout(); fig.savefig('phys_omega_3panels.png', dpi=150)

# --- График 2: alpha ---
fig2, ax2 = plt.subplots(figsize=(8, 4))
ax2.plot(sol.t, alpha); ax2.set_xlabel('t, с'); ax2.set_ylabel('α'); ax2.grid()
ax2.set_title(f'α(t): max = {np.max(np.abs(alpha)):.2e}')
fig2.tight_layout(); fig2.savefig('phys_alpha.png', dpi=150)

# --- Оценка частот (нулевые переходы) ---
def count_zero_crossings(signal):
    return np.sum(np.diff(np.sign(signal)) != 0)

n_alpha = count_zero_crossings(alpha - np.mean(alpha))
n_w1 = count_zero_crossings(w1 - np.mean(w1))
print(f"Нулевых переходов: alpha = {n_alpha}, w1 = {n_w1}, отношение = {n_alpha/max(n_w1,1):.2f}")

print("Сохранено: phys_omega_3panels.png, phys_alpha.png")
print("Готово.")

Код Python 3

*может требоваться правка синтаксиса кода через нейросеть перед запуском.