Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import numpy as np
- import matplotlib.pyplot as plt
- from IPython.display import Latex
- from numpy.linalg import eigvals
- def system(x, z, s, r, b):
- dx_dt = s * x * (r - z - 1)
- dz_dt = x * x * r - x * x * z - b * z
- return dx_dt, dz_dt
- def runge_kutta(x, z, s, r, b, dt):
- k1x, k1z = system(x, z, s, r, b)
- k2x, k2z = system(x + 0.5 * k1x * dt, z + 0.5 * k1z * dt, s, r, b)
- k3x, k3z = system(x + 0.5 * k2x * dt, z + 0.5 * k2z * dt, s, r, b)
- k4x, k4z = system(x + k3x * dt, z + k3z * dt, s, r, b)
- x_new = x + (k1x + 2 * k2x + 2 * k3x + k4x) * dt / 6
- z_new = z + (k1z + 2 * k2z + 2 * k3z + k4z) * dt / 6
- return x_new, z_new
- def solve(s, r, b, x0, z0, dt, num_steps):
- x_values = [x0]
- z_values = [z0]
- for _ in range(num_steps):
- x, z = x_values[-1], z_values[-1]
- x_new, z_new = runge_kutta(x, z, s, r, b, dt)
- x_values.append(x_new)
- z_values.append(z_new)
- return x_values, z_values
- def jacobian(x, z, s, r, b):
- J = np.array([
- [s*(r-z-1), -s*x],
- [2*r*x - 2*x*z, x*x - b]
- ])
- return J
- def eigenvalues_at_point(x, z, s, r, b):
- Jacobian = jacobian(x, z, s, r, b)
- eig_vals = eigvals(Jacobian)
- return eig_vals
- def determine_point_type(eigvals):
- real_parts = np.real(eigvals)
- imag_parts = np.imag(eigvals)
- if all(rp < 0 for rp in real_parts) and all(ip == 0 for ip in imag_parts):
- return "Stable\ node"
- elif all(rp > 0 for rp in real_parts) and all(ip == 0 for ip in imag_parts):
- return "Unstable\ node"
- 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):
- return "Saddle\ point"
- elif all(rp == 0 for rp in real_parts) and all(ip != 0 for ip in imag_parts):
- return "Center"
- elif all(rp < 0 for rp in real_parts) and all(ip != 0 for ip in imag_parts):
- return "Stable\ Focus"
- elif all(rp > 0 for rp in real_parts) and all(ip != 0 for ip in imag_parts):
- return "Unstable\ Focus"
- else:
- return "Unknown"
- def plot_phase_portrait(ax, s, r, b, x_range=(-4, 4), z_range=(-1, 2)):
- x_vals = np.linspace(x_range[0], x_range[1], 100)
- z_vals = np.linspace(z_range[0], z_range[1], 100)
- X, Z = np.meshgrid(x_vals, z_vals)
- DX, DZ = system(X, Z, s, r, b)
- ax.streamplot(X, Z, DX, DZ, color='black')
- ax.set_xlabel(r'$x$')
- ax.set_ylabel(r'$z$')
- ax.set_title(fr'$Phase\ Portrait\ (r={r},\ b={b})$')
- x0 = 1
- z0 = 1
- s = 10
- r = 28
- b = 8/3
- dt = 0.01
- num_steps = 10000
- x_values, z_values = solve(s, r, b, x0, z0, dt, num_steps)
- s = 1
- b = 1
- points = [
- (0, 0),
- (np.sqrt(b * (r - 1)), r - 1),
- (-np.sqrt(b * (r - 1)), r - 1)
- ]
- r_values = np.linspace(0, 4, 100)
- plt.figure(figsize=(8, 5))
- plt.title(r"$Графік\ залежності\ x\ від\ r$")
- plt.xlabel(r"$r$")
- plt.ylabel(r"$x$")
- plt.grid()
- for point in points:
- x_values = []
- for r in r_values:
- x0, z0 = point
- dt = 0.01
- num_steps = 5000
- x, z = solve(1, r, 1, x0, z0, dt, num_steps)
- x_values.append(x[-1])
- plt.plot(r_values, x_values, label=fr'$x_0={round(point[0], 4)}, z_0={point[1]}$')
- plt.legend()
- plt.tight_layout()
- plt.show()
- r_values = np.linspace(0, 4, 100)
- plt.figure(figsize=(8, 5))
- plt.title(r"$Графік\ залежності\ z(r)$")
- plt.xlabel(r"$r$")
- plt.ylabel(r"$z$")
- plt.grid()
- for point in points:
- z_values = []
- for r in r_values:
- x0, z0 = point
- dt = 0.01
- num_steps = 5000
- x, z = solve(1, r, 1, x0, z0, dt, num_steps)
- z_values.append(z[-1])
- plt.plot(r_values, z_values, label=fr'$x_0={round(point[0], 4)}, z_0={point[1]}$')
- plt.legend()
- plt.tight_layout()
- plt.show()
- s = 1
- r = 0.5
- b = 1
- points = [
- (0, 0),
- (np.sqrt(b * (r - 1)), r - 1),
- (-np.sqrt(b * (r - 1)), r - 1)
- ]
- for point in points:
- x, z = point
- try:
- eig_vals = eigenvalues_at_point(x, z, s, r, b)
- point_type = determine_point_type(eig_vals)
- print()
- display(Latex(fr'$Point\ at\ x = {x},\ z = {z}\ is\ a\ {point_type}$'))
- for val in eig_vals:
- if np.iscomplex(val):
- display(Latex(fr'$Complex\ eigenvalue:\ {val}$'))
- real_part = np.real(val)
- imag_part = np.imag(val)
- if np.isclose(real_part, 0):
- display(Latex(fr"$Real\ part\ {real_part}\ is\ close\ to\ 0$"))
- else:
- display(Latex(fr"$Real\ part\ {real_part}\ is\ not\ close\ to\ 0$"))
- if np.sign(real_part) == 1:
- display(Latex(fr"$Real\ part\ {real_part}\ is\ positive$"))
- else:
- display(Latex(fr"$Real\ part\ {real_part}\ is\ negative$"))
- if np.isclose(imag_part, 0):
- display(Latex(fr"$Imaginary\ part\ {imag_part}\ is\ close\ to\ 0$"))
- else:
- display(Latex(fr"$Imaginary\ part\ {imag_part}\ is\ not\ close\ to\ 0$"))
- if np.sign(imag_part) == 1:
- display(Latex(fr"$Imaginary\ part\ {imag_part}\ is\ positive$"))
- else:
- display(Latex(fr"$Imaginary\ part\ {imag_part}\ is\ negative$"))
- else:
- display(Latex(fr"$Real\ eigenvalue:\ {val}$"))
- if val > 0:
- display(Latex(fr"$Real\ eigenvalue\ {val}\ is\ positive$"))
- else:
- display(Latex(fr"$Real\ eigenvalue\ {val}\ is\ negative$"))
- except:
- continue
- s = 1
- r = 1.1
- b = 1.5
- points = [
- (0, 0),
- (np.sqrt(b * (r - 1)), r - 1),
- (-np.sqrt(b * (r - 1)), r - 1)
- ]
- for point in points:
- x, z = point
- eig_vals = eigenvalues_at_point(x, z, s, r, b)
- point_type = determine_point_type(eig_vals)
- print()
- display(Latex(fr"$Point\ at\ x = {round(x, 4)},\ z = {round(z, 4)}\ is\ a\ {point_type}$"))
- for val in eig_vals:
- if np.iscomplex(val):
- display(Latex(fr"$Complex\ eigenvalue:\ {round(val, 4)}$"))
- real_part = np.real(val)
- imag_part = np.imag(val)
- if np.isclose(real_part, 0):
- display(Latex(fr"$Real\ part\ {round(real_part, 4)}\ is\ close\ to\ 0$"))
- else:
- display(Latex(fr"$Real\ part\ {round(real_part, 4)}\ is\ not\ close\ to\ 0$"))
- if np.sign(real_part) == 1:
- display(Latex(fr"$Real\ part\ {round(real_part, 4)}\ is\ positive$"))
- else:
- display(Latex(fr"$Real\ part\ {round(real_part, 4)}\ is\ negative$"))
- if np.isclose(imag_part, 0):
- display(Latex(fr"$Imaginary\ part\ {round(imag_part, 4)}\ is\ close\ to\ 0$"))
- else:
- display(Latex(fr"$Imaginary\ part\ {round(imag_part, 4)}\ is\ not\ close\ to\ 0$"))
- if np.sign(imag_part) == 1:
- display(Latex(fr"$Imaginary\ part\ {round(imag_part, 4)}\ is\ positive$"))
- else:
- display(Latex(fr"$Imaginary\ part\ {round(imag_part, 4)}\ is\ negative$"))
- else:
- display(Latex(fr"$Real\ eigenvalue:\ {round(val, 4)}$"))
- if val > 0:
- display(Latex(fr"$Real\ eigenvalue\ {round(val, 4)}\ is\ positive$"))
- else:
- display(Latex(fr"$Real\ eigenvalue\ {round(val, 4)}\ is\ negative$"))
- s = 1
- r = 1.5
- b = 1
- points = [
- (0, 0),
- (np.sqrt(b * (r - 1)), r - 1),
- (-np.sqrt(b * (r - 1)), r - 1)
- ]
- for point in points:
- x, z = point
- eig_vals = eigenvalues_at_point(x, z, s, r, b)
- point_type = determine_point_type(eig_vals)
- print()
- display(Latex(fr"$Point\ at\ x = {round(x, 4)},\ z = {round(z, 4)}\ is\ a\ {point_type}$"))
- for val in eig_vals:
- if np.iscomplex(val):
- display(Latex(fr"$Complex\ eigenvalue:\ {np.round(val, 4)}$"))
- real_part = np.real(val)
- imag_part = np.imag(val)
- if np.isclose(real_part, 0):
- display(Latex(fr"$Real\ part\ {round(real_part, 4)}\ is\ close\ to\ 0$"))
- else:
- display(Latex(fr"$Real\ part\ {round(real_part, 4)}\ is\ not\ close\ to\ 0$"))
- if np.sign(real_part) == 1:
- display(Latex(fr"$Real\ part\ {round(real_part, 4)}\ is\ positive$"))
- else:
- display(Latex(fr"$Real\ part\ {round(real_part, 4)}\ is\ negative$"))
- if np.isclose(imag_part, 0):
- display(Latex(fr"$Imaginary\ part\ {round(imag_part, 4)}\ is\ close\ to\ 0$"))
- else:
- display(Latex(fr"$Imaginary\ part\ {round(imag_part, 4)}\ is\ not\ close\ to\ 0$"))
- if np.sign(imag_part) == 1:
- display(Latex(fr"$Imaginary\ part\ {round(imag_part, 4)}\ is\ positive$"))
- else:
- display(Latex(fr"$Imaginary\ part\ {round(imag_part, 4)}\ is\ negative$"))
- else:
- display(Latex(fr"$Real\ eigenvalue:\ {round(val, 4)}$"))
- if val > 0:
- display(Latex(fr"$Real\ eigenvalue\ {round(val, 4)}\ is\ positive$"))
- else:
- display(Latex(fr"$Real\ eigenvalue\ {round(val, 4)}\ is\ negative$"))
- s = 1
- fig, axs = plt.subplots(2, 3, figsize=(15, 10))
- r = 0.9
- b = 0.7
- plot_phase_portrait(axs[0, 0], s, r, b)
- r = 1
- b = 1.5
- plot_phase_portrait(axs[0, 1], s, r, b)
- r = 1.1
- b = 1.5
- plot_phase_portrait(axs[0, 2], s, r, b)
- r = 1.2
- b = 1.8
- plot_phase_portrait(axs[1, 0], s, r, b)
- r = 1.5
- b = 1
- plot_phase_portrait(axs[1, 1], s, r, b)
- r = 2.5
- b = 1.5
- plot_phase_portrait(axs[1, 2], s, r, b)
- plt.tight_layout()
- plt.show()
- x = 1
- z = 1
- s = 1
- r_values = np.linspace(0, 10, 400)
- plt.figure(figsize=(10, 7))
- plt.fill_between(r_values, 0, (8 * (r_values - 1))/r_values**2, color='lightblue', label=r'$Сідло\ та\ 2\ стійких\ фокуса$')
- 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\ стійких\ вузли$')
- 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'$Один\ стійкий\ вузол$')
- r1 = 2
- b1 = 1.5
- plt.scatter(r1, b1, color='black', label=fr'$Point ({r1},\ {b1})$')
- r1 = 0.5
- b1 = 1
- plt.scatter(r1, b1, color='blue', label=fr'$Point ({r1},\ {b1})$')
- r1 = 1.2
- b1 = 1.75
- plt.scatter(r1, b1, color='red', label=fr'$Point ({r1},\ {b1})$')
- plt.xlabel(r'$r$')
- plt.ylabel(r'$b$')
- plt.title(r'$Dependence\ of\ b\ on\ r\ for\ the\ System$')
- plt.xlim(0, 3)
- plt.ylim(0, 2)
- plt.legend(loc=(0.69, 0))
- plt.grid(True)
- plt.show()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement