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

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

def euler_method_with_delay(T_initial, T_ambient, k, delta_t, num_iterations, cream_effect=0, delay_time=0):

    time = [0]

    if cream_effect != 0 and delay_time == 0:
        temperature = [T_initial - cream_effect]
    else:
        temperature = [T_initial]

    relaxation_time = None

    for i in range(1, num_iterations + 1):
        if (i * delta_t == delay_time) and cream_effect != 0:
            temperature[-1] -= cream_effect

        if i * delta_t >= delay_time and delay_time != 0:

            dT_dt = -k * (temperature[-1] - T_ambient)
            temp = (temperature[-1] + dT_dt * delta_t)

        else:
            dT_dt = -k * (temperature[-1] - T_ambient)
            temp = (temperature[-1] + dT_dt * delta_t)

        temperature.append(temp)
        time.append(i * delta_t)

        if relaxation_time is None and abs(temp - T_ambient) <= abs(T_initial - T_ambient) * 0.37:

            relaxation_time = i * delta_t


    return time, temperature, relaxation_time

def exact_solution(T_initial, T_ambient, k, time, cream_effect=0):
    C = T_ambient + (T_initial - T_ambient - cream_effect) * np.exp(-k * time)
    return C

# Початкові умови
T_initial = 83
T_ambient = 18
k = 0.039

dt = 0.01
dt_values = [1, 0.1, 0.01, 0.001, 0.0001]
num_iterations = 30

cream_effect = 5

plt.figure(figsize=(10, 6))
plt.xlabel('Час, хв')
plt.ylabel('Температура, градуси')
plt.title('Охолодження кави з різними значеннями k та ефектом вершків')
plt.grid()

num_iterations_dt = int(30 / dt)
time_euler, temperature_euler, _ = euler_method_with_delay(T_initial, T_ambient, k, dt, num_iterations_dt)
time_exact = np.arange(0, (num_iterations_dt + 1) * dt, dt)
temperature_exact = exact_solution(T_initial, T_ambient, k, time_exact)

plt.plot(time_euler, np.array(temperature_euler), label='Без вершків')
plt.plot(time_euler, np.array(temperature_exact), linestyle='--', label='Точне рішення без вершків')

time_cream, temperature_cream, _ = euler_method_with_delay(T_initial, T_ambient, k, dt, num_iterations_dt, cream_effect)
temperature_exact_cream = exact_solution(T_initial, T_ambient, k, time_exact, cream_effect)
plt.plot(time_exact, np.array(temperature_exact_cream), linestyle='-', label='З вершками')
plt.legend()
plt.show()


plt.figure(figsize=(10, 6))
plt.xlabel('Час, хв')
plt.ylabel('Абсолютна похибка')
for dt in dt_values:
    num_iterations_dt = int(30 / dt)
    time_euler, temperature_euler, _ = euler_method_with_delay(T_initial, T_ambient, k, dt, num_iterations_dt)
    time_exact = np.arange(0, (num_iterations_dt + 1) * dt, dt)
    temperature_exact = exact_solution(T_initial, T_ambient, k, time_exact)
    absolute_error = np.abs(np.array(temperature_exact) - np.array(temperature_euler))
    cutoff_index = find_cutoff_index(time_euler, absolute_error)

    start_index = 0
    for i in range(len(time_euler)):
        if time_euler[i] >= 1:
            start_index = i
            break

    plt.plot(time_euler[start_index:cutoff_index], absolute_error[start_index:cutoff_index], label=f'dt = {dt}')

plt.xlabel('Час, хв')
plt.ylabel('Абсолютна похибка')
plt.title('Графіки абсолютної похибки від часу для різних значень dt')
plt.yscale('log')
plt.legend()
plt.grid()
plt.show()




real_time_values = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15])
real_temperature_values = np.array([83, 77.7, 75.1, 73, 71.1, 69.4, 67.8, 66.4, 64.7, 63.4, 62.1, 61, 59.9, 58.7, 57.8, 56.6])
temp_df = pd.DataFrame({'Time': real_temperature_values, 'T': real_time_values})

