Методика комп’ютерного експерименту ЛР№2

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.optimize import root, minimize


def equation(a, x, eps):
    return a * (1 - x) - x * np.exp(-eps * x)

def jacobian(a, x, eps):
    df_dx = -a - np.exp(-eps * x) * (1 - eps * x)
    return df_dx

def is_stable(stationary_points, stability):
    stable = []
    not_stable = []
    for i in range(len(stationary_points)):

        if stability[i] == 'b':
            stable.append(stationary_points[i])
            not_stable.append(np.nan)
        else:
            not_stable.append(stationary_points[i])
            stable.append(np.nan)
    stable = np.array(stable)
    not_stable = np.array(not_stable)

    return stable, not_stable


def find_min(x_values, alpha, epsilon):
    min_ = []
    for i in range(1, len(x_values)-2):
        d2V_dx2_prev = (V(x_values[i-1], alpha, epsilon) - 2*V(x_values[i], alpha, epsilon) + V(x_values[i+1], alpha, epsilon))/(dx*dx)
        d2V_dx2_next = (V(x_values[i], alpha, epsilon) - 2*V(x_values[i+1], alpha, epsilon) + V(x_values[i+2], alpha, epsilon))/(dx*dx)
        if d2V_dx2_prev * d2V_dx2_next < 0:
            min_.append(x_values[i])
    return min_


def V(x, alpha, epsilon):
    eq = -alpha*(x - 0.5*x**2) + (-x*epsilon*np.exp(-epsilon*x) - np.exp(-epsilon*x))/(epsilon**2)
    return eq


def alpha_st(x, epsilon):
    eq = -(x*np.exp(-epsilon*x))/(-1.0+x)
    return eq


epsilons = np.array([0.01, 0.1, 1, 2, 4, 6])
X = np.linspace(0, 1, 1000)


colors = [
    "#FF5733", "#33FF57", "#5733FF", "#FFFF33", "#33FFFF",
    "#FF33FF", "#FF5733", "#33FF57", "#5733FF", "#FFFF33",
    "#33FFFF", "#FF33FF", "#FF5733", "#33FF57", "#5733FF",
    "#FFFF33", "#33FFFF", "#FF33FF", "#FF5733", "#33FF57"
]

Stability = []
stationary_points = []
for eps in epsilons:
    stationary_points_eps = []
    Stability_eps = []
    for x in X:
        result = root(equation, x0=[0.1], args=(x, eps))
        if result.success:
            stationary_point = result.x
            stability = 'b' if jacobian(stationary_point, x, eps) < 0 else 'r'
            Stability_eps.append(stability)
            stationary_point = result.x
            stationary_points_eps.append(stationary_point[0])
        else:
            stationary_points_eps.append(np.nan)
            Stability_eps.append('white')
    stationary_points.append(stationary_points_eps)
    Stability.append(Stability_eps)

stationary_points = np.array(stationary_points)
Stability = np.array(Stability)

plt.figure(figsize=(15, 8))

for i, eps in enumerate(epsilons):
    plt.plot(stationary_points[i], X, color=colors[i], label=f'eps={np.round(eps, 2)}')


plt.xlabel(r'$\alpha$')
plt.ylabel('X')
plt.title(r'$\alpha (1 - x) x \exp(- \epsilon x) = 0$', size=20)
plt.legend()
plt.xlim(0, 0.5)
plt.grid(True)
plt.show()


plt.figure(figsize=(15, 8))
plt.plot(stationary_points[3], X, color=colors[3], lw=1.75, label=f'eps={np.round(epsilons[3], 2)}')


stable, not_stable = is_stable(stationary_points[3], Stability[3])
plt.plot(stable, X, color='b', lw=1, label=f'eps={np.round(epsilons[3], 2)}, стікий стан')
plt.plot(not_stable, X, color='red', lw=1, label=f'eps={np.round(epsilons[3], 2)}, не стійкий стан')
plt.xlabel(r'$\alpha$')
plt.ylabel('X')
plt.title(r'$\alpha (1 - x) x \exp(- \epsilon x) = 0$', size=20)
plt.legend()
plt.xlim(0, 0.5)
plt.grid(True)
plt.show()


