Интерполяция и линейные системы

Эта тема про две частые инженерные задачи: «достроить» функцию между известными точками и решить систему линейных уравнений.

Интерполяция

Есть значения функции в нескольких точках — замеры, табличные данные. Нужно оценить её между ними. Это и есть интерполяция.

Линейная: `interp1d`

Соседние точки соединяются отрезками:

from scipy.interpolate import interp1d

hours = [0, 6, 12, 18]
temp = [10, 8, 20, 14]
f = interp1d(hours, temp)       # f — новая функция
print(f(3))                     # 9.0
print(f([9, 15]))               # [14. 17.]

interp1d возвращает функцию, которую можно вызывать от числа или от списка. Точки x должны быть разными; если они не отсортированы, функция отсортирует их сама.

Параметр kind выбирает способ: "linear" (по умолчанию), "nearest" (значение ближайшей точки), "cubic" и другие. За пределами данных по умолчанию будет ошибка; продолжить функцию можно параметром fill_value="extrapolate".

Для линейного случая есть и короткая функция numpy: np.interp(t, x, y).

Сплайны: `CubicSpline`

Ломаная имеет «изломы» в точках. Если нужна гладкая кривая (траектория, плавная анимация, профиль дороги), используют кубический сплайн: между каждыми соседними точками — свой многочлен третьей степени, и куски состыкованы так, что совпадают значения, первые и вторые производные.

from scipy.interpolate import CubicSpline

cs = CubicSpline([0, 1, 2, 3], [0, 1, 0, 1], bc_type="natural")
print(cs(1.5))       # значение
print(cs(1.5, 1))    # первая производная (скорость изменения)
print(cs(1.5, 2))    # вторая производная

У сплайна остаётся свобода на краях, и её фиксируют граничными условиями bc_type:

"not-a-knot" (по умолчанию) — первые два и последние два куска — один и тот же многочлен;
"natural" — вторая производная на концах равна нулю: кривая «выпрямляется» на краях.

Линейные системы: `scipy.linalg`

Система

2x + y = 3 x + 3y = 5

в матричном виде: A x = b, где A = [[2, 1], [1, 3]], b = [3, 5].

from scipy import linalg

A = [[2, 1], [1, 3]]
b = [3, 5]
print(linalg.solve(A, b))   # [0.8 1.4]
print(linalg.det(A))        # 5.0 — определитель
print(linalg.inv(A))        # обратная матрица

строки A — уравнения, столбцы — неизвестные;
solve надёжнее и быстрее, чем inv(A) @ b;
если определитель равен нулю, решения нет или их бесконечно много — solve выдаст ошибку LinAlgError.

scipy.linalg похож на numpy.linalg, но умеет больше: разложения матриц (lu, qr, cholesky), матричную экспоненту expm и т. д.

Итог

interp1d(x, y) — функция-ломаная; kind="nearest", "cubic";
CubicSpline(x, y, bc_type="natural") — гладкая кривая; cs(t, 1) — производная;
linalg.solve(A, b) — решение A x = b; det, inv — определитель и обратная матрица.

Интерполяция и линейные системы