Практические кейсы: остывание, RC-цепь, химия

Четыре маленькие задачи из физики и химии, в каждой из которых дифференциальное уравнение описывает что-то из реальной жизни — и решается в несколько строк кода.

Уравнение первого порядка с релаксацией — модель вида y' = -k(y - y_цель), описывающая экспоненциальное приближение величины к равновесному значению; встречается в остывании, заряде конденсатора и распаде.

Многие явления вокруг нас подчиняются одному и тому же математическому шаблону: величина меняется тем быстрее, чем дальше она от равновесия. Остывающий кофе, заряжающийся конденсатор, исчезающее в реакции вещество — всё это разные лица одного уравнения. Разберём четыре кейса, в каждом запустим численное решение и сверим его с аналитикой.

Кейс А: остывание кофе (закон Ньютона)

Закон охлаждения Ньютона: скорость остывания пропорциональна разности температур тела и окружающей среды. T' = -k(T - T_env). Здесь T — температура кофе, T_env — комнатная температура, k — коэффициент теплоотдачи. Аналитическое решение известно точно: T(t) = T_env + (T0 - T_env) * exp(-k*t). Это даёт нам эталон для проверки численного метода.

import math

T_env = 20.0
k = 0.1
T0 = 90.0

def dT(T):
    return -k * (T - T_env)

def rk4_scalar(f, y, h):
    k1 = f(y)
    k2 = f(y + 0.5*h*k1)
    k3 = f(y + 0.5*h*k2)
    k4 = f(y + h*k3)
    return y + (h/6.0)*(k1 + 2*k2 + 2*k3 + k4)

T = T0
h = 1.0
print("мин   RK4      точно")
for minute in range(0, 31):
    if minute % 5 == 0:
        exact = T_env + (T0 - T_env)*math.exp(-k*minute)
        print(str(minute).rjust(3), "  ", round(T,4), "  ", round(exact,4))
    T = rk4_scalar(dT, T, h)

Вывод:

мин   RK4      точно
  0    90.0    90.0
  5    62.4636    62.4636
 10    45.7392    45.7392
 15    35.5837    35.5837
 20    29.4153    29.4153
 25    25.6688    25.6688
 30    23.394    23.394

RK4 совпадает с аналитикой до четвёртого знака — для гладкого экспоненциального решения метод почти идеален. Кофе с 90°C за полчаса остывает примерно до 23°C, асимптотически приближаясь к комнатным 20°C, но никогда их точно не достигая.

Кейс Б: зарядка конденсатора в RC-цепи

Если подать напряжение Vin на цепь из резистора R и конденсатора C, конденсатор заряжается по уравнению V' = (Vin - V) / (R*C). Произведение R*C называют постоянной времени tau — за время tau напряжение достигает約 63% от Vin. Аналитика: V(t) = Vin * (1 - exp(-t / (R*C))).

import math

Vin = 5.0
R = 1000.0
C = 0.001
tau = R * C

def dV(V):
    return (Vin - V) / tau

def rk4_scalar(f, y, h):
    k1 = f(y)
    k2 = f(y + 0.5*h*k1)
    k3 = f(y + 0.5*h*k2)
    k4 = f(y + h*k3)
    return y + (h/6.0)*(k1 + 2*k2 + 2*k3 + k4)

V = 0.0
h = 0.1
print("tau =", tau, "сек")
print("  t      V(RK4)   V(точно)")
t = 0.0
for step in range(0, 51):
    if step % 10 == 0:
        exact = Vin * (1 - math.exp(-t / tau))
        print(str(round(t,1)).rjust(4), "  ", round(V,4), "  ", round(exact,4))
    V = rk4_scalar(dV, V, h)
    t += h

Вывод:

tau = 1.0 сек
  t      V(RK4)   V(точно)
 0.0    0.0    0.0
 1.0    3.1606    3.1606
 2.0    4.3233    4.3233
 3.0    4.7511    4.7511
 4.0    4.9084    4.9084
 5.0    4.9663    4.9663

