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

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.integrate import solve_ivp
from joblib import Parallel, delayed
from multiprocessing import current_process



def system(x, y, z, a, b, c):
    dx_dt = -y - z
    dy_dt = x + a * y
    dz_dt = b + z * (x - c)
    return dx_dt, dy_dt, dz_dt

def runge_kutta(x, y, z, a, b, c, dt):
    k1x, k1y, k1z = system(x, y, z, a, b, c)
    k2x, k2y, k2z = system(x + 0.5 * k1x * dt, y + 0.5 * k1y * dt, z + 0.5 * k1z * dt, a, b, c)
    k3x, k3y, k3z = system(x + 0.5 * k2x * dt, y + 0.5 * k2y * dt, z + 0.5 * k2z * dt, a, b, c)
    k4x, k4y, k4z = system(x + k3x * dt, y + k3y * dt, z + k3z * dt, a, b, c)
    x_new = x + (k1x + 2 * k2x + 2 * k3x + k4x) * dt / 6
    y_new = y + (k1y + 2 * k2y + 2 * k3y + k4y) * dt / 6
    z_new = z + (k1z + 2 * k2z + 2 * k3z + k4z) * dt / 6
    return x_new, y_new, z_new

def solve(a, b, c, x0, y0, z0, dt, num_steps):
    x_values = [x0]
    y_values = [y0]
    z_values = [z0]
    for _ in range(num_steps):
        x, y, z = x_values[-1], y_values[-1], z_values[-1]
        x_new, y_new, z_new = runge_kutta(x, y, z, a, b, c, dt)
        x_values.append(x_new)
        y_values.append(y_new)
        z_values.append(z_new)
    return x_values, y_values, z_values

x0 = -1
y0 = -1
z0 = 1

a = 0.2
b = 0.2
c = 14

dt = 0.01
num_steps = 100000



# Розв'язання системи Реслера
x_values, y_values, z_values = solve(a, b, c, x0, y0, z0, dt, num_steps)

# Побудова тривимірного графіка
fig = plt.figure(figsize=(12, 10))
ax = fig.add_subplot(221, projection='3d')
ax.plot(x_values, y_values, z_values, '-', linewidth=0.5)
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$y$')
ax.set_zlabel(r'$z$')
ax.set_title(r"$Фазовий\ портрет\ системи\ Реслера$")

# Графіки x(t), y(t), z(t)
ax_x = fig.add_subplot(222)
ax_x.plot(np.arange(len(x_values)) * dt, x_values, '-', linewidth=0.5)
ax_x.set_xlabel(r'$t$')
ax_x.set_ylabel(r'$x$')
ax_x.set_title(r"$Залежність\ x(t)$")

ax_y = fig.add_subplot(223)
ax_y.plot(np.arange(len(y_values)) * dt, y_values, '-', linewidth=0.5)
ax_y.set_xlabel(r'$t$')
ax_y.set_ylabel(r'$y$')
ax_y.set_title(r"$Залежність\ y(t)$")

ax_z = fig.add_subplot(224)
ax_z.plot(np.arange(len(z_values)) * dt, z_values, '-', linewidth=0.5)
ax_z.set_xlabel(r'$t$')
ax_z.set_ylabel(r'$z$')
ax_z.set_title(r"$Залежність\ z(t)$")

plt.tight_layout()
plt.show()


def ressler(t, y, a, b, c):
    dydt = [
        -y[1] - y[2],
        y[0] + a * y[1],
        b + y[2] * (y[0] - c)
    ]
    return dydt

def lyapunov(a, b, c, p=0.5, duration=200, delta_t=0.01, num_iterations=10):
    np.random.seed(0)  # Для відтворюваності результатів
    y0 = np.random.rand(3) * 20 - 10  # Випадкові початкові умови

    L_sum = 0
    t_span = [0, duration]
    t_eval = np.arange(0, duration, delta_t)

    for _ in range(num_iterations):
        y = y0
        L = 0

        sol = solve_ivp(ressler, t_span, y, args=(a, b, c), t_eval=t_eval)

        for i in range(len(sol.t) - 1):
            t0, t1 = sol.t[i], sol.t[i + 1]
            y0, y1 = sol.y[:, i], sol.y[:, i + 1]

            # Розв'язуємо систему з невеликим зсувом
            y0_perturbed = y0 + p * np.random.rand(3)
            sol_perturbed = solve_ivp(ressler, [t0, t1], y0_perturbed, args=(a, b, c))
            y1_perturbed = sol_perturbed.y[:, -1]

            # Обчислення ляпуновського показника
            delta = np.linalg.norm(y1 - y1_perturbed)
            L += np.log(delta / p)

        L_sum += L / duration

    return L_sum / num_iterations

a = 0.3
b = 0.3
c = 28

