velikol.ru
  1 ... 22 23 24 25
^

Тема 11. Полиномиальная аппроксимация.




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

Я несколько лет работал в фирме Intel в составе группы, составлявшей библиотеки программ. Мне, в частности, пришлось заниматься такими программами, как exp(x), Г(x), lГ(x), cPower(z), функция Бесселя. Библиотека обычно включала функции форматов single, double, long double. С математической точки зрения работа заключалась в подборе одного или нескольких отрезков аппроксимации, определении минимальной допустимой степени полинома на выбранном отрезке, нахождении коэффициентов полинома. В качестве инструментария использовалась система Maple. Работа отнимала много времени, метод Ремеза для нахождения полинома наилучшего работал медленно, иногда давал отказы или зависал.

В какой-то момент я предложил создать свой инструментарий и составил демонстрационную версию программы. К сожалению, в напряженном плане работ Intel не нашлось времени для реализации моего проекта. А жаль, это сэкономило бы много времени в будущем. Здесь я привожу описание этой демоверсии. Основные идеи, заложенные в ней, состоят в следующем:



  1. Использование интерактивного режима, что позволяет часто получить нужную точность без построения полиномов наилучшего приближения,

  2. Использование на стадии вычисления представления полинома в форме Лагранжа, что позволяет использовать эффективный итерационный процесс нахождения полинома наилучшего приближения. Этот процесс в интерактивном режиме можно прервать в любой момент, когда погрешность достаточно мала, но еще не является равномерной.



В демоверсии аппроксимация строится для формата single. Определение степени полинома и вычисление его коэффициентов занимает несколько минут. Подобный инструментарий может быть полезен также экспериментаторам и другим специалистам, работающим с приближенными функциями. Несколько раньше я уже делал программу для построения дробно – полиномиальных аппроксимационных выражений уравнений состояния металлов и других конденсированных веществ.


AppLab




Лаборатория AppLab предназначена для аппроксимации моделируемой функции Fm с помощью полиномиальной функции Fa, зависящей от ряда параметров. Аппроксимация производится на заданном интервале аргумента функции, минимизируется относительная погрешность. Величина погрешности выдается в единицах ULP, что соответствует последнему двоичному разряду мантиссы, для формата single – это примерно 10-7.

Интерфейс и функциональное наполнение комплекса ориентированы на быстрый перебор ряда возможных вариантов выбора Fa с отсеиванием неперспективных случаев на ранних стадиях, что обеспечивает в итоге поиск приемлемого варианта с минимальными затратами времени разработчика. При этом также имеется в виду, что вариант, обеспечивающий к примеру точность 0.1 ULP, может оказаться приемлемым и затраты времени на дальнейшую оптимизацию, например до 0.08 ULP, не оправданы. Важным фактором ускорения принятия решений является автоматическая графическая поддержка всех промежуточных операций.

Аппроксимирующая функция Fa ищется в виде


Fa = F0 + Pr(X) * Pa(X-Xp).


Здесь F0 - константа, Pr - префиксная функция, Pa - полином заданной степени с варьируемыми коэффициентами. Чаще всего префиксные функции не нужны, поэтому по умолчании полагается F0=0, Pr=1. Иногда бывает необходимо выделить некоторую точку X0. Для обеспечения более точной аппроксимации в окрестности этой точки достаточно взять F0=Fm(X0), Pr(X)=X-X0. AppLab предусматривает задание двух названных префиксов, их задание осуществляется достаточно просто и обеспечивает основные потребности. В исключительных случаях может потребоваться задание более сложного префикса, что обеспечивается функцией Pr(X), ее ввод, если она отсутствует в поставляемой библиотеке, требует программирования, о чем будет сказано ниже. В качестве примера можно указать функцию erfc(X), ее аппроксимация для больших значений X значительно облегчается, если взять Pr(X) = exp(-X^2)/X.


^

Инструкция пользователю




В распоряжении пользователя имеется три набора функций: Fm, Pr и Fa. Выбор осуществляется с помощью раскрывающихся списков. При этом Fm, как правило, выбирается из числа имеющихся, Pr требуется, как отмечалось выше, крайне редко и тоже выбирается из списка. Fa, напротив, обычно строится заново, но может быть сохранено в списке. Число функций в наборе ограничено (16 в данной реализации), но все три набора можно сохранять в файле (меню File - Save As…), число файлов не ограничено.

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

  1. Вызвать проект AppLab, в модуле Fm вызвать функцию Fm и вставить в нее под строкой Case S строку Fm = G(X). Здесь в качестве G(X) программируется требуемая функция. Если тело фукции достаточно велико, его можно оформить в виде отдельной процедуры в том же модуле, а в качестве G(X) указать имя этой процедуры,

  2. В процедуре LetArrFm под тем же номером S задать имя функции,

  3. Создать EXE - файл.

