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

import numpy as np
import matplotlib.pyplot as plt
import time
from IPython.display import Latex



def my_time_(s):
    h = s // 3600
    time = s % 3600
    m = time // 60
    s = time % 60
    if h < 10:
        h = '0' + str(int(h))
    if m < 10:
        m = '0' + str(int(m))
    if s < 10:
        s = '0' + str(int(s))

    return fr'{h}:{m}:{s}'


def phi(r, sigma=SIGMA, epsilon=EPSILON):
    sigma_r_pow_6 = (sigma / r) ** 6
    return 4 * epsilon * (sigma_r_pow_6 ** 2 - sigma_r_pow_6)
def init_state(number_of_pieces=NUMBER_OF_PIECES, mass=1):
    pieces = np.zeros((number_of_pieces, 5))
    pieces[:, MASS] = np.full(number_of_pieces, mass)  # mass
    pieces[:, X] = np.random.rand(number_of_pieces) * BOX_WIDTH  # x
    pieces[:, Y] = np.random.rand(number_of_pieces) * BOX_HEIGHT  # y
    pieces[:, VX] = (2 * np.random.rand(number_of_pieces) - 1) * V_MAX  # vx
    pieces[:, VY] = (2 * np.random.rand(number_of_pieces) - 1) * V_MAX  # vy

    pieces[:, VX] -= np.mean(pieces[:, VX])
    pieces[:, VY] -= np.mean(pieces[:, VY])

    while True:
        flag = False
        for i in range(0, number_of_pieces):
            for j in range(i + 1, number_of_pieces):
                rx = pieces[i, X] - pieces[j, X]
                ry = pieces[i, Y] - pieces[j, Y]
                r = (rx ** 2 + ry ** 2) ** 0.5
                if r < 1:
                    pieces[j, X] = np.random.rand(1)[0] * BOX_WIDTH
                    flag = True
        if not flag:
            break

    return pieces
class FCalc:
    def __init__(self, number_of_pieces=NUMBER_OF_PIECES, epsilon=EPSILON, sigma=SIGMA):
        self.number_of_pieces = number_of_pieces
        self.epsilon = epsilon
        self.sigma = sigma
        self.cache = np.full((number_of_pieces, number_of_pieces), None)

    def get_from_cache(self, i, j):
        return self.cache[i][j]

    def save_to_cache(self, Fr, i, j):
        self.cache[i][j] = Fr


    def f(self, r):
        sigma_r = self.sigma / r
        sigma_r_pow_6 = sigma_r ** 6
        return 24 * self.epsilon * sigma_r * sigma_r_pow_6 * (2 * sigma_r_pow_6 - 1)

    def fi(self, i, state, axis):
        forces = np.zeros(self.number_of_pieces)
        for j in range(0, self.number_of_pieces):
            if i == j:
                Fr = 0
            else:
                x = state[i, X] - state[j, X]
                y = state[i, Y] - state[j, Y]
                r = (y ** 2 + x ** 2) ** 0.5

                Fr_from_cache = self.get_from_cache(i, j)
                if Fr_from_cache is not None:
                    Fr = Fr_from_cache
                else:
                    Fr = self.f(r)
                    self.save_to_cache(Fr, i, j)

                if axis == X:
                    Fr = x / r * Fr
                elif axis == Y:
                    Fr = y / r * Fr

            forces[j] = Fr

        return np.sum(forces)


def mean_kinetic_energy(state):
    mass = state[:, MASS]
    vx = state[:, VX]
    vy = state[:, VY]
    v2 = vx ** 2 + vy ** 2
    kinetic_energy = 0.5 * mass * v2
    return np.mean(kinetic_energy)


def rt(t0_state, t_state):
    xt0 = t0_state[:, X]
    yt0 = t0_state[:, Y]
    xt = t_state[:, X]
    yt = t_state[:, Y]
    return np.hypot((xt - xt0), (yt - yt0))


def mean_potential_energy(state):
    number_of_pieces = state.shape[0]

    counter = 0
    phis_count = int(number_of_pieces * (number_of_pieces + 1) / 2)
    phis = np.zeros(phis_count)
    for i in range(0, number_of_pieces):
        for j in range(i + 1, number_of_pieces):
            x = state[i, X] - state[j, X]
            y = state[i, Y] - state[j, Y]
            r = (y ** 2 + x ** 2) ** 0.5
            phis[counter] = phi(abs(r))
            counter += 1

    return np.sum(phis) / number_of_pieces


def Rt2(state, t0_state):
    return np.mean(rt(state, t0_state) ** 2)


def handle_corner_cases(x, y, vx, vy):
    if x < 0:
        x = -x
        vx = -vx
    elif x > BOX_WIDTH:
        x = BOX_WIDTH - (x - BOX_WIDTH)
        vx = -vx
    if y < 0:
        y = -y
        vy = -vy
    elif y > BOX_HEIGHT:
        y = BOX_HEIGHT - (y - BOX_HEIGHT)
        vy = -vy
    return x, y, vx, vy


# Ейлер
def calc_next_state(state, epsilon=EPSILON, sigma=SIGMA):
    next_state = np.copy(state)
    f_calc = FCalc(NUMBER_OF_PIECES, epsilon, sigma)

    for i, piece_state in enumerate(state):
        Fix = f_calc.fi(i, state, X)
        Fiy = f_calc.fi(i, state, Y)
        axi = Fix / piece_state[MASS]
        ayi = Fiy / piece_state[MASS]
        vxi = piece_state[VX] + axi * DELTA_T
        vyi = piece_state[VY] + ayi * DELTA_T

        xi = piece_state[X] + vxi * DELTA_T
        yi = piece_state[Y] + vyi * DELTA_T
        xi, yi, vxi, vyi = handle_corner_cases(xi, yi, vxi, vyi)

        next_state[i][X] = xi
        next_state[i][Y] = yi
        next_state[i][VX] = vxi
        next_state[i][VY] = vyi

    return next_state





NUMBER_OF_PIECES = 100
EPSILON = 1
SIGMA = 1
BOX_WIDTH = 20
BOX_HEIGHT = 20
WEGHT = 1
DELTA_T = 0.001
V_MAX = 1

MASS, X, Y, VX, VY = [0, 1, 2, 3, 4]


start_time = time.time()

initial_state = init_state()
x_values = [initial_state[-1, X]]
y_values = [initial_state[-1, Y]]

current_state = initial_state
current_time = 0
total_time = 100

while current_time < total_time:
    current_state = calc_next_state(current_state)
    print(f'current_time={current_time}\n')
    x_values.append(current_state[-1, X])
    y_values.append(current_state[-1, Y])
    current_time += DELTA_T

plt.plot(x_values, y_values)
plt.scatter(x_values[0], y_values[0], c='red', label=r'$start$')
plt.scatter(x_values[-1], y_values[-1], c='blue', label=r'$finish$')
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')
plt.title(r'$Фазовий\ портрет\ частинки$')
plt.legend()
plt.grid(True)
plt.show()

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

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








start_time = time.time()

initial_state = init_state()
current_state = initial_state
kinetic_energies = []
potential_energies = []
total_energies = []

plt.figure(figsize=(10, 6))

current_time = 0
while current_time < total_time:
    print(f'current_time={current_time}\n')
    kinetic_energy_value = mean_kinetic_energy(current_state)
    potential_energy_value = mean_potential_energy(current_state)
    total_energy_value = kinetic_energy_value + potential_energy_value
    kinetic_energies.append(kinetic_energy_value)
    potential_energies.append(potential_energy_value)
    total_energies.append(total_energy_value)
    current_state = calc_next_state(current_state)
    current_time += DELTA_T


time_points = np.linspace(0, 100, len(kinetic_energies))

plt.plot(time_points, kinetic_energies, label=r'$Kinetic\ Energy$')
plt.plot(time_points, potential_energies, label=r'$Potential\ Energy$')
plt.plot(time_points, total_energies, label=r'$Total\ Energy$')
plt.xlabel(r'$Time$')
plt.ylabel(r'$Energy$')
plt.title(r'$Залежність\ середніх\ значень\ енергій\ від\ часу$')
plt.legend()
plt.grid(True)
plt.show()

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

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






start_time = time.time()


initial_state = init_state()
current_state = initial_state
displacements = []
time_points = []

current_time = 0
while current_time < total_time:
    print(f'current_time={current_time}\n')
    rt2_value = Rt2(current_state, initial_state)
    displacements.append(rt2_value)
    time_points.append(current_time)

    current_state = calc_next_state(current_state)
    current_time += DELTA_T


plt.figure(figsize=(10, 6))
plt.loglog(time_points, displacements, c='red', marker='o', linestyle='', markersize=5)
plt.xlabel(r'$Time$')
plt.ylabel(r'$Mean\ Square\ Displacement$')
plt.title(r'$Середньоквадратичне\ зміщення\ для\ всіх\ частинок$')
plt.grid(True)
plt.show()

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

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