Интегралы и дифференциальные уравнения

Модуль scipy.integrate решает две задачи: считает интегралы (площадь под графиком, накопленная величина) и решает дифференциальные уравнения (как величина меняется со временем).

Интеграл функции: quad

import math
from scipy.integrate import quad

value, err = quad(lambda x: math.sin(x), 0, math.pi)
print(value)   # 2.0
print(err)     # оценка погрешности, около 2e-14

quad(f, a, b) возвращает пару: значение интеграла и оценку его погрешности. Обычно нужно только значение — quad(...)[0] или распаковка, как выше.

Границы могут быть бесконечными:

print(quad(lambda x: math.exp(-x), 0, math.inf)[0])   # 1.0

math.inf, float("inf") и numpy.inf — одна и та же бесконечность.

Функция может зависеть от параметров — их берут из окружающих переменных или передают через args:

f = lambda x, p: x ** p
print(quad(f, 0, 1, args=(3,))[0])   # 0.25

Интеграл по таблице значений: trapezoid

Если функции нет, а есть только измерения в точках, интеграл приближают формулой трапеций: соседние точки соединяются отрезками, и считается площадь под ломаной.

from scipy.integrate import trapezoid

x = [0, 1, 2, 4]
y = [0, 1, 4, 16]
print(trapezoid(y, x))   # 0.5 + 2.5 + 20 = 23.0

Внимание на порядок: сначала значения y, потом точки x. Точки не обязаны идти с равным шагом. Есть и более точная для гладких функций формула Симпсона — simpson(y, x=x).

Дифференциальные уравнения: solve_ivp

Дифференциальное уравнение описывает, как быстро меняется величина. Например, остывание чая: скорость падения температуры пропорциональна разнице с комнатой, T' = -0.1 (T - 20).

solve_ivp нужны:

1. функция f(t, y), которая по времени и текущему состоянию возвращает производную;
2. отрезок времени (t0, t1);
3. начальное состояние — список, даже если величина одна.

from scipy.integrate import solve_ivp

sol = solve_ivp(lambda t, y: [-0.1 * (y[0] - 20)], (0, 30), [90], t_eval=[0, 10, 20, 30])
print(sol.t)      # [ 0. 10. 20. 30.]
print(sol.y[0])   # температура в эти моменты

• t_eval — моменты, в которые нужен ответ (внутри отрезка времени);
• sol.y — таблица: строка — компонента состояния, столбец — момент времени;
• rtol, atol — допустимая относительная и абсолютная погрешность. По умолчанию они довольно грубые (rtol=1e-3); для точных ответов ставят, например, rtol=1e-9, atol=1e-9.

Уравнения высокого порядка

Если в уравнении есть вторая производная, в состояние добавляют первую. Например, движение тела на пружине x'' = -x: состояние [x, v], где v = x', а производная состояния — [v, -x]:

def rhs(t, s):
    x, v = s
    return [v, -x]

sol = solve_ivp(rhs, (0, 10), [1, 0], t_eval=[0, 1, 2], rtol=1e-9, atol=1e-9)
print(sol.y[0])   # ≈ cos(t): [1.  0.5403  -0.4161]

Итог

• quad(f, a, b) → (значение, погрешность); границы могут быть math.inf;
• trapezoid(y, x) — интеграл по точкам, сначала y, потом x;
• solve_ivp(f, (t0, t1), y0, t_eval=..., rtol=..., atol=...) — решение уравнения y' = f(t, y); ответ в sol.y;
• уравнение второго порядка → система первого порядка: [x, x'].