Advertisement
mirosh111000

Мірошниченко_Артем_ПМ-11_Лр№5_Мет_комп_експ

May 4th, 2024
33
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 11.23 KB | None | 0 0
  1. import numpy as np
  2. import matplotlib.pyplot as plt
  3. from IPython.display import Latex
  4. from numpy.linalg import eigvals
  5.  
  6. def system(x, z, s, r, b):
  7.     dx_dt = s * x * (r - z - 1)
  8.     dz_dt = x * x * r - x * x * z - b * z
  9.     return dx_dt, dz_dt
  10.  
  11. def runge_kutta(x, z, s, r, b, dt):
  12.     k1x, k1z = system(x, z, s, r, b)
  13.     k2x, k2z = system(x + 0.5 * k1x * dt,  z + 0.5 * k1z * dt, s, r, b)
  14.     k3x, k3z = system(x + 0.5 * k2x * dt,  z + 0.5 * k2z * dt, s, r, b)
  15.     k4x, k4z = system(x + k3x * dt,  z + k3z * dt, s, r, b)
  16.     x_new = x + (k1x + 2 * k2x + 2 * k3x + k4x) * dt / 6
  17.     z_new = z + (k1z + 2 * k2z + 2 * k3z + k4z) * dt / 6
  18.     return x_new, z_new
  19.  
  20. def solve(s, r, b, x0, z0, dt, num_steps):
  21.     x_values = [x0]
  22.     z_values = [z0]
  23.     for _ in range(num_steps):
  24.         x, z = x_values[-1], z_values[-1]
  25.         x_new,  z_new = runge_kutta(x, z, s, r, b, dt)
  26.         x_values.append(x_new)
  27.         z_values.append(z_new)
  28.     return x_values,  z_values
  29.  
  30. def jacobian(x, z, s, r, b):
  31.     J = np.array([
  32.         [s*(r-z-1), -s*x],
  33.         [2*r*x - 2*x*z, x*x - b]
  34.     ])
  35.     return J
  36.  
  37. def eigenvalues_at_point(x, z, s, r, b):
  38.     Jacobian = jacobian(x, z, s, r, b)
  39.     eig_vals = eigvals(Jacobian)
  40.     return eig_vals
  41.  
  42. def determine_point_type(eigvals):
  43.     real_parts = np.real(eigvals)
  44.     imag_parts = np.imag(eigvals)
  45.     if all(rp < 0 for rp in real_parts) and all(ip == 0 for ip in imag_parts):
  46.         return "Stable\ node"
  47.     elif all(rp > 0 for rp in real_parts) and all(ip == 0 for ip in imag_parts):
  48.         return "Unstable\ node"
  49.     elif any(rp > 0 for rp in real_parts) and any(rp < 0 for rp in real_parts) and all(ip == 0 for ip in imag_parts):
  50.         return "Saddle\ point"
  51.     elif all(rp == 0 for rp in real_parts) and all(ip != 0 for ip in imag_parts):
  52.         return "Center"
  53.     elif  all(rp < 0 for rp in real_parts) and all(ip != 0 for ip in imag_parts):
  54.         return "Stable\ Focus"
  55.     elif  all(rp > 0 for rp in real_parts) and all(ip != 0 for ip in imag_parts):
  56.         return "Unstable\ Focus"
  57.     else:
  58.         return "Unknown"
  59.    
  60.    
  61. def plot_phase_portrait(ax, s, r, b, x_range=(-4, 4), z_range=(-1, 2)):
  62.     x_vals = np.linspace(x_range[0], x_range[1], 100)
  63.     z_vals = np.linspace(z_range[0], z_range[1], 100)
  64.     X, Z = np.meshgrid(x_vals, z_vals)
  65.     DX, DZ = system(X, Z, s, r, b)
  66.  
  67.     ax.streamplot(X, Z, DX, DZ, color='black')
  68.     ax.set_xlabel(r'$x$')
  69.     ax.set_ylabel(r'$z$')
  70.     ax.set_title(fr'$Phase\ Portrait\ (r={r},\ b={b})$')
  71.  
  72.  
  73.  
  74. x0 = 1
  75. z0 = 1
  76.  
  77. s = 10
  78. r = 28
  79. b = 8/3
  80.  
  81. dt = 0.01
  82. num_steps = 10000
  83. x_values, z_values = solve(s, r, b, x0,  z0, dt, num_steps)
  84. s = 1
  85. b = 1
  86.  
  87. points = [
  88.     (0, 0),
  89.     (np.sqrt(b * (r - 1)), r - 1),
  90.     (-np.sqrt(b * (r - 1)), r - 1)
  91. ]
  92.  
  93. r_values = np.linspace(0, 4, 100)
  94. plt.figure(figsize=(8, 5))
  95. plt.title(r"$Графік\ залежності\ x\ від\ r$")
  96. plt.xlabel(r"$r$")
  97. plt.ylabel(r"$x$")
  98. plt.grid()
  99.  
  100. for point in points:
  101.     x_values = []
  102.     for r in r_values:
  103.         x0,  z0 = point
  104.         dt = 0.01
  105.         num_steps = 5000
  106.         x, z = solve(1, r, 1, x0, z0, dt, num_steps)
  107.         x_values.append(x[-1])
  108.     plt.plot(r_values, x_values, label=fr'$x_0={round(point[0], 4)}, z_0={point[1]}$')
  109.    
  110. plt.legend()
  111. plt.tight_layout()
  112. plt.show()
  113.  
  114. r_values = np.linspace(0, 4, 100)
  115. plt.figure(figsize=(8, 5))
  116. plt.title(r"$Графік\ залежності\ z(r)$")
  117. plt.xlabel(r"$r$")
  118. plt.ylabel(r"$z$")
  119. plt.grid()
  120.  
  121. for point in points:
  122.     z_values = []
  123.     for r in r_values:
  124.         x0, z0 = point
  125.         dt = 0.01
  126.         num_steps = 5000
  127.         x, z = solve(1, r, 1, x0, z0, dt, num_steps)
  128.         z_values.append(z[-1])
  129.     plt.plot(r_values, z_values, label=fr'$x_0={round(point[0], 4)}, z_0={point[1]}$')
  130.  
  131. plt.legend()
  132. plt.tight_layout()
  133. plt.show()
  134.  
  135.  
  136.  
  137.  
  138.  
  139.  
  140.  
  141.  
  142.  
  143. s = 1
  144. r = 0.5
  145. b = 1
  146.  
  147. points = [
  148.     (0, 0),
  149.     (np.sqrt(b * (r - 1)), r - 1),
  150.     (-np.sqrt(b * (r - 1)), r - 1)
  151. ]
  152.  
  153. for point in points:
  154.     x, z = point
  155.     try:
  156.         eig_vals = eigenvalues_at_point(x, z, s, r, b)
  157.         point_type = determine_point_type(eig_vals)
  158.         print()
  159.         display(Latex(fr'$Point\ at\ x = {x},\ z = {z}\ is\ a\ {point_type}$'))
  160.  
  161.         for val in eig_vals:
  162.             if np.iscomplex(val):
  163.                 display(Latex(fr'$Complex\ eigenvalue:\ {val}$'))
  164.                 real_part = np.real(val)
  165.                 imag_part = np.imag(val)
  166.                 if np.isclose(real_part, 0):
  167.                     display(Latex(fr"$Real\ part\ {real_part}\ is\ close\ to\ 0$"))
  168.                 else:
  169.                     display(Latex(fr"$Real\ part\ {real_part}\ is\ not\ close\ to\ 0$"))
  170.  
  171.                 if np.sign(real_part) == 1:
  172.                     display(Latex(fr"$Real\ part\ {real_part}\ is\ positive$"))
  173.                 else:
  174.                     display(Latex(fr"$Real\ part\ {real_part}\ is\ negative$"))
  175.  
  176.                 if np.isclose(imag_part, 0):
  177.                     display(Latex(fr"$Imaginary\ part\ {imag_part}\ is\ close\ to\ 0$"))
  178.                 else:
  179.                     display(Latex(fr"$Imaginary\ part\ {imag_part}\ is\ not\ close\ to\ 0$"))
  180.  
  181.                 if np.sign(imag_part) == 1:
  182.                     display(Latex(fr"$Imaginary\ part\ {imag_part}\ is\ positive$"))
  183.                 else:
  184.                     display(Latex(fr"$Imaginary\ part\ {imag_part}\ is\ negative$"))
  185.             else:
  186.                 display(Latex(fr"$Real\ eigenvalue:\ {val}$"))
  187.                 if val > 0:
  188.                     display(Latex(fr"$Real\ eigenvalue\ {val}\ is\ positive$"))
  189.                 else:
  190.                     display(Latex(fr"$Real\ eigenvalue\ {val}\ is\ negative$"))
  191.     except:
  192.         continue
  193.        
  194.        
  195.        
  196.        
  197.        
  198.  
  199.        
  200.        
  201.        
  202. s = 1
  203. r = 1.1
  204. b = 1.5
  205.  
  206. points = [
  207.     (0, 0),
  208.     (np.sqrt(b * (r - 1)), r - 1),
  209.     (-np.sqrt(b * (r - 1)), r - 1)
  210. ]
  211.  
  212. for point in points:
  213.     x, z = point
  214.     eig_vals = eigenvalues_at_point(x, z, s, r, b)
  215.     point_type = determine_point_type(eig_vals)
  216.     print()
  217.     display(Latex(fr"$Point\ at\ x = {round(x, 4)},\ z = {round(z, 4)}\ is\ a\ {point_type}$"))
  218.     for val in eig_vals:
  219.         if np.iscomplex(val):
  220.             display(Latex(fr"$Complex\ eigenvalue:\ {round(val, 4)}$"))
  221.             real_part = np.real(val)
  222.             imag_part = np.imag(val)
  223.             if np.isclose(real_part, 0):
  224.                 display(Latex(fr"$Real\ part\ {round(real_part, 4)}\ is\ close\ to\ 0$"))
  225.             else:
  226.                 display(Latex(fr"$Real\ part\ {round(real_part, 4)}\ is\ not\ close\ to\ 0$"))
  227.  
  228.             if np.sign(real_part) == 1:
  229.                 display(Latex(fr"$Real\ part\ {round(real_part, 4)}\ is\ positive$"))
  230.             else:
  231.                 display(Latex(fr"$Real\ part\ {round(real_part, 4)}\ is\ negative$"))
  232.  
  233.             if np.isclose(imag_part, 0):
  234.                 display(Latex(fr"$Imaginary\ part\ {round(imag_part, 4)}\ is\ close\ to\ 0$"))
  235.             else:
  236.                 display(Latex(fr"$Imaginary\ part\ {round(imag_part, 4)}\ is\ not\ close\ to\ 0$"))
  237.  
  238.             if np.sign(imag_part) == 1:
  239.                 display(Latex(fr"$Imaginary\ part\ {round(imag_part, 4)}\ is\ positive$"))
  240.             else:
  241.                 display(Latex(fr"$Imaginary\ part\ {round(imag_part, 4)}\ is\ negative$"))
  242.         else:
  243.             display(Latex(fr"$Real\ eigenvalue:\ {round(val, 4)}$"))
  244.             if val > 0:
  245.                 display(Latex(fr"$Real\ eigenvalue\ {round(val, 4)}\ is\ positive$"))
  246.             else:
  247.                 display(Latex(fr"$Real\ eigenvalue\ {round(val, 4)}\ is\ negative$"))
  248.                
  249.                
  250.                
  251.  
  252.                
  253.                
  254.                
  255.                
  256.  
  257. s = 1
  258. r = 1.5
  259. b = 1
  260.  
  261.  
  262. points = [
  263.     (0, 0),
  264.     (np.sqrt(b * (r - 1)), r - 1),
  265.     (-np.sqrt(b * (r - 1)), r - 1)
  266. ]              
  267.                
  268. for point in points:
  269.     x, z = point
  270.     eig_vals = eigenvalues_at_point(x, z, s, r, b)
  271.     point_type = determine_point_type(eig_vals)
  272.     print()
  273.     display(Latex(fr"$Point\ at\ x = {round(x, 4)},\ z = {round(z, 4)}\ is\ a\ {point_type}$"))
  274.     for val in eig_vals:
  275.         if np.iscomplex(val):
  276.             display(Latex(fr"$Complex\ eigenvalue:\ {np.round(val, 4)}$"))
  277.             real_part = np.real(val)
  278.             imag_part = np.imag(val)
  279.             if np.isclose(real_part, 0):
  280.                 display(Latex(fr"$Real\ part\ {round(real_part, 4)}\ is\ close\ to\ 0$"))
  281.             else:
  282.                 display(Latex(fr"$Real\ part\ {round(real_part, 4)}\ is\ not\ close\ to\ 0$"))
  283.  
  284.             if np.sign(real_part) == 1:
  285.                 display(Latex(fr"$Real\ part\ {round(real_part, 4)}\ is\ positive$"))
  286.             else:
  287.                 display(Latex(fr"$Real\ part\ {round(real_part, 4)}\ is\ negative$"))
  288.  
  289.             if np.isclose(imag_part, 0):
  290.                 display(Latex(fr"$Imaginary\ part\ {round(imag_part, 4)}\ is\ close\ to\ 0$"))
  291.             else:
  292.                 display(Latex(fr"$Imaginary\ part\ {round(imag_part, 4)}\ is\ not\ close\ to\ 0$"))
  293.  
  294.             if np.sign(imag_part) == 1:
  295.                 display(Latex(fr"$Imaginary\ part\ {round(imag_part, 4)}\ is\ positive$"))
  296.             else:
  297.                 display(Latex(fr"$Imaginary\ part\ {round(imag_part, 4)}\ is\ negative$"))
  298.         else:
  299.             display(Latex(fr"$Real\ eigenvalue:\ {round(val, 4)}$"))
  300.             if val > 0:
  301.                 display(Latex(fr"$Real\ eigenvalue\ {round(val, 4)}\ is\ positive$"))
  302.             else:
  303.                 display(Latex(fr"$Real\ eigenvalue\ {round(val, 4)}\ is\ negative$"))
  304.                
  305.                
  306.                
  307.                
  308.                
  309.                
  310.                
  311.                
  312.                
  313. s = 1
  314.  
  315. fig, axs = plt.subplots(2, 3, figsize=(15, 10))
  316.  
  317. r = 0.9
  318. b = 0.7
  319. plot_phase_portrait(axs[0, 0], s, r, b)
  320.  
  321. r = 1
  322. b = 1.5
  323. plot_phase_portrait(axs[0, 1], s, r, b)
  324.  
  325. r = 1.1
  326. b = 1.5
  327. plot_phase_portrait(axs[0, 2], s, r, b)
  328.  
  329. r = 1.2
  330. b = 1.8
  331. plot_phase_portrait(axs[1, 0], s, r, b)
  332.  
  333. r = 1.5
  334. b = 1
  335. plot_phase_portrait(axs[1, 1], s, r, b)
  336.  
  337. r = 2.5
  338. b = 1.5
  339. plot_phase_portrait(axs[1, 2], s, r, b)
  340. plt.tight_layout()
  341. plt.show()
  342.  
  343.  
  344.  
  345.  
  346.  
  347.  
  348.  
  349.  
  350.  
  351. x = 1
  352. z = 1
  353. s = 1
  354.  
  355. r_values = np.linspace(0, 10, 400)
  356. plt.figure(figsize=(10, 7))
  357. plt.fill_between(r_values, 0, (8 * (r_values - 1))/r_values**2, color='lightblue', label=r'$Сідло\ та\ 2\ стійких\ фокуса$')
  358. plt.fill_between(r_values, (8 * (r_values - 1))/r_values**2, np.max((8 * (r_values - 1))/r_values**2), color='lightgreen', label=r'$Сідло\ та\ 2\ стійких\ вузли$')
  359. plt.fill_between(r_values, np.max((8 * (r_values - 1))/r_values**2), 0, where=(r_values < 1), color='red', alpha=0.3, label=r'$Один\ стійкий\ вузол$')
  360. r1 = 2
  361. b1 = 1.5
  362. plt.scatter(r1, b1, color='black', label=fr'$Point ({r1},\ {b1})$')
  363.  
  364. r1 = 0.5
  365. b1 = 1
  366. plt.scatter(r1, b1, color='blue', label=fr'$Point ({r1},\ {b1})$')
  367.  
  368.  
  369. r1 = 1.2
  370. b1 = 1.75
  371. plt.scatter(r1, b1, color='red', label=fr'$Point ({r1},\ {b1})$')
  372.  
  373.  
  374. plt.xlabel(r'$r$')
  375. plt.ylabel(r'$b$')
  376. plt.title(r'$Dependence\ of\ b\ on\ r\ for\ the\ System$')
  377. plt.xlim(0, 3)
  378. plt.ylim(0, 2)
  379. plt.legend(loc=(0.69, 0))
  380. plt.grid(True)
  381. plt.show()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement