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

import numpy as np
import matplotlib.pyplot as plt

# Параметри пластини
Lx = 100
Ly = 100
time_max = 10
T1 = 0
T2 = 1
D = 10
dt = 0.001

# Параметри джерела
rect_x = Lx // 4
rect_y = Ly // 4
rect_width = Lx // 2
rect_height = Ly // 2

T = np.ones((Ly + 2, Lx + 2)) * T1
T[rect_y:rect_y + rect_height, rect_x:rect_x + rect_width] = T2

def calc_next(T, D, dt):
    TempMas = np.copy(T)
    for j in range(1, Ly + 1):
        for i in range(1, Lx + 1):
            TempMas[j, i] = T[j, i] + D * dt * (
                T[j, i + 1] + T[j, i - 1] + T[j + 1, i] + T[j - 1, i] - 4 * T[j, i]
            )
    return TempMas


T_list = []
time_points = []
diff_history = []
temperature_history = []
t = dt
while t <= time_max:
    T = calc_next(T, D, dt)
    average_temp = np.mean(T)
    temperature_history.append(average_temp)
    time_points.append(t)
    diff = np.mean(T**2) - average_temp**2
    diff_history.append(diff)
    if round(t * 1000) % 1000 == 0 or t == dt:
        T_list.append(np.copy(T))
    t += dt


fig, axs = plt.subplots(nrows=len(T_list)//3+1, ncols=3, figsize=(15, 5*(len(T_list)//3)+1))
for idx, t_value in enumerate(np.linspace(0, time_max, len(T_list))):
    ax = axs[idx//3, idx%3]
    im = ax.imshow(T_list[idx], cmap='inferno', interpolation='nearest', vmin=T1, vmax=T2)
    ax.set_title(fr"$Time: {round(t_value, 6)}$")
    fig.colorbar(im, ax=ax)
    ax.set_xticks([])
    ax.set_yticks([])
    ax.grid(False)


plt.figure(figsize=(10, 5))
plt.plot(time_points, temperature_history, color='red')
plt.xlabel(fr'$Час$')
plt.ylabel(fr'$Середня\ температура$')
plt.title(fr'$Залежність\ середньої\ температури\ від\ часу$')
plt.grid(True)
plt.show()


plt.figure(figsize=(10, 5))
plt.plot(time_points, diff_history)
plt.xlabel(fr'$Час$')
plt.ylabel(fr'$<T^2>-<T>^2$')
plt.title(fr"$Залежність\ <T^2>-<T>^2\ від\ часу$")
plt.grid(True)
plt.tight_layout()
plt.show()









T = np.ones((Ly + 2, Lx + 2)) * T1
T[rect_y:rect_y + rect_height, rect_x:rect_x + rect_width] = T2

def calc_next2(T, D, dt):
    TempMas = np.copy(T)
    for j in range(1, Ly + 1):
        for i in range(1, Lx + 1):
            if not (rect_y <= j < rect_y + rect_height and rect_x <= i < rect_x + rect_width):
                TempMas[j, i] = T[j, i] + D * dt * (T[j, i + 1] + T[j, i - 1] + T[j + 1, i] + T[j - 1, i] - 4 * T[j, i])
    return TempMas

T_list = []
time_points = []
diff_history = []
temperature_history = []
t = dt
while t <= time_max:
    T = calc_next2(T, D, dt)
    average_temp = np.mean(T)
    temperature_history.append(average_temp)
    time_points.append(t)
    diff = np.mean(T**2) - np.mean(T)**2
    diff_history.append(diff)
    if round(t * 1000) % 1000 == 0 or t == dt:
        T_list.append(T)
    t += dt



fig, axs = plt.subplots(nrows=len(T_list)//3+1, ncols=3, figsize=(15, 5*(len(T_list)//3+1)))
for idx, t_value in enumerate(np.linspace(0, time_max, len(T_list))):
    ax = axs[idx//3, idx%3]
    im = ax.imshow(T_list[idx], cmap='inferno', interpolation='nearest', vmin=T1, vmax=T2)
    ax.set_title(fr"$Time: {round(t_value, 6)}$")
    fig.colorbar(im, ax=ax)
    ax.set_xticks([])
    ax.set_yticks([])
    ax.grid(False)



plt.figure(figsize=(10, 5))
plt.plot(time_points, temperature_history)
plt.xlabel(fr'$Час$')
plt.ylabel(fr'$Середня\ температура$')
plt.title(fr"$Залежність\ середньої\ температури\ від\ часу$")
plt.grid(True)

plt.figure(figsize=(10, 5))
plt.plot(time_points, diff_history)
plt.xlabel(fr'$Час$')
plt.ylabel(fr'$<T^2>-<T>^2$')
plt.title(fr"$Залежність\ <T^2>-<T>^2 \ від\ часу$")
plt.grid(True)

plt.tight_layout()
plt.show()