API tools faq
paste
mirosh111000

Мірошниченко_КМЗПМ_ЛР№4

mirosh111000
Oct 30th, 2025
1,086
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 8.14 KB | None | 0 0
raw download clone embed print report
  1. from __future__ import annotations
  2. import numpy as np
  3. import matplotlib.pyplot as plt
  4. from matplotlib.patches import ConnectionPatch
  5.  
  6.  
  7.  
  8. BETA = 0.9
  9. EPS = 4.0
  10. RHO = 0.25
  11.  
  12.  
  13. CROSSES = [
  14.     (0.1, 0.75), (0.1, 1.3), (0.1, 1.8),
  15.     (0.2, 0.75), (0.2, 1.3), (0.2, 1.8),
  16. ]
  17.  
  18.  
  19. P_MIN, P_MAX, P_POINTS = 0.0, 0.5, 220
  20. T_MIN, T_MAX, T_POINTS = 0.05, 2.00, 150
  21. K_MAX, K_POINTS = 3.0, 400
  22.  
  23.  
  24.  
  25. def mu_of(c: np.ndarray | float) -> np.ndarray | float:
  26.  
  27.     return c * (1.0 - c)
  28.  
  29.  
  30. def f_c(c: np.ndarray | float, P: float, T: float, beta: float = BETA, eps: float = EPS, Dpar: float = 0.1) -> np.ndarray | float:
  31.  
  32.     return (
  33.         P * (1.0 - c) * (1.0 - beta * c)
  34.         - c * (1.0 - beta * c) * np.exp(-2.0 * eps * c / T)
  35.         - Dpar * c * (1.0 - beta)
  36.     )
  37.  
  38.  
  39. def df_dc(c: float, P: float, T: float, beta: float = BETA, eps: float = EPS, Dpar: float = 0.1) -> float:
  40.  
  41.  
  42.     term1 = P * (-(1.0 + beta) + 2.0 * beta * c)
  43.  
  44.     g  = c - beta * c * c  
  45.     gp = 1.0 - 2.0 * beta * c    
  46.     h  = np.exp(-2.0 * eps * c / T)
  47.     hp = (-2.0 * eps / T) * h
  48.     term2 = -gp * h + (2.0 * eps / T) * g * h
  49.  
  50.     term3 = -Dpar * (1.0 - beta)
  51.     return float(term1 + term2 + term3)
  52.  
  53.  
  54.  
  55.  
  56. def find_c_st(P: float, T: float, beta: float = BETA, eps: float = EPS, Dpar: float = 0.1,
  57.                tol: float = 1e-10, maxit: int = 200) -> float:
  58.  
  59.     a, b = 0.0, 1.0
  60.     fa = f_c(a, P, T, beta, eps, Dpar)
  61.     fb = f_c(b, P, T, beta, eps, Dpar)
  62.     if fa * fb > 0:  
  63.         xs = np.linspace(0.0, 1.0, 2001)
  64.         vals = f_c(xs, P, T, beta, eps, Dpar)
  65.         s = np.sign(vals)
  66.         idx = np.where(s[:-1] * s[1:] <= 0)[0]
  67.         if len(idx) == 0:
  68.  
  69.             return float(xs[np.argmin(np.abs(vals))])
  70.         a, b = float(xs[idx[0]]), float(xs[idx[0] + 1])
  71.         fa, fb = float(vals[idx[0]]), float(vals[idx[0] + 1])
  72.     for _ in range(maxit):
  73.         m = 0.5 * (a + b)
  74.         fm = f_c(m, P, T, beta, eps, Dpar)
  75.         if abs(fm) < tol or (b - a) < tol:
  76.             return float(m)
  77.         if fa * fm <= 0:
  78.             b, fb = m, fm
  79.         else:
  80.             a, fa = m, fm
  81.     return float(0.5 * (a + b))
  82.  
  83.  
  84.  
  85.  
  86. def lambda_k(k: np.ndarray, P: float, T: float, Dpar: float = 0.1,
  87.               beta: float = BETA, eps: float = EPS, rho: float = RHO) -> np.ndarray:
  88.     T_eff = max(T, 1e-9)
  89.     cst = find_c_st(P, T_eff, beta=beta, eps=eps, Dpar=Dpar)
  90.     lam0 = df_dc(cst, P, T_eff, beta=beta, eps=eps, Dpar=Dpar)
  91.     mu   = mu_of(cst)
  92.     return lam0 - k**2 * (1.0 - 2.0 * (eps / T_eff) * mu * (1.0 - (rho**2) * k**2))
  93.  
  94.  
  95.  
  96.  
  97. def max_lambda_over_k(P: float, T: float, Dpar: float) -> float:
  98.     ks = np.linspace(0.0, K_MAX, K_POINTS)
  99.     lams = lambda_k(ks, P, T, Dpar=Dpar)
  100.     return float(np.max(lams))
  101.  
  102.  
  103. def Tcrit_numeric_for_P(P: float, Dpar: float, t_min: float = T_MIN, t_max: float = T_MAX,
  104.                         nT: int = T_POINTS) -> float | float('nan'):
  105.  
  106.     Ts = np.linspace(t_min, t_max, nT)
  107.     phi = np.array([max_lambda_over_k(P, T, Dpar) for T in Ts])
  108.     sgn = np.sign(phi)
  109.     roots = []
  110.     for i in range(len(Ts) - 1):
  111.         if sgn[i] == 0:
  112.             roots.append(Ts[i]); continue
  113.         if sgn[i] * sgn[i + 1] < 0:
  114.             a, b = Ts[i], Ts[i + 1]
  115.             fa, fb = phi[i], phi[i + 1]
  116.  
  117.             for _ in range(40):
  118.                 m = 0.5 * (a + b)
  119.                 fm = max_lambda_over_k(P, m, Dpar)
  120.                 if abs(fm) < 1e-6 or (b - a) < 1e-3:
  121.                     roots.append(m); break
  122.                 if fa * fm <= 0:
  123.                     b, fb = m, fm
  124.                 else:
  125.                     a, fa = m, fm
  126.     if len(roots) == 0:
  127.         return float('nan')
  128.     return float(np.max(roots))
  129.  
  130.  
  131. def Tcrit_curve(P_grid: np.ndarray, Dpar: float) -> tuple[np.ndarray, np.ndarray]:
  132.     T_vals = np.array([Tcrit_numeric_for_P(P, Dpar) for P in P_grid])
  133.     mask = np.isfinite(T_vals)
  134.     return P_grid[mask], T_vals[mask]
  135.  
  136.  
  137.  
  138.  
  139. def plot_PT_main():
  140.     P_grid = np.linspace(P_MIN, P_MAX, P_POINTS)
  141.     P_b, T_b = Tcrit_curve(P_grid, Dpar=0.1)
  142.     P_r, T_r = Tcrit_curve(P_grid, Dpar=0.5)
  143.  
  144.     fig, ax = plt.subplots(figsize=(9.5, 5.5), dpi=130)
  145.     ax.plot(P_b, T_b, 'k-', lw=2.2, label=r'$D_\parallel=0.1$')
  146.     ax.plot(P_r, T_r, 'r--', lw=2.2, label=r'$D_\parallel=0.5$')
  147.     for (p, t) in CROSSES:
  148.         ax.plot(p, t, 'kx', ms=8, mew=1.8)
  149.     ax.set_xlim(P_MIN, P_MAX)
  150.     ax.set_ylim(0.0, T_MAX)
  151.     ax.set_xlabel('P'); ax.set_ylabel('T')
  152.     ax.text(0.04, 1.72, 'B', fontsize=12)
  153.     ax.text(0.44, 0.32,  'A', fontsize=12)
  154.     ax.legend(loc='upper right', frameon=False)
  155.     ax.set_title('Умови формування просторових структур')
  156.     plt.tight_layout(); plt.show()
  157.  
  158.  
  159. def plot_lambda_small(P_val: float, T_val: float):
  160.     ks = np.linspace(0.0, K_MAX, K_POINTS)
  161.     lam_b = lambda_k(ks, P_val, T_val, Dpar=0.1)
  162.     lam_r = lambda_k(ks, P_val, T_val, Dpar=0.5)
  163.     fig, ax = plt.subplots(figsize=(5.3, 4.0), dpi=130)
  164.     ax.plot(ks, lam_b, 'k-', lw=2, label=r'$D_\parallel=0.1$')
  165.     ax.plot(ks, lam_r, 'r--', lw=2, label=r'$D_\parallel=0.5$')
  166.     ax.axhline(0.0, color='gray', lw=1)
  167.     ax.set_xlim(0.0, K_MAX)
  168.     ax.set_xlabel('k'); ax.set_ylabel(r'$\lambda(k)$')
  169.     ax.set_title(f'λ(k): P={P_val:.3f}, T={T_val:.3f}')
  170.     ax.legend(loc='upper right', frameon=False, fontsize=9)
  171.     plt.tight_layout(); plt.show()
  172.  
  173.  
  174. def lambda_max_point(P_val: float, T_val: float, Dpar: float) -> tuple[float, float]:
  175.     ks = np.linspace(0.0, K_MAX, K_POINTS)
  176.     lam = lambda_k(ks, P_val, T_val, Dpar=Dpar)
  177.     i = int(np.argmax(lam))
  178.     return float(ks[i]), float(lam[i])
  179.  
  180.  
  181. def plot_composite():
  182.  
  183.     P_grid = np.linspace(P_MIN, P_MAX, P_POINTS)
  184.     P_b, T_b = Tcrit_curve(P_grid, Dpar=0.1)
  185.     P_r, T_r = Tcrit_curve(P_grid, Dpar=0.5)
  186.  
  187.  
  188.     fig = plt.figure(figsize=(13.5, 8.4), dpi=120)
  189.     gs = fig.add_gridspec(3, 3, width_ratios=[1.3, 2.2, 1.3], wspace=0.35, hspace=0.35)
  190.  
  191.     axC = fig.add_subplot(gs[:, 1])
  192.     axC.plot(P_b, T_b, 'k-', lw=2.2, label=r'$D_\parallel=0.1$')
  193.     axC.plot(P_r, T_r, 'r--', lw=2.2, label=r'$D_\parallel=0.5$')
  194.     for (p, t) in CROSSES:
  195.         axC.plot(p, t, 'k*', ms=9)
  196.     axC.set_xlim(P_MIN, P_MAX); axC.set_ylim(0.0, T_MAX)
  197.     axC.set_xlabel('P'); axC.set_ylabel('T')
  198.     axC.legend(loc='upper right', frameon=False)
  199.     axC.set_title('Умови формування просторових структур')
  200.  
  201.  
  202.     left_points  = CROSSES[2::-1]   # (1.8, 1.3, 0.75) при P=0.1
  203.     right_points = CROSSES[5:2:-1]  # (1.8, 1.3, 0.75) при P=0.2
  204.  
  205.     def add_lambda_panel(ax, P_val, T_val):
  206.         ks = np.linspace(0.0, K_MAX, K_POINTS)
  207.         lam_b = lambda_k(ks, P_val, T_val, Dpar=0.1)
  208.         lam_r = lambda_k(ks, P_val, T_val, Dpar=0.5)
  209.         ax.plot(ks, lam_b, 'k-', lw=2, label=r'$D_\parallel=0.1$')
  210.         ax.plot(ks, lam_r, 'r--', lw=2, label=r'$D_\parallel=0.5$')
  211.         ax.axhline(0.0, color='gray', lw=1)
  212.         ax.set_xlim(0.0, K_MAX)
  213.         ax.set_title(f'P={P_val:.3f}, T={T_val:.3f}', fontsize=10)
  214.         ax.set_xlabel('k'); ax.set_ylabel(r'$\lambda(k)$')
  215.  
  216.         km, lm = lambda_max_point(P_val, T_val, Dpar=0.5)
  217.         ax.plot(km, lm, 'o', color='#2b6fff', ms=5)
  218.         return km, lm
  219.  
  220.     left_axes = [fig.add_subplot(gs[i, 0]) for i in range(3)]
  221.     for ax, (p, t) in zip(left_axes, left_points):
  222.         km, lm = add_lambda_panel(ax, p, t)
  223.         con = ConnectionPatch(xyA=(km, lm), coordsA=ax.transData,
  224.                               xyB=(p, t),   coordsB=axC.transData,
  225.                               arrowstyle='->', lw=1.6, color='#2b6fff')
  226.         fig.add_artist(con)
  227.  
  228.     right_axes = [fig.add_subplot(gs[i, 2]) for i in range(3)]
  229.     for ax, (p, t) in zip(right_axes, right_points):
  230.         km, lm = add_lambda_panel(ax, p, t)
  231.         con = ConnectionPatch(xyA=(km, lm), coordsA=ax.transData,
  232.                               xyB=(p, t),   coordsB=axC.transData,
  233.                               arrowstyle='->', lw=1.6, color='#2b6fff')
  234.         fig.add_artist(con)
  235.  
  236.     plt.show()
  237.     plt.close(fig)
  238.  
  239.  
  240.  
  241. if __name__ == '__main__':
  242.  
  243.     plot_PT_main()
  244.  
  245.     for i, (p, t) in enumerate(CROSSES, start=1):
  246.         plot_lambda_small(p, t)
  247.  
  248.     plot_composite()
  249.  
Advertisement
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!