При tau=1 секунда напряжение к t=1 действительно достигает около 3.16 В — это и есть те самые 63% от 5 В. К t=5 (пять постоянных времени) конденсатор заряжен практически полностью. Та же экспоненциальная релаксация, что и у кофе, только теперь величина растёт к Vin, а не падает.

Кейс В: химическая кинетика A в B

Реакция первого порядка: вещество A превращается в B со скоростью, пропорциональной концентрации A. Уравнения для концентраций: [A]' = -k[A], [B]' = k[A]. Это система двух уравнений, но они связаны: сколько A исчезло, столько B появилось. Сумма [A] + [B] сохраняется. Аналитика: [A](t) = A0 * exp(-k*t), [B](t) = A0 * (1 - exp(-k*t)).

import math

k = 0.3
A0 = 1.0

def rates(state):
    A, B = state
    return [-k * A, k * A]

def rk4_step(f, state, h):
    k1 = f(state)
    s2 = [state[i] + 0.5*h*k1[i] for i in range(2)]
    k2 = f(s2)
    s3 = [state[i] + 0.5*h*k2[i] for i in range(2)]
    k3 = f(s3)
    s4 = [state[i] + h*k3[i] for i in range(2)]
    k4 = f(s4)
    return [state[i] + (h/6.0)*(k1[i] + 2*k2[i] + 2*k3[i] + k4[i]) for i in range(2)]

state = [A0, 0.0]
h = 0.5
print("  t      [A]      [B]    сумма")
t = 0.0
for step in range(0, 21):
    if step % 4 == 0:
        print(str(round(t,1)).rjust(4), "  ", round(state[0],4), "  ", round(state[1],4), "  ", round(state[0]+state[1],4))
    state = rk4_step(rates, state, h)
    t += h

Вывод:

  t      [A]      [B]    сумма
 0.0    1.0    0.0    1.0
 2.0    0.5488    0.4512    1.0
 4.0    0.3012    0.6988    1.0
 6.0    0.1653    0.8347    1.0
 8.0    0.0907    0.9093    1.0

Концентрация A падает экспоненциально, B растёт зеркально, а сумма остаётся равной 1.0 на каждом шаге — закон сохранения вещества соблюдён численно. Это хороший тест корректности: если сумма «поплыла», значит, в коде ошибка.

Как работает под капотом

Все четыре кейса — это релаксация к равновесию. Математически они идентичны уравнению y' = -k(y - y_цель), решение которого — экспонента. Для кофе цель T_env, для конденсатора Vin, для вещества A — ноль. Разница только в знаке и в том, растёт величина или убывает. Поэтому один и тот же численный метод (RK4) с одинаковой лёгкостью решает их все, а для гладких экспонент он даёт практически точный результат даже при крупном шаге.

Частые ошибки

• Неверный знак k. В остывании и распаде производная отрицательна (величина убывает). Перепутанный знак даёт экспоненциальный рост вместо затухания.
• Забытая постоянная времени. В RC-цепи делить надо именно на R*C, а не на R или C по отдельности; единицы измерения должны сходиться.
• Нарушение закона сохранения. В химии [A]+[B] обязано сохраняться. Проверяйте это как контроль ошибок.
• Слишком большой шаг при малом tau. Если постоянная времени мала, а шаг велик, решение может осциллировать или взорваться — об этом подробнее в теме жёсткости.

Зачем это нужно

Эти простые кейсы показывают, что дифференциальные уравнения — рабочий инструмент инженера и учёного, а не абстракция. Один шаблон релаксации описывает тепло, электричество и химию. Научившись распознавать его, вы будете моделировать самые разные процессы одним и тем же кодом. А привычка сверять численное решение с аналитикой и проверять законы сохранения — основа доверия к любой симуляции.

• Закон Ньютона, RC-цепь и реакция первого порядка — все это релаксация y' = -k(y - цель).
• RK4 для гладких экспонент совпадает с аналитикой до многих знаков.
• Постоянная времени tau = R*C задаёт скорость зарядки конденсатора (63% за tau).
• Сохранение суммы [A]+[B] — встроенный тест корректности численной схемы.