Řešení rovnice ochlazování

Naivní metoda, využívající Eulerovu metodu, spočívá v tom, že čas rozdělíme na malé intervaly délky \(h\) a v každém kroku aktualizujeme teplotu kávy podle vzorce

\[T(t+h) = T(t) - k \cdot h \cdot (T(t) - T_{\infty})\]

kde \(T(t)\) je teplota kávy v aktuálním čase, \(T(t+h)\) je teplota kávy v následujícím časovém okamžiku, \(k\) je konstanta ochlazování, \(h\) je délka časového kroku a \(T_{\infty}\) je teplota okolí. Tento vzorec je odvozen z původní diferenciální rovnice

\[\frac{\mathrm dT}{\mathrm dt} = -k (T - T_{\infty})\]

pomocí aproximaci derivace dopřednou diferencí, tj. z rovnice

\[\frac{T(t+h) - T(t)}{h} \approx -k (T(t) - T_{\infty}).\]

import numpy as np
import matplotlib.pyplot as plt
# Parametry
T0 = 100  # Počáteční teplota kávy
T_inf = 20  # Teplota okolí
k = 0.1  # Koeficient ochlazování
h = 0.1  # Krok časového intervalu
max_time = 60  # Maximální čas pro simulaci
# Inicializace
times = np.arange(0, max_time, h)
T = np.zeros_like(times)
T[0] = T0
# Numerické řešení
for i in range(1, len(times)):
    T[i] = T[i-1] - k * h * (T[i-1] - T_inf)
# Graf
plt.plot(times, T)
plt.xlabel('Čas (min)')
plt.ylabel('Teplota kávy (°C)')
plt.title('Ochlazování kávy podle Newtonova zákona')
plt.grid()
plt.show()