Математика / Интерполяция

Интерполяция Акимы: расчёт наклона по четырём секущим

Метод Акимы подбирает наклон кубика в узле по четырём соседним секущим. Вес каждой из двух центральных секущих зависит от изменения наклона с противоположной стороны.

Опубликовано: Обновлено:

Формула

$$m_i=\frac{|\delta_{i+1}-\delta_i|\delta_{i-1}+|\delta_{i-1}-\delta_{i-2}|\delta_i}{|\delta_{i+1}-\delta_i|+|\delta_{i-1}-\delta_{i-2}|}$$
Схема Секущие вокруг двух узлов Акимы
01234505101521yxСекущие

Через шесть исходных узлов проведена ломаная. Два рассчитанных наклона m₂ и m₃ равны 4, поэтому интерполянт на отрезке от 2 до 3 совпадает с его центральным звеном.

Обозначения

$mᵢ$
производная в выбранном узле, ед. y / ед. x
$δⱼ$
секущий наклон между xⱼ и xⱼ₊₁, ед. y / ед. x
$i$
номер узла, в котором находят наклон, —
$j$
номер одного из четырёх интервалов, —

Условия применения

Формула записана для исходных весов Акимы и внутреннего узла, имеющего два промежутка по обе стороны. Абсциссы строго возрастают. Дробь применяется при ненулевой сумме весов; при двух нулевых весах здесь используют полусумму центральных секущих.

Ограничения

На краях нужны дополнительные правила продолжения последовательности секущих. У модифицированной схемы makima другие веса, и приведённое выражение к ней не относится. Метод может создать немонотонный участок даже при монотонных данных.

Подробное объяснение

В узле известны два ближайших секущих наклона, но для выбора касательной метод смотрит ещё на два соседних отрезка. Если справа две секущие совпадают, а слева меняются, правый наклон получает весь вес. Так прямой участок может сохранить свой наклон у стыка с более изменчивым окружением.

Левая центральная секущая умножается на абсолютную разность правых наклонов. Правая умножается на абсолютную разность левых. Перекрёстное расположение весов отражает идею локальной устойчивости: чем резче меняется наклон на одной стороне, тем меньше доверия к её ближайшей секущей.

При ненулевой сумме весов наклон лежит между центральными секущими. Это свойство относится только к производной в узле. Для монотонности кубика приходится учитывать ещё наклон на другом его конце и величину приращения на всём интервале.

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

Как пользоваться формулой

  1. Выберите узел с двумя интервалами по обе стороны и посчитайте четыре секущие с учётом длин шагов.
  2. Вычтите два правых наклона и возьмите модуль, затем проделайте то же для двух левых наклонов.
  3. Умножьте левую центральную секущую на правый вес, а правую секущую умножьте на левый вес.
  4. Разделите сумму произведений на сумму весов или обработайте вырождение по выбранному правилу.

Историческая справка

В 1970 году Хироси Акима опубликовал в Journal of the ACM работу о построении гладкой интерполяционной кривой локальными процедурами. В постановке с табличными данными локальный расчёт наклонов позволял строить кубические куски, не связывая вычисление каждой производной с системой по всей таблице.

Название Akima сохранилось в современных библиотеках. Документация SciPy для одного класса отдельно показывает исходную формулу и модификацию makima с дополнительными слагаемыми в весах. Для чтения старых алгоритмов и воспроизведения современных расчётов это различие существенно: при тех же четырёх секущих варианты могут вернуть разные производные.

Историческая линия формулы

Локальный метод кубической интерполяции связан с работой Хироси Акимы 1970 года. В формуле здесь используются исходные веса по абсолютным разностям соседних секущих, а современный вариант makima указан отдельно, поскольку изменяет эти веса.

Пример

На равномерной сетке x = 0, 1, 2, 3, 4, 5 значения y равны 0, 1, 3, 7, 11, 21. Секущие последовательно равны 1, 2, 4, 4 и 10. При x=2 правый вес |4−4|=0, левый |2−1|=1, поэтому m₂=(0·2+1·4)/1=4. При x=3 веса |10−4|=6 и |4−2|=2, поэтому m₃=(6·4+2·4)/8=4. На отрезке [2; 3] обе касательные совпадают с секущей. Кубик Эрмита здесь сводится к прямой y=3+4(x−2), поэтому при x=2,5 выходит y=5. Контроль: прямая возвращает 3 и 7 на концах, а её наклон равен четырём единицам y на единицу x. Отсчёт 2,5 находится посередине между двумя выбранными узлами.

Частая ошибка

Подстановка приращений y вместо наклонов даёт ошибку на сетке с разными шагами. Прямое умножение левой секущей на левый вес разрушает выбор в пользу более стабильной стороны. Если все четыре секущие равны, результат вычисляют без деления. Свойства makima нельзя переносить на исходные веса Акимы.

Практика

Задачи с решением

Ненулевые веса с обеих сторон

Условие. Четыре секущие вокруг узла равны 0, 2, 5 и 6. Найдите производную.

Решение. Правый вес |6−5|=1, левый |2−0|=2. Производная (1·2+2·5)/(1+2)=12/3.

Ответ. m = 4.

Два нулевых веса

Условие. Четыре секущие равны 1, 1, 3, 3. Примените правило полусуммы.

Решение. Разности справа и слева нулевые. Центральные значения равны 1 и 3, их полусумма (1+3)/2 равна 2.

Ответ. m = 2.

Дополнительные источники

  • https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.Akima1DInterpolator.html
  • https://dl.acm.org/doi/10.1145/321607.321609

Связанные формулы

Математика

Кубическая интерполяция Эрмита по значениям и наклонам

$H(x)=h_{00}(s)y_0+h_{10}(s)\ell m_0+h_{01}(s)y_1+h_{11}(s)\ell m_1,\quad s=\frac{x-x_0}{\ell},\quad \ell=x_1-x_0$

Кубический полином Эрмита согласует значения и первые производные на двух концах интервала. Четыре заданных условия определяют единственную кубическую или более простую кривую между узлами.

Математика

Линейная интерполяция: значение между двумя узлами

$P(x)=y_0+\frac{y_1-y_0}{x_1-x_0}(x-x_0)=(1-t)y_0+ty_1,\quad t=\frac{x-x_0}{x_1-x_0}$

Линейная интерполяция оценивает значение между двумя заданными точками по соединяющей их прямой. Доля положения внутри интервала задаёт веса концов; оценка учитывает наклон, но не кривизну исходного графика.

Математика

Внутренний наклон PCHIP через гармоническое среднее

$m_i=\begin{cases}0,&\delta_{i-1}\delta_i\le0,\\\dfrac{w_1+w_2}{w_1/\delta_{i-1}+w_2/\delta_i},&\delta_{i-1}\delta_i>0,\end{cases}\quad w_1=2h_i+h_{i-1},\quad w_2=h_i+2h_{i-1}$

PCHIP вычисляет производную во внутреннем узле по двум секущим. При общем знаке он использует взвешенное гармоническое среднее, а у плато или разворота задаёт нуль.