Недостатки классической интерполяции Лагранжа

Как мы уже изучали в разделе об аппроксимации функций, классический интерполяционный многочлен Лагранжа строится так, чтобы его график строго проходил через заданный набор дискретных точек (узлов). Однако в инженерном проектировании, компьютерной графике (CAD) и вычислительной гидродинамике простого прохождения кривой через точки часто бывает недостаточно. Представьте, что вы проектируете профиль автомобильной дороги или траекторию движения резца станка с ЧПУ. Если вы используете полином Лагранжа, кривая пройдет через все контрольные точки, но в самих этих точках ее наклон (первая производная) или кривизна (вторая производная) могут изменяться непредсказуемым и диким образом, порождая рывки и изломы скорости.

Инженеру необходимо жестко контролировать не только то, ГДЕ должна оказаться кривая в данный момент времени, но и КУДА она должна быть направлена (какова ее скорость и ускорение). Эту фундаментальную задачу решает метод интерполяции Шарля Эрмита. Алгоритм Эрмита обобщает интерполяцию Лагранжа на случай, когда в узлах сетки известны (и должны строго соблюдаться) не только значения самой функции, но и значения ее производных до заданного порядка.

Разделенные разности с кратными узлами

Построение многочлена Эрмита математически наиболее элегантно реализуется с использованием аппарата разделенных разностей (формула Ньютона для интерполяции). Гениальность подхода заключается в фиктивном (виртуальном) дублировании узлов сетки. Если в каком-то узле x_i нам задано не только значение функции, но и значение ее первой производной, мы записываем этот узел в таблицу разделенных разностей дважды (как кратный узел).

Когда алгоритм пытается вычислить разделенную разность первого порядка (скорость изменения функции) между двумя одинаковыми продублированными узлами, в знаменателе дроби возникает ноль: (x_i - x_i). Но мы знаем из математического анализа, что предел отношения приращения функции к приращению аргумента — это и есть производная! Поэтому алгоритм просто заменяет эту неопределенность 0/0 на заданное нам известное значение производной в этой точке. Аналогично, если задана вторая производная (ускорение), узел дублируется трижды. В результате мы получаем единый полином высокой степени, который идеально ложится на все точки и векторы касательных.

Кубические сплайны Эрмита и сплайны Акимы

Глобальный многочлен Эрмита высокой степени (как и многочлен Лагранжа) страдает от феномена Рунге — катастрофических осцилляций на краях отрезка. Поэтому в индустрии повсеместно применяются кусочно-полиномиальные кривые — кубические сплайны Эрмита.

На каждом интервале между двумя точками строится отдельный полином 3-й степени, который жестко фиксируется 4 условиями: двумя значениями функции по краям и двумя значениями производных по краям. В отличие от обычного кубического сплайна, где производные вычисляются алгоритмом глобально (через решение СЛАУ для обеспечения гладкости второй производной во всех точках), в сплайнах Эрмита производные задаются локально и вручную инженером (или локальными эвристиками). Одной из самых знаменитых реализаций является Сплайн Акимы (Hiroshi Akima, 1970). Алгоритм Акимы вычисляет производные в узлах локально, опираясь на геометрию нескольких соседних точек. В результате сплайн Акимы практически не подвержен паразитным выбросам (overshoots) и изломам в местах резкого изменения данных, являясь золотым стандартом для интерполяции физических табличных данных (например, термодинамических свойств газов) без искусственных осцилляций.

Фундамент итерационной математики

В самых первых статьях нашего цикла мы рассматривали метод дихотомии и метод Ньютона для поиска корней нелинейных уравнений вида f(x) = 0. Метод Ньютона быстр, но требует вычисления производных и часто расходится при плохом выборе начального приближения. Существует еще один, невероятно элегантный и концептуально прозрачный класс алгоритмов — метод простой итерации (или метод последовательных приближений). Его глубокая математическая суть базируется не на геометрии касательных, а на фундаментальных понятиях топологии и функционального анализа.

Чтобы применить этот метод, исходное уравнение f(x) = 0 необходимо алгебраически преобразовать к эквивалентному виду: x = phi(x). Это можно сделать бесконечным числом способов (например, просто прибавив x к обеим частям, или выразив x из сложного логарифма). Теперь наша задача поиска корня превратилась в задачу поиска неподвижной точки отображения phi (fixed point). Неподвижная точка — это такое значение x, которое при подстановке в функцию phi выдает на выходе в точности то же самое значение x. Вычислительный алгоритм становится тривиальным: мы берем начальное случайное приближение x0, вычисляем x1 = phi(x0), затем x2 = phi(x1) и так далее до бесконечности.

Сжимающие отображения и теорема Банаха

Главный вопрос вычислительной математики: при каких условиях эта бесконечная последовательность x_n действительно сойдется к истинному корню, а не улетит в бесконечность и не зациклится? Ответ дает знаменитая теорема Банаха о неподвижной точке (принцип сжимающих отображений), доказанная польским математиком Стефаном Банахом в 1922 году.

Теорема гласит: если функция phi(x) является сжимающим отображением на замкнутом интервале (то есть расстояние между любыми двумя точками после применения функции строго уменьшается с коэффициентом сжатия q < 1), то на этом интервале существует ровно одна неподвижная точка, и метод простой итерации гарантированно сойдется к ней из абсолютно любого начального приближения! С точки зрения дифференциального исчисления, условие сжатия означает, что модуль производной функции phi(x) в окрестности корня должен быть строго меньше единицы. Чем ближе производная к нулю, тем меньше коэффициент сжатия q, и тем феноменально быстрее алгоритм будет сходиться. Если же производная по модулю больше единицы (отображение растягивающее), алгоритм неминуемо разойдется, как бы близко к корню мы ни стартовали. Искусство применения метода простой итерации заключается именно в таком хитром алгебраическом преобразовании f(x) = 0 в x = phi(x), чтобы искусственно сделать производную phi(x) максимально близкой к нулю.