lyap_max = lyapunov(a, b, c)
display(Latex((fr"$Старший\ показник\ Ляпунова\ для\ системи\ Ресслера:\ {lyap_max}$"))












a_range = np.linspace(0.2, 0.3, 100)
c_range = np.linspace(20, 45, 100)

lyap_values = np.zeros((len(a_range), len(c_range)))

start_time = time.time()

idx, values = [], []
for i, a in enumerate(a_range):
    for j, c in enumerate(c_range):
        idx.append((i, j))
        values.append([a, a, c])

results = Parallel(n_jobs=-1)(
            delayed(lyapunov)(
                values[i][0], values[i][1], values[i][2], delta_t=5, index=i) for i in range(len(idx)))

for i in range(len(results)):
    lyap_values[idx[results[i][0]]] = results[i][1]


end_time = time.time()
my_time = end_time - start_time

plt.figure(figsize=(10, 6))
plt.imshow(lyap_values, extent=[c_range[0], c_range[-1], a_range[0], a_range[-1]], aspect='auto', cmap='jet')
plt.colorbar(label=r'$Старший\ показник\ Ляпунова$')
plt.xlabel(r'$Параметр\ c$')
plt.ylabel(r'$Параметр\ a$')
plt.title(r'$Карта\ старших\ показників\ Ляпунова\ для\ системи\ Ресслера$')
plt.show()

display(Latex(fr'$Час\ моделювання: \quad {my_time_(round(my_time))}.$'))












a_values = [0.2, 0.2, 0.24, 0.29, 0.26, 0.3]
b_values = [0.2, 0.2, 0.24, 0.29, 0.26, 0.3]
c_values = [30, 43, 35, 32, 26, 27]

dt = 0.01
num_steps = 5000

fig = plt.figure(figsize=(18, 15))

for i in range(len(a_values)):
    a = a_values[i]
    b = b_values[i]
    c = c_values[i]

    x0 = -1
    y0 = -1
    z0 = 1

    x_values, y_values, z_values = solve(a, b, c, x0, y0, z0, dt, num_steps)

    ax = fig.add_subplot(3, 3, i + 1, projection='3d')
    ax.plot(x_values, y_values, z_values, 'r-', linewidth=0.5)
    ax.set_xlabel(r'$x$')
    ax.set_ylabel(r'$y$')
    ax.set_zlabel(r'$z$')
    ax.set_title(fr"$a={a},\ b={b},\ c={c}$")

plt.tight_layout()
plt.show()









# Знаходження комбінацій параметрів з найвищим, найнижчим і середнім значенням ляпуновського показника
max_idx = np.unravel_index(np.argmax(lyap_values), lyap_values.shape)
min_idx = np.unravel_index(np.argmin(lyap_values), lyap_values.shape)
mean_idx = np.unravel_index(np.nanargmin(np.abs(lyap_values - np.nanmean(lyap_values))), lyap_values.shape)

max_a, max_c = a_range[max_idx[0]], c_range[max_idx[1]]
min_a, min_c = a_range[min_idx[0]], c_range[min_idx[1]]
mean_a, mean_c = a_range[mean_idx[0]], c_range[mean_idx[1]]


display(Latex(fr'$Combination\ with\ highest\ Lyapunov\ exponent:\ {max_a}\quad {max_c}\quad {np.max(lyap_values)}$'))
display(Latex(fr"$Combination\ with\ lowest\ Lyapunov\ exponent:\ {min_a}\quad {min_c}\quad {np.min(lyap_values)}$"))
display(Latex(fr"$Combination\ with\ mean\ Lyapunov\ exponent:\ {mean_a}\quad {mean_c}\quad {np.mean(lyap_values)}$"))


def plot_3d_phase_portrait(a, b, c, ax):
    y0 = np.random.rand(3) * 20 - 10  # Випадкові початкові умови
    sol = solve_ivp(ressler, [0, 200], y0, args=(a, b, c), dense_output=True)
    t = np.linspace(0, 200, 10000)
    y = sol.sol(t)
    ax.plot(y[0], y[1], y[2], label=fr'$a={round(a, 6)},\ c={round(c, 6)}$')
    ax.set_xlabel(fr'$x$')
    ax.set_ylabel(fr'$y$')
    ax.set_zlabel(fr'$z$')
    ax.legend()


# Побудова 3D фазових портретів для комбінацій параметрів з найвищим, найнижчим і середнім значенням ляпуновського показника
fig = plt.figure(figsize=(15, 5))

ax1 = fig.add_subplot(131, projection='3d')
plot_3d_phase_portrait(max_a, max_a, max_c, ax1)
ax1.set_title(fr'$Highest\ Lyapunov\ Exponent$')

ax2 = fig.add_subplot(132, projection='3d')
plot_3d_phase_portrait(min_a, min_a, min_c, ax2)
ax2.set_title(fr'$Lowest\ Lyapunov\ Exponent$')

ax3 = fig.add_subplot(133, projection='3d')
plot_3d_phase_portrait(mean_a, mean_a, mean_c, ax3)
ax3.set_title(fr'$Mean\ Lyapunov\ Exponent$')

plt.tight_layout()
plt.show()