Аналогично редактируется или создается новый префикс, функция записывается в процедуру Pr, имя - в LetArrPr модуля AppLab.


Пример: допустим, что в стандартном наборе функций, поставляемом с AppLab нет гиперболического синуса, но есть функция Expm1(X) = Exp(X) - 1. Тогда Вам достаточно задать имя, например, Sinh и запомнить его в свободной строке, допустим 12 - ой. При этом в функции Fm Вам нужно ввести строки (первая из них имеется, вторую нужно изменить):


Case 12

Fm = 0.5 * (Expm1(X) - Expm1( - X))


Заметим, что экономичность данной функции Вас не должна беспокоить, но она должна обеспечивать требуемую точность ( чего, кстати, не дает то же выражение с обычной экспонентой).

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

В диалоге можно задать имя функции, если ее предполагается запоминать. Затем определяются границы интервала (Xleft, Xright). Щелчком на соответствующем идентификаторе точку можно пометить - при этом круглая скобка замениться квадратной. Повторный щелчок отменяет метку. Помеченные точки в режиме Lagr (об этом ниже) считаются точнее, что позволяет 'склеивать' аппроксимацию на соседних отрезках, однако средняя точность на отрезке при этом несколько снижается.

Необходимо также задать точки X0 и Xp. Полная погрешность Fa складывается из погрешности аппроксимации, которая собственно минимизируется в лаборатории, и погрешности округлений, которая определяется величиной коэффициентов полинома Pa(X-Xp) и величиной степеней аргумента. Для уменьшения погрешности желательно брать Xp ближе к интервалу или внутри, например, посередине. Это не всегда удобно для последующего использования, но в ряде случаев является единственным способом избавиться от неприемлемой погрешности округлений.

Точку X0 также можно пометить, при этом перед идентификатором появляется звездочка. Известно, что функцию, имеющую корень на заданном интервале, лучше аппроксимировать с префиксом, например, Sin(X) = X * P(X). Вместо этого можно задать в качестве X0 корень ( для синуса X0 = 0) и пометить эту точку. В режиме Lagr это обеспечит точное совпадение Fa =Fm в этой точке и позволит обойтись без префикса.

Затем задается число параметров Npar. Полином Pa(X-Xp) может быть задан абсциссами Npar лагранжевых точек (ординаты в этом случае определяются функцией Fm и префиксами) или Npar коэффициентами. Степень его в любом случае равна Npar-1. После этого выбором из списка задается способ построения полинома. Первые два способа определяют лагранжевы точки ( равномерно или в корнях полинома Чебышева). Помеченные точки включаются в число этих точек. В двух других случаях находятся коэффициенты полинома путем разложения Fm в ряд по полиномам Лежандра или Чебышева, помеченные точки при этом не учитываются. Преимущество такого способа в том, что построенные полиномы зачастую имеют удовлетворительную точность и не требуют дальнейших итераций. Завершается построение нажатием клавиши OK.

Сохранить все введенные параметры можно нажав клавишу Save и указав место в открывшемся списке. Это можно сделать и позднее, после оптимизации параметров.

Сразу после построения начального приближения Fa на графике строятся обе функции Fm и Fa, соответсвенно красным и синим цветом. Если они на этой стадии заметно отличаются, то скорее всего была допущена ошибка и нужно проверить все параметры или уменьшить интервал (либо увеличить Npar). В случае визуального совпадения графиков следует построить график относительной погрешности и в завимости от величины этой погрешности либо приступать к оптимизации, либо снова вернуться к параметрам Fa. При желании можно посмотреть и абсолютную погрешность.

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

На основной панели располжено также окно отображения полинома Pa. Режим отображения (Lagr, Coef, Fact) определяет также текущее представление полинома, с которым работают оптимизирующие программы и для которого строится график. Это позволяет, в частности, контролировать потерю точности, которая в некоторых случаях возможна при преобразованиях полинома. Преобразования осуществляются в одном направлении: Lagr à Coef à Fact и выполняются кнопками со стрелками. Это соответствует логике действий пользователя: вначале полином удобнее задать в лагранжевой форме, для которой существуют более надежные и разнообразные методы оптимизации, затем перевести в ту форму, которая нужна для использования, подкорректировать, если нужно, точность и записать в файл.