T_desired = 50
k_values = [0.03, 0.038, 0.05]


time_values = np.arange(0, 30, 1)
temperature_values = []


for k in k_values:
    temperature = [T_initial]
    for t in range(1, 30):
        dT_dt = -k * (temperature[-1] - T_ambient)
        temperature.append(temperature[-1] + dT_dt)

    temperature_values.append(temperature)


plt.figure(figsize=(10, 6))
for i, k in enumerate(k_values):
    plt.plot(time_values, temperature_values[i], label=f'r = {k}')

plt.scatter(real_time_values, real_temperature_values, color='black', marker='s', label='Експеремент')
plt.axhline(y=T_desired, color='r', linestyle='--', label='Комфортна температура')
plt.xlabel('Час, хв')
plt.ylabel('Температура, градуси')
plt.title('Охолодження кави з різними значеннями r')
plt.legend()
plt.grid()
plt.show()

delay_time = 3

plt.figure(figsize=(10, 9))
plt.xlabel('Час, хв')
plt.ylabel('Температура, градуси')
plt.title('Охолодження кави', size=24)
plt.grid()

num_iterations_dt = int(30 / dt)
time_euler, temperature_euler, _ = euler_method_with_delay(T_initial, T_ambient, k, dt, num_iterations_dt)
plt.plot(time_euler, np.array(temperature_euler), label='Без вершків')

temperature_exact_cream = exact_solution(T_initial, T_ambient, k, time_exact, cream_effect)
plt.plot(time_exact, np.array(temperature_exact_cream), linestyle='-', label='З вершками')

time_cream, temperature_cream, _ = euler_method_with_delay(T_initial, T_ambient, k, dt, num_iterations_dt, cream_effect=cream_effect, delay_time=delay_time)
plt.plot(time_cream, np.array(temperature_cream), linestyle='--', label=f'З вершками (затримка {delay_time} хв)')

time_cup, temperature_cup, _ = euler_method_with_delay(T_initial-10, T_ambient, k, dt, num_iterations_dt)
plt.plot(time_cup, np.array(temperature_cup), linestyle='-', label=f'Товща чашка(кава без вершків)')

time_s, temperature_s, _ = euler_method_with_delay(T_initial, 30, k, dt, num_iterations_dt)
plt.plot(time_s, np.array(temperature_s), linestyle='-', label=f'Кава без вершків в літній час')



plt.legend()
plt.show()



T_initial_values = np.linspace(80, 100, 11)
plt.figure(figsize=(10, 5))
for r in [0.038, 0.05, 0.03]:
    relaxation_times = []
    for T_initial in T_initial_values:
        _, _, relaxation_time = euler_method_with_delay(T_initial, T_ambient, r, dt, num_iterations_dt)
        relaxation_times.append(relaxation_time)
    plt.plot(T_initial_values, relaxation_times, linestyle='-', marker='o', label=f'Час релаксації(r={r})')
plt.xlabel('Початкова температура, градуси')
plt.ylabel('Час релаксації, хв')
plt.title('Залежність часу релаксації від початкової температури')
plt.grid()
plt.legend()
plt.show()

T_ambient_values = np.linspace(0, 50, 11)
plt.figure(figsize=(10, 5))
for r in [0.038, 0.05, 0.03]:
    relaxation_times = []
    for T_ambient in T_ambient_values:
        _, _, relaxation_time = euler_method_with_delay(T_initial, T_ambient, r, dt, num_iterations_dt)
        relaxation_times.append(relaxation_time)
    plt.plot(T_ambient_values, relaxation_times, linestyle='-', marker='o', label=f'Час релаксації(r={r})')


plt.xlabel('Температура навколишнього середовища, градуси')
plt.ylabel('Час релаксації, хв')
plt.title('Залежність часу релаксації від температури навколишнього середовища')
plt.grid()
plt.legend()
plt.show()