Эта статья не про NURBS. И не про методы Лагранжа и Ньютона. Про это и так написано достаточно. Эта статья про одну интересную технику интерполяции c несколькими нетривиальными примерами ее применения при разработке игр.
Введение
Говоря «полином», обычно подразумевают алгебраический полином, то есть функцию вещественного числа вида:
Существуют другие функции называемые полиномами, например, полиномы Фурье, но в этой статье они рассматриваться не будут. Равно как и полиномы комплексной переменной и нескольких переменных. Хотя эти математические объекты тоже очень интересны и заслуживают отдельного внимания.
Полиномы в программировании удобны благодаря нескольким своим свойствам. Во-первых, они легко записываются. Полином можно представить как упорядоченный список его коэффициентов. Операторы дифференцирования и интегрирования на таком объекте делаются очень просто, ведь производная полинома – это тоже полином.
Во-вторых, они легко считаются. По схеме Горнера, например: . Ровно два сложения и два умножения. Понятно, что какая степень полинома, столько умножений со сложениями и потребуется.
В-третьих, полиномы – функции непрерывно определенные на всей числовой прямой. То есть полином обязательно ставит в соответствие любому числу какое-то другое число. Нам не придется следить, например, чтобы аргумент был определен только на каком-то определенном интервале, как в случае с логарифмами, или избегал каких-то конкретных точек, как в случае с тангенсом.
Что еще надо знать о полиномах. Коэффициенты полинома – это значения его соответствующих производных в нуле. Логично же: при дифференцировании свободный член исчезает, на его место приходит коеффициент при аргументе первой степени, а сам аргумент в нуле в любой степени конечно же тоже ноль, значит значение полинома в нуле и есть его свободный член.
Ну а еще полиномы просто складываются, вычитаються и умножаются на скаляр. И все это очень легко программировать.
Теперь о интерполяции. Интерполяция – это такой способ аппроксимации, то есть нахождения функции приближенной к некоторой искомой, при котором интерполянт в исходных точках принимает исходные значения. То есть: имеем множество пар <точка, значение> , таких что интерполянт для всех этих пар. Это множество пар называется каркасом интерполяции.
Интерполяция может быть и функции нескольких переменных. При этом точки пар каркаса определяются не на числовой прямой, а на числовом пространстве , где k – число переменных функции.
Полиномиальной интерполяцией часто называют интерполяцию функции одной переменной полиномом Лагранжа. Это такой полином, который являясь интерполяционным для n пар сам имеет n-1 степень. Например, уравнение прямой для двух точек со значениями на числовой прямой, парабола для трех. Но не для всех трех. Если, например, из трех точек две совпадают, а их значения нет – найти такой полином не получится. По определению функции. А если две пары совпадают полностью, то для них можно найти бесконечное множество парабол, а не какую-то одну.
Но в остальном да: n точек – n-1 степень. Такая интерполяция в свое время использовалась довольно широко, но сейчас на практике ее применять можно ну очень ограниченно. Самая главная причина этому – склонность полинома Лагранжа к осцилляциям при увеличении числа точек каркаса интерполяции.
Для того, чтобы найти полином Лагранжа, используют в частности одноименный Лагранжа метод. Суть его в том, что интерполяционный полином представляется в виде суммы из n полиномов Pi(x), которые имеют корни во всех точках каркаса интерполяции, кроме xi, а в xi равняються yi.
Довольно просто: корень полинома – это ноль значения функции. Складываем n-1 нулей с yi – получаем yi. Следовательно, полином интерполяционный.
На самом деле, эту схему можно модифицировать для получения не лагранжевых полиномов, но пока подробно это рассматривать не стоит. А вот другую схему, которую можно использовать для получения, как Лагранжевых, так и не лагранжевых полиномов, пожалуй, стоит и рассмотреть.
Допустим, каркас включает в себя три различные точки (x1, x2, x3) и три соответствующих значения функции (y1, y2, y3) в них. Сам факт того, что интерполяционный полином P(x) определен в этих точках, означает:
P(x1)=y1
P(x2)=y2
P(x3)=y3
Или же:
Тут намеренно коэффициенты записаны после иксов, чтобы показать, что получается система линейных алгебраических уравнений. Относительно не иксов, конечно же, а коэффициентов. При этом все то, что мы знаем про СЛАУ, объясняет существование полинома Лагранжа, и отлично при этом иллюстрируется геометрически.
Через n различных точек на плоскости можно провести один и только один график полинома c n коэффициентами. Обратите внимание, это не обязательно лагранжев полином. Для нашего примера мы можем взять, например такие:
И даже:
Кроме того, если две точки каркаса совпадают, а значения нет – это означает, что СЛАУ не имеет решения, а если точки совпадают вместе со значениями, то система будет содержать два линейно зависимых, а точнее тождественных уравнения и, соответственно иметь бесконечное множество решений. И действительно, имея на плоскости две точки мы можем провести через них бесконечное множество парабол.
Но самое замечательное, что мы можем обогащать нашу систему новыми уравнениями, повышая, естественно, степень полинома, при этом оставляя полином интерполяционным. Такими уравнениями может быть, например, условие заданной наперед производной. Мы можем в нашем примере затребовать, чтобы в средней точке график касательная графика функции была параллельной оси Ox, то есть производная функции равнялась бы нулю:
Еще мы можем потребовать, например, чтобы интеграл интерполянта на некотором промежутке был определен числом. Или чтобы вторая производная в нескольких точках была определена одинаково, что может быть полезным для построения графиков с заданной кривизной. Практически любое условие, касающееся аналитических свойств функции, мы можем записать уравнением и добавить в систему.
Есть еще совсем экзотическое приложение полиномиальной интерполяции: мы можем потребовать, чтобы в некотором числе точкек соблюдалось определенное дифференциальное уравнение. Конечно, мы не можем требовать такого на всей области определения, соответственно, точно решать диффуры полиномами тоже не можем. Но приближенно и с небольшой точностью — вполне.
Синус
Задача: найти такое приближение синуса на интервале , которое удовлетворительно бы описывало собственно функцию sin(x) для использования при рассчетах поворотов игровых объектов.
Зачем это делать?
Тот синус, который на x86 считается одной инструкцией сопроцессора FSIN, на самом деле довольно точный и, соответственно, медленный. Мы можем определить такой полином, который на некотором интервале вел бы себя как синус, работая для нас его удовлетворительным приближением, а считался бы быстрее, чем FSIN.
Для того, чтобы получить даже не просто синус, а весь набор тригонометрических функций, самое меньшее, что нам надо — что-то похожее на синус на . Потом из этого можно формулами приведения сделать весь косинус, тангенс и котангенс на всей числовой прямой.
Самый простой и общеизвестный способ, которым можно это сделать – интерполяция полиномом Лагранжа по сетке Чебышева. Такое решение дает наилучшую точность на всем интервале, но не обеспечивает гладкости функции при приведении. Это хороший способ приближения функций, но про него и так достаточно написано, поэтому рассмотрим, как это еще можно сделать по-другому.
Что мы знаем про синус. В нуле он равен нулю. В — единице. В некоторой малой окрестности нуля он хорошо приближается прямой: y=x. Единица — максимальное значение синуса. Вот из этого фактажа мы уже можем вылепить себе подходящую систему уравнений.
То, что график функции проходит через точку (x,y), означает, что F(x)=y. Для полиномов, естественно, это тоже справедливо. То есть искомый полином должен удовлетворять двум условиям:
То, что синус в окрестности нуля ведет себя как y=x, означает, что производная его в нуле равна единице. То что в гладкая функция достигает своего максимума, а синус, конечно же, гладкая функция, — что в ее производная равна нулю. То есть:
Итого: четыре уравнения. Мы можем собрать из них систему и решить ее. Это будет система с четыремя неизвестными, соответственно, минимальная степень полинома, который мы так можем получить — 3. . Запишем систему с нашими числами:
В принципе, это уже можно решать. Но тут уже видно, что два уравнения фактичеси решены. Правильно. Если график полинома проходит через (0,0), то его a0 = 0. Обратное верно. Ведь любой член, кроме нулевого, в нуле сам умножается на нулевой икс и превращается в ноль.
То же с производной в нуле. Все члены производной кроме a1 — нули. Стало быть, если производная — единица, то и a1 — единица.
То есть, систему, на самом деле, надо решать не с четырьмя неизвестными, а всего с двумя. Вот ее численное решение:
a3 = -0.110739
a2 = -0.057385
Ну и да:
a1 = 1
a0 = 0
Полином можно записать по схеме Горнера, чтобы его было легче вычислять, и на x86 он действительно посчитается быстрее настоящего синуса. Хотя и не очень точно. Точность обеспечится только до третьего знака, что для функции, не превышающей единицу, довольно грубо. Попробуем добавить еще одно свойство.
Мы знаем, что синус антисимметричен. Если мы хотим определить антисимметричный полином, то надо просто избавить его от всех четных членов. Четные члены всегда симметричны, нечетные — антисимметричны. Антисимметричные функции всегда проходят через ноль, если они в принципе там определены. То есть нам надо определить полином пятой степени, но без четных членов.
То есть:
Решение:
a5 = 0.007403
a3 = -0.165538
a1 = 1
Можно еще попробовать поднять степень и обозначить полином в . Так, чтобы и формулы приведения упростилить. Но тогда придется, ради гладкости, уточнять и производную. А это означает — полином девятой степени. Уже можно натолкнуться на феномен Рунге, хотя сейчас речь даже не о нем. Мы же хотим, чтобы функция была заметно быстрее синуса сопроцессора. Полином девятой степени уже, увы, считается едва ли не медленнее.
На всякий случай же стоит отметить, что так можно строить не только синусы-косинусы. Можно придумать себе любую функцию, определить в каких-то точках ее значение и сколько угодно производных. Для моделирования каких-то простых процессов, движения в какой-нибудь анимации, например, такая интерполяция — очень удобный инструмент.
Анимации на плоскости как раз и посвящен следующий пример.
Полет игрового объекта
Задача. Есть траектория, заданная точками. Игровой объект должен пройти по заданной траектории так, чтобы выйти из исходной точки и прийти в конечную под заранее заданными углами.
Траектория – это кривая. В двухменом пространстве кривую можно задать парой вещественночисловых функций от некоторой переменной. Для удобства назовем эту переменную t – чтобы она была похожа на время, ведь у нас все таки речь идет о полете. Координаты точки в любой момент «времени» опишем так: x(t), y(t). Такой способ задания кривой называется параметрическим.
Очевидно, что мы вполне можем записать известные точки траектории в виде табличной функции и получить гладкие x(t), y(t) интерполяцией. В том числе и полиомиальной.
Но что делать с началом и концом траектории? Нам же надо, чтобы она начиналась и заканчивалась под заданным углом. Наклон касательной графика задается ее производной. В применении к параметрической кривой – парой производных, а точнее даже их отношением .
Вот это отношение само по себе мы не можем, конечно, подставить в систему линейных уравнений. Ведь икс и игрек считаются отдельно. Но нам же нужно задать только отношение. То есть мы можем просто принять, скажем, y'(t) равным единице и посчитать для него x'(t). Вот это уже можно смело подставлять.
Остался невыясненным один вопрос: как поворачивать объект во время полета. Но ответ на него очевиден. Если мы можем использовать угол и соответствующие производные для составления системы уравнений, то и сам полнином для получения угла использовать можно и нужно. Снова таки: , ну или $a = atan({x'(t)}/{y'(t)})$ за исключением двух случаев, когда y'(t) = 0.
Вообще, такой способ позволяет задать траекторию только для небольшого числа точек каркаса. При построении кривых из большого числа заданных точек приходится пользоваться сплайнами. В том числе и полиномиальными.
Но есть одна задача, с которой сплайны тоже справляются не очень хорошо. И она тоже связана с полетами и траекториями. Рассмотрим ее в следующей главе.
Поход персонажа
Задача. Есть игровой объект, которым управляет пользователь — «персонаж». Надо смоделировать его гладкое перемешение по сцене, так, чтобы пользователь мог в любой момент выбрать для него новую точку отправления.
Очевидно, что нужно использовать какой-то локальный сплайн. Так как траектория должна быть гладкой, но больше никаких требований к ней не предъявлено, то минимальная степень такого сплайна – 2.
Итак. Допустим, объект начал двигаться. Сначала мы о нем знаем только точку отправления и точку назначения на первом промежутке. Ну, этого маловато даже для построения определенной параболы. Поэтому первый свой промежуток объект пройдет по прямой. Однако, уже на следующем промежутке мы знаем не просто две точки, но и направление, по которому только что двигался объект. Условие гладкости говорит, что направление движения объекта на новом промежутке должно начинаться таким же. Что на первый взгляд означает равенство производных функции движения объекта в начальной точке второго промежутке и конечной первого.
Однако, если попробовать так и решить задачу, мы получим, конечно, искомую параболу, однако от точки к точке ее траектория будет все более непредсказуемой. Старые производные слишком сильно влияют на то, куда повернет кривая. Надо бы с этим что-то сделать.
На самом деле, для гладкости параметрической прямой не обязательно равенство производных в точке, достаточно равенства отношения этих самых производных. Что означает, мы можем свободно делить обе производные, оставшиеся нам в наследство от первого промежутка на любое удобное число перед тем, как использовать для составления параболы на втором. Графики функций x(t) и y(t) при этом гладкость, конечно же, потеряют, а параметрическая кривая — нет.
Естественно, такой сплайн не обязательно должен быть квадратичным. Теми же соображениями можно пользоваться для построения сплайна любой степени.
Выводы
Конечно, каждую из описанных задач можно решить как-нибудь иначе. Первую, например, можно вообще не решать. Ее постановка была актуальна лет пять назад, вполне вероятно, что на современных процессорах никакого смысла облегчать синус и вовсе нет. Вторую – полиномом Бернштейна нужной степени. Третью – как настоящую кинематическую задачу – решением дифференциального уравнения второго порядка. Описываемая в статье техника – это просто расширение хорошо известной полиномиальной интерполяции, никаких серебряных пуль. Тем не менее, базовые принципы ее не сложны, а лишний способ решения задач в голове лишним не будет.
Автор: akalenuk
