НПтаМ_ПР№2_Мірошниченко

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

def system(t, y, s, tau):
    eta, S = y
    return [-eta * (1 - s * S), (1 - S * (1 + eta**2)) / tau]

def lambda_D(s, tau):
    a = (s - 1) - 1 / tau
    disc = 1 + 4 * (1 / tau) * (s - 1) / (a * a)
    r = math.sqrt(disc)
    return 0.5 * a * (1 + r), 0.5 * a * (1 - r)

def lambda_O(s, tau):
    disc = 1 - 8 * tau * (s - 1) / (s * s)
    pref = -s / (2 * tau)
    if disc >= 0:
        r = math.sqrt(disc)
        return pref * (1 + r), pref * (1 - r)
    r = math.sqrt(-disc)
    return complex(pref, -pref * r), complex(pref, pref * r)

def plot_phase_portraits(s, taus, ordered, out_path):
    fig, axes = plt.subplots(1, 3, figsize=(14, 4.5), constrained_layout=True)
    labels = ["а", "б", "в"]

    for ax, tau, lab in zip(axes, taus, labels):
        if ordered:
            Smin, Smax = 0.0, 1.5
            etamin, etamax = -1.3, 1.3
        else:
            Smin, Smax = 0.0, 2.0
            etamin, etamax = -1.1, 1.1

        Sg = np.linspace(Smin, Smax, 120)
        Eg = np.linspace(etamin, etamax, 120)
        X, Y = np.meshgrid(Sg, Eg)

        U = (1 - X * (1 + Y**2)) / tau
        V = -Y * (1 - s * X)
        speed = np.sqrt(U * U + V * V)

        ax.streamplot(
            Sg, Eg, U / (speed + 1e-9), V / (speed + 1e-9),
            density=1.0, linewidth=0.75, arrowsize=0.7, color="0.25"
        )

        eta_nc = np.linspace(etamin, etamax, 400)
        ax.plot(1 / (1 + eta_nc**2), eta_nc, "--", lw=1, color="0.5")
        ax.axhline(0, ls=":", lw=1, color="0.5")
        ax.axvline(1 / s, ls=":", lw=1, color="0.7")

        if ordered:
            eta0 = math.sqrt(s - 1)
            inits = [
                (1.1, 0.1), (-1.1, 0.1), (1.1, 1.4), (-1.1, 1.4),
                (0.0, 1.4), (0.4, 0.5), (-0.4, 0.5),
                (0.9, 1.0), (-0.9, 1.0),
            ]
        else:
            inits = [
                (1.0, 0.1), (-1.0, 0.1), (1.0, 1.9), (-1.0, 1.9),
                (0.0, 1.8), (0.6, 0.5), (-0.6, 0.5),
                (0.6, 1.2), (-0.6, 1.2),
            ]

        t_end = 25 if tau == 0.01 else (40 if tau == 1 else 500)
        max_step = 0.05 if tau == 0.01 else (0.1 if tau == 1 else 1.0)

        for e0, s0 in inits:
            sol = solve_ivp(
                system, [0, t_end], [e0, s0],
                args=(s, tau), method="Radau",
                max_step=max_step, rtol=1e-5, atol=1e-8
            )
            ax.plot(sol.y[1], sol.y[0], color="k", lw=0.8)

        ax.plot([1], [0], "ko", ms=3)
        ax.text(1.02, -0.08, "D", fontsize=10)

        if ordered:
            ax.plot([1 / s, 1 / s], [eta0, -eta0], "ko", ms=3)
            ax.text(1 / s + 0.03, eta0 + 0.03, "O+", fontsize=10)
            ax.text(1 / s + 0.03, -eta0 - 0.12, "O-", fontsize=10)

        ax.set_title(f"{lab}) τ = {tau:g}")
        ax.set_xlim(Smin, Smax)
        ax.set_ylim(etamin, etamax)
        ax.set_xlabel("S")
        ax.set_ylabel("η")

    fig.savefig(out_path, dpi=220, bbox_inches="tight")
    plt.close(fig)

if __name__ == "__main__":
    os.makedirs("phase_output", exist_ok=True)
    taus = [0.01, 1, 100]
    plot_phase_portraits(0.7, taus, ordered=False, out_path="phase_output/part1.png")
    plot_phase_portraits(1.5, taus, ordered=True, out_path="phase_output/part2.png")
    print("λ_D for s = 0.7:")
    for tau in taus:
        print(tau, lambda_D(0.7, tau))
    print("λ_O for s = 1.5:")
    for tau in taus:
        print(tau, lambda_O(1.5, tau))