НПтаМ_Ат_Мірошниченко

from __future__ import annotations

import math
from pathlib import Path

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

S_VALUE = 0.7
CASES = [
    ('a', 0.01, r'$\tau=0.01$'),
    ('b', 1.0,  r'$\tau=1$'),
    ('c', 100.0, r'$\tau=100$'),
]


def rhs(_, z, s, tau):
    eta, h = z
    return np.array([
        -eta + h,
        (s * eta - h * (1.0 + eta**2)) / tau,
    ])


def lyapunov_roots(s, tau):
    disc = (tau + 1.0) ** 2 - 4.0 * tau * (1.0 - s)
    root = math.sqrt(disc)
    lam1 = (-(tau + 1.0) + root) / (2.0 * tau)
    lam2 = (-(tau + 1.0) - root) / (2.0 * tau)
    return lam1, lam2


def build_initial_conditions():
    initials = []
    for eta0 in np.linspace(-1.2, 1.2, 7):
        for h0 in np.linspace(-1.0, 1.0, 6):
            if abs(eta0) < 0.08 and abs(h0) < 0.08:
                continue
            initials.append((eta0, h0))
    return initials


def make_plot(output, s=S_VALUE):
    etas = np.linspace(-1.4, 1.4, 250)
    hs = np.linspace(-1.15, 1.15, 250)
    E, H = np.meshgrid(etas, hs)
    initials = build_initial_conditions()

    fig, axes = plt.subplots(1, 3, figsize=(13.5, 4.3), constrained_layout=True)

    for ax, (_, tau, title) in zip(axes, CASES):
        U = -E + H
        V = (s * E - H * (1.0 + E**2)) / tau
        ax.streamplot(etas, hs, U, V, density=1.1, color='0.75', linewidth=0.75, arrowsize=0.75)

        t_max = 80 if tau < 10 else 250
        t_eval = np.linspace(0, t_max, 2200)
        for idx, z0 in enumerate(initials):
            sol = solve_ivp(rhs, [0, t_max], z0, t_eval=t_eval, args=(s, tau), rtol=1e-7, atol=1e-9)
            lw = 1.0 if idx % 5 else 1.2
            ax.plot(sol.y[0], sol.y[1], color='black', linewidth=lw, alpha=0.85)

        ax.plot(etas, etas, linestyle='--', linewidth=1.4, color='black', label=r'$\dot{\eta}=0$')
        ax.plot(etas, s * etas / (1 + etas**2), linestyle=':', linewidth=1.8, color='black', label=r'$\dot{h}=0$')
        ax.scatter([0], [0], s=28, color='black', zorder=5)
        ax.text(0.05, 0.05, 'D', fontsize=11, transform=ax.transAxes)
        ax.set_title(title, fontsize=13, pad=8)
        ax.set_xlabel(r'$\eta$', fontsize=12)
        if ax is axes[0]:
            ax.set_ylabel(r'$h$', fontsize=12)
        ax.set_xlim(-1.35, 1.35)
        ax.set_ylim(-1.05, 1.05)
        ax.grid(alpha=0.18)

    handles = [
        plt.Line2D([0], [0], linestyle='--', linewidth=1.4, color='black', label=r'$\dot{\eta}=0$'),
        plt.Line2D([0], [0], linestyle=':', linewidth=1.8, color='black', label=r'$\dot{h}=0$'),
        plt.Line2D([0], [0], color='black', linewidth=1.1, label='траєкторії'),
    ]
    fig.legend(handles=handles, loc='lower center', ncol=3, frameon=False, bbox_to_anchor=(0.5, -0.03), fontsize=11)
    fig.savefig(output, dpi=220, bbox_inches='tight')
    plt.close(fig)


def main():
    output = Path('phase_portraits_task2.png')
    make_plot(output)

    print('Завдання №2: неупорядкована фаза для випадку τ_s << τ_η, τ_h; s =', S_VALUE)
    for label, tau, _ in CASES:
        lam1, lam2 = lyapunov_roots(S_VALUE, tau)
        print(f'{label}: tau={tau:>6g}; lambda1={lam1: .6f}; lambda2={lam2: .6f}; тип точки D -> стійкий вузол')
    print(f'Зображення збережено у файл: {output.resolve()}')


if __name__ == '__main__':
    main()