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

import numpy as np, math
import matplotlib.pyplot as plt
import mpmath as mp
from pathlib import Path



k = 1.0
eta_t = 0.5
theta = eta_t**2

# ---------- Part 1 ----------
etas = np.linspace(-2.5, 2.5, 600)
Se_vals = [0.8, 1.0, 1.6, 2.5]

fig = plt.figure(figsize=(8, 5))
for Se in Se_vals:
    S = Se / (1 + etas**2)
    h = Se * etas / (1 + etas**2)
    plt.plot(etas, S, label=fr'$S(\eta),\ S_e={Se}$')
    plt.plot(etas, h, linestyle='--', label=fr'$h(\eta),\ S_e={Se}$')
plt.axhline(0, linewidth=0.8, color='black')
plt.axvline(0, linewidth=0.8, color='black')
plt.xlabel(r'$\eta$')
plt.ylabel('Значення')
plt.title('Частина 1: залежності $S(\eta)$ та $h(\eta)$')
plt.legend(fontsize=8, ncol=2)
plt.tight_layout()
plt.savefig('part1_dependencies.png', dpi=200)
plt.close(fig)

fig = plt.figure(figsize=(8, 5))
for Se in [0.8, 1.0, 1.3, 2.0]:
    V = 0.5 * etas**2 - 0.5 * Se * np.log(1 + etas**2)
    plt.plot(etas, V, label=fr'$S_e={Se}$')
plt.axhline(0, linewidth=0.8, color='black')
plt.axvline(0, linewidth=0.8, color='black')
plt.xlabel(r'$\eta$')
plt.ylabel(r'$V(\eta)$')
plt.title('Частина 1: синергетичний потенціал')
plt.legend()
plt.tight_layout()
plt.savefig('part1_potential.png', dpi=200)
plt.close(fig)

Se_grid = np.linspace(0, 3, 601)
eta_pos = np.where(Se_grid >= 1, np.sqrt(Se_grid - 1), np.nan)
eta_neg = np.where(Se_grid >= 1, -np.sqrt(Se_grid - 1), np.nan)
fig = plt.figure(figsize=(8, 5))
plt.plot(Se_grid, np.zeros_like(Se_grid), label=r'$\eta_0=0$')
plt.plot(Se_grid, eta_pos, label=r'$\eta_0=+\sqrt{S_e-1}$')
plt.plot(Se_grid, eta_neg, label=r'$\eta_0=-\sqrt{S_e-1}$')
plt.axvline(1, linestyle='--', color='gray', label=r'$S_c=1$')
plt.xlabel(r'$S_e$')
plt.ylabel(r'$\eta_0$')
plt.title('Частина 1: стаціонарний параметр порядку')
plt.legend()
plt.tight_layout()
plt.savefig('part1_stationary.png', dpi=200)
plt.close(fig)

# ---------- Part 2 ----------
def V2(eta, Se):
    return 0.5*eta**2 + 0.5*k*eta_t**2*np.log(1 + (eta/eta_t)**2) - 0.5*Se*np.log(1 + eta**2)

def stationary_x_roots(Se):
    B = (1+k)*eta_t**2 + (1-Se)
    C = eta_t**2*(1+k-Se)
    disc = B*B - 4*C
    if disc < -1e-12:
        return []
    disc = max(disc, 0.0)
    roots = [(-B - math.sqrt(disc))/2, (-B + math.sqrt(disc))/2]
    return sorted([float(r) for r in roots if r >= -1e-12])

def d2V(eta, Se):
    return 1 + k*(1-(eta/eta_t)**2)/(1+(eta/eta_t)**2)**2 - Se*(1-eta**2)/(1+eta**2)**2

Se_plateau = 1 - theta + k*theta + 2*math.sqrt(k*theta*(1-theta))
Se_spinodal0 = 1 + k

def coexistence_eq(y):
    Se = (1+y)*(1 + k/(1+y/theta))
    return y + k*theta*mp.log(1+y/theta) - Se*mp.log(1+y)

ys = np.linspace(1e-5, 10, 5000)
vals = [float(coexistence_eq(y)) for y in ys]
root_y = None
for i in range(len(ys)-1):
    if vals[i] == 0 or vals[i+1] == 0 or vals[i]*vals[i+1] < 0:
        root_y = float(mp.findroot(lambda z: coexistence_eq(z), (ys[i], ys[i+1])))
        break
Se_coexist = float((1+root_y)*(1 + k/(1+root_y/theta)))


fig = plt.figure(figsize=(8, 5.3))
for Se, label in [(1.7, r'$S_e<S_e^{pl}$'), (Se_plateau, r'$S_e=S_e^{pl}$'), (Se_coexist, r'$S_e=S_e^*$'), (2.0, r'$S_e=S_e^{(0)}$')]:
    plt.plot(etas, V2(etas, Se), label=label)
plt.axhline(0, linewidth=0.8, color='black')
plt.axvline(0, linewidth=0.8, color='black')
plt.xlabel(r'$\eta$')
plt.ylabel(r'$V(\eta)$')
plt.title('Частина 2: синергетичний потенціал за різних $S_e$')
plt.legend(fontsize=9)
plt.tight_layout()
plt.savefig('part2_potential.png', dpi=200)
plt.close(fig)


Se_grid = np.linspace(1.4, 2.2, 500)
stable_pos, stable_neg, unstable_pos, unstable_neg = [], [], [], []
for Se in Se_grid:
    roots = stationary_x_roots(Se)
    st_pos = st_neg = un_pos = un_neg = np.nan
    for x in roots:
        if x < 1e-10:
            continue
        eta = math.sqrt(x)
        if d2V(eta, Se) > 0:
            st_pos, st_neg = eta, -eta
        else:
            un_pos, un_neg = eta, -eta
    stable_pos.append(st_pos); stable_neg.append(st_neg)
    unstable_pos.append(un_pos); unstable_neg.append(un_neg)
fig = plt.figure(figsize=(8, 5.3))
zero_stable = np.where(Se_grid <= Se_spinodal0, 0.0, np.nan)
zero_unstable = np.where(Se_grid >= Se_spinodal0, 0.0, np.nan)
plt.plot(Se_grid, zero_stable, label='стійка нульова гілка')
plt.plot(Se_grid, stable_pos, label='стійка ненульова гілка')
plt.plot(Se_grid, stable_neg)
plt.plot(Se_grid, zero_unstable, linestyle='--', label='нестійка нульова гілка')
plt.plot(Se_grid, unstable_pos, linestyle='--', label='нестійка ненульова гілка')
plt.plot(Se_grid, unstable_neg, linestyle='--')
plt.axvline(Se_plateau, linestyle=':', color='gray', label=fr'$S_e^{{pl}}\approx {Se_plateau:.3f}$')
plt.axvline(Se_coexist, linestyle='-.', color='gray', label=fr'$S_e^*\approx {Se_coexist:.3f}$')
plt.axvline(Se_spinodal0, linestyle='--', color='black', label=fr'$S_e^{{(0)}}={Se_spinodal0:.1f}$')
plt.xlabel(r'$S_e$')
plt.ylabel(r'$\eta_0$')
plt.title('Частина 2: стаціонарні стани і область гістерезису')
plt.legend(fontsize=8, ncol=2)
plt.tight_layout()
plt.savefig('part2_stationary.png', dpi=200)
plt.close(fig)

print('done', Se_plateau, Se_coexist, Se_spinodal0)