Оптимизация полинома осуществляется вручную или с помощью командных кнопок, расположенных под графическим окном. Прежде за счет выбора параметров (длины интервала, степени полинома, вида префиксов) необходимо, конечно, добиться, чтобы погрешность превышала требуемую величину не более, чем на один - два порядка. При выдаче погрешности лагранжевого полинома абсциссы его точек выделены зеленым или - для помеченных точек - красным цветом. Зеленую точку можно захватить мышкой и двигать влево или вправо, при этом график погрешности будет меняться (красные точки не двигаются). Таким способом иногда можно изменить ситуацию, что не всегда по силам регулярным методам. Чаще, конечно, этот прием служит для понимания характера погрешности и выбора новой стратегии. Автоматический сдвиг точек осуществляется некоторым методом, который условно назван методом покоординатного спуска ( Cdescent ) и вызывается командной кнопкой с этим названием. При этом выполняется фиксированное число итераций, после чего выдается новый график погрешности. При необходимости итерации можно повторить. Метод продемонстрировал достаточную надежность ( предположительно более высокую, чем метод Ремеза ), вторым важным преимуществом является возможность выделения 'помеченных' точек, что необходимо для описания сингулярностей и 'склейки' интервалов.

[В дальнейшем предполагается разработать аналог этого метода для полиномов, заданных коэффициентами и в факторизованной форме. Для обеспечения свободы выбора пользователя необходима также реализация нескольких традиционных методов: Валле Пуссена и др.].


^

Программная реализация




Комплекс содержит формы AppLab, Fm, Fa, Pref и модули AppLab, Fm, Fa, Pref, Polynom, Draw.


Пример.

В качестве примера рассмотрим аппроксимацию функции exp(x) на отрезке ( -1, 1 ) полиномом 8 степени. Начальная аппроксимация строится лагранжевым полиномом на равномерной сетке.





Погрешность составляет 5 ULP. Захватим мышкой первую зеленую точку и потянем ее влево.




Погрешность уменьшилась до 1 ULP. Можно и дальше действовать мышкой, но эффективнее использовать метод спуска.

Нажмем на клавишу метода спуска.




Погрешность уменьшилась до 0.15 ULP.

Нажмем еще раз на клавишу метода спуска.




Погрешность стала меньше 0.1 ULP. Полученный полином еще не является полиномом наилучшего приближения, однако для практической аппроксимации такой точности достаточно. Все манипуляции заняли несколько минут. Точность после преобразования полинома к каноническому виду в данном случае не уменьшилась (на последнем графике полином задан коэффициентами).


Литература




  1. Загускин В.Л. Справочник по численным методам решения алгебраических и трансцендентных уравнений. //Под редакцией А.М.Лопшица. –М.: –ФИЗМАТГИЗ, – 1960.

  2. Рыбаков Л.М. Метод последовательного вычисления всех действительных корней уравнения. //Математическое просвещение. – 1961. №6. – с.262–263.

  3. Загускин В.Л. Вычислительная схема для ускоренного вычисления всех корней уравнения. //Математическое просвещение. – 1961. №6. – с.263–265.

  4. В.В. Воеводин. Программа определения корней алгебраического многочлена. Сб. Вычислительные методы и программирование. МГУ, 1962, с.253–265.

  5. Загускин В.Л., Харитонов А.В. Итерационный метод решения задачи об устойчивости. //ЖВМ и МФ. – 1973, т. 3. №2. – с. 361–364.

  6. Уилкинсон Дж.Х. Алгебраическая проблема собственных значений.: Пер. с англ. –М.: Наука, –1970.

  7. Загускин В.Л. Численные методы решения плохо обусловленных задач. –РГУ, –1976.

  8. Курош А.Г. Курс высшей алгебры. –Москва, Ленинград: – ГИТТЛ, –1949.

  9. Окунев Л.Я. Высшая алгебра. –Москва, Ленинград: – ГИТТЛ, –1949.

  10. Демидович Б.П. и др. Численные методы анализа. –М.: –ФИЗМАТГИЗ, – 1962.

  11. Westerfield E.C. A new bound for the zeros of polynomials.// Amer. Math. Monthly., –1933. –v.40. –№ 1. –p.18–23.

  12. Bairstow L. Investigations Relating to the Stability of the Aeroplane.// Rep. Memor. Adv. Comm. Aero. – Lond., –1914. №154. –p.51–64.

  13. Загускин В.Л. Универсальные методы для вычисления корней многочленов и решения задачи об устойчивости.// Доклад на 1 Всесоюзном математическом съезде, секция "Вычислительная математика". – 1961.




<< предыдущая страница