plt.figure(figsize=(15, 8))
plt.plot(stationary_points[-1], X, color=colors[5], lw=1.75, label=f'eps={np.round(epsilons[-1], 2)}')
stable, not_stable = is_stable(stationary_points[-1], Stability[-1])
plt.plot(stable, X, color='b', lw=1, label=f'eps={np.round(epsilons[-1], 2)}, стікий стан')
plt.plot(not_stable, X, color='red', lw=1, label=f'eps={np.round(epsilons[-1], 2)}, не стійкий стан')
plt.plot()

plt.xlabel(r'$\alpha$')
plt.ylabel('X')
plt.title(r'$\alpha (1 - x) x \exp(- \epsilon x) = 0$', size=20)
plt.legend()
plt.xlim(0, 0.1)
plt.grid(True)
plt.show()





dx = 10**-4
x_values = np.arange(0, 1, dx)
alphas = [0.01, 0.025, 0.04, 0.06, 0.08, 0.09]

num_plots = len(alphas)
num_rows = 2
num_cols = 3

fig, axes = plt.subplots(num_rows, num_cols, figsize=(12, 10))

for i, alpha in enumerate(alphas):
    row = i // num_cols
    col = i % num_cols
    ax = axes[row, col]

    V_values = [V(x, alpha, 6) for x in x_values]
    ax.plot(x_values, V_values, label=f'alpha = {alpha}')
    ax.set_xlabel('$x$')
    ax.set_ylabel('$V(x)$')
    ax.set_title(f'$V(x): \epsilon=6$')
    ax.set_xlim(0, 1)
    ax.legend()
    ax.grid(True)

plt.tight_layout()
plt.show()











dx = 10**-4
x_values = np.arange(0, 1, dx)

min_per_epsilon = {}


for epsilon in np.arange(3, 6.1, 0.1):
    alpha_values = []
    for alpha in np.arange(0, 0.2, 0.01):
        V_values = [V(x, alpha, epsilon) for x in x_values]
        min_ = find_min(x_values, alpha, epsilon)
        if len(min_) == 2:
            V_min = [V(x, alpha, epsilon) for x in min_]
            if np.abs(V_min[0] - V_min[1]) < 0.001:
                alpha_values.append(alpha)


    min_per_epsilon[epsilon] = alpha_values


average_epsilon_per_alpha = {}
for epsilon, alpha_values in min_per_epsilon.items():
    for alpha in alpha_values:
        if alpha not in average_epsilon_per_alpha:
            average_epsilon_per_alpha[alpha] = []
        average_epsilon_per_alpha[alpha].append(epsilon)


for alpha, epsilons in average_epsilon_per_alpha.items():
    average_epsilon = np.mean(epsilons)






dx = 0.00001
epsilon_1_list = []
a2_1_list = []
epsilon_2_list = []
a2_2_list = []

for epsilon in np.arange(3, 6.1, 0.1):
    for x in np.arange(0.00001, 0.99999, dx):
        a1 = alpha_st(x - dx, epsilon)
        a2 = alpha_st(x, epsilon)
        a3 = alpha_st(x + dx, epsilon)
        if a1 < a2 and a3 < a2:
            epsilon_1_list.append(epsilon)
            a2_1_list.append(a2)
        if a1 > a2 and a3 > a2:
            epsilon_2_list.append(epsilon)
            a2_2_list.append(a2)
            break


alphas = []
average_epsilons = []

for alpha, epsilons in average_epsilon_per_alpha.items():
    alphas.append(alpha)
    average_epsilons.append(np.mean(epsilons))

plt.figure(figsize=(15, 8))
plt.plot(alphas, average_epsilons, marker='o', linestyle='-')
plt.plot(a2_1_list, epsilon_1_list, 'o', linestyle='-', label='Не стійкий стан')
plt.plot(a2_2_list, epsilon_2_list, 'o', linestyle='-', label='Стійкий стан')
plt.legend()
plt.xlabel(r'$\alpha$')
plt.ylabel(r'$\epsilon$')
plt.title('Фазова діаграма')
plt.grid(True)
plt.show()