CC BY
10ISSN 0868-5886
НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2023, том 33, № 3, c. 74-91
МАТЕМАТИЧЕСКИЕ МЕТОДЫ И МОДЕЛИРОВАНИЕ В ПРИБОРОСТРОЕНИИ
УДК537.213, 537.612, 517.5 © А. С. Бердников, С. В. Масюкевич, 2023
ЧИСЛЕННЫЙ АЛГОРИТМ ДЛЯ МИНИМАКСНОЙ ПОЛИНОМИАЛЬНОЙ АППРОКСИМАЦИИ ФУНКЦИЙ
С ЗАДАННЫМ ВЕСОМ
В статье рассматривается быстросходящийся численный алгоритм для определения полиномов заданной степени, который обеспечивает на заданном интервале оптимальное приближение заданной функции в минимаксной норме с заданным весом при условии, что весовая функция не обращается в ноль на рассматриваемом интервале, за исключением, быть может, начальной и/или конечной точек интервала.
Кл. сл.: минимаксная норма, полиномы Чебышева, оптимальная аппроксимация, интерполяция
ВВЕДЕНИЕ
Применение аппроксимаций функций, оптимальных по минимаксной (равномерной) норме, имеет существенные преимущества по сравнению с более простой аппроксимацией функций по методу наименьших квадратов [1, 2]. В качестве полезного инструмента при конструировании таких аппроксимаций используются полиномы Чебышева первого рода, наименее отклоняющиеся от нуля с весовой функцией, равной единице [3, 4]. Разложение функции в усеченный ряд, состоящий из полиномов, наименее отклоняющихся от нуля (в случае единичной весовой функции — полиномов Чебышева), является одним из способов приближенного конструирования полиномиальных аппроксимаций, оптимальных по минимаксной норме [4, 5]. Также для конструирования полиномиальных аппроксимаций, которые являются хорошим приближением к оптимальной минимаксной аппроксимации, используют аппроксимацию по дискретному набору коллокационных точек, соответствующих нулям полинома Чебышева со степенью на единицу больше, чем степень аппроксимирующего полинома [1, 2, 4, 5].
Такие способы характеризуются удобством и легкостью конструирования, а также определенными эмпирическими предпосылками получить на выходе достаточно хорошую точность. Однако упрощенные способы аппроксимации отличаются от истинной аппроксимации, оптимальной по минимаксной норме, и, следовательно, обеспечивают неконтролируемо большую ошибку аппроксимации, чем точное решение соответствующей минимаксной оптимизационной задачи. Использование аппроксимаций, обеспечивающих для заданной функции точное решение минимаксной оптимиза-
ционной задачи (когда это возможно), определенно более предпочтительно.
Множество функций, для которых имеется явное алгебраическое (аналитическое) представление для полиномиальной аппроксимации, являющейся оптимальной по минимаксной норме, не слишком велико (практически все такие случаи приводятся в [6-8]). В данной работе исследуется уточненный вариант быстро сходящегося численного алгоритма для определения коэффициентов полиномов, наименее отклоняющихся от заданной функции на заданном интервале с заданным весом, который ранее был рассмотрен в [9].
МИНИМАКСНАЯ АППРОКСИМАЦИЯ ФУНКЦИЙ И ТЕОРИЯ ЧЕБЫШЕВА
Постановка задачи
Пусть имеется непрерывная функция Дх) и конечный интервал [а, Ь], для которого задана непрерывная весовая функция ю(х), являющаяся строго положительной на этом интервале, за исключением, быть может, концов интервала, для которых она может обращаться в ноль.1 Требуется найти полином р(х) степени не выше п с заранее не-
1 Для весовой функции с нулями внутри интервала [а, Ь] потребуется, чтобы максимумы и минимумы Чебышева (и пробные точки для численного алгоритма) не совпадали с нулями весовой функции, а знаки попеременных максимумов и минимумов изменялись в соответствии со знаком весовой функции. Также интервал [а, Ь] может быть бесконечным справа и/или слева, но при этом и аппроксимируемая функция Дх), и весовая функция q(x) должны стремиться к нулю на бесконечности не медленнее, чем степенные функции \/хк, причем с одинаковой скоростью [6-8].
определенными коэффициентами a0, a1, a2, ..., an:
p(x) = a0 + a1x + a2x2 + ... + an-1xn-1 + anxn, (1)
который является решением оптимизационной задачи
max |w(x) (fx) - p(x))| ^ min, (2)
где максимизация выполняется по переменной x е [a, b], а минимизация выполняется по коэффициентам a0, a1, a2, ..., an.
Без ограничения общности можно считать, что функция fx) не является полиномом степени n или ниже, — в противном случае оптимальное решение для задачи (2) заранее известно и совпадает с функцией fx).
Критерий Чебышева для минимаксной полиномиальной аппроксимации
Для того чтобы полином p(x) степени n был решением задачи (2), необходимо и достаточно, чтобы существовали такой набор из n + 2 точек x0 < x1< x2 < ... < xn + 1, принадлежащих интервалу [a, b], и такое число в (положительное либо отрицательное), чтобы были выполнены условия:
—|в| < w(x) (f(x) —p(x)) < |в| при x е [a, b], (3) ®(x,) (f(xk) — p(xk)) = (—1)k в при к = 0, 1, .., n + 1. (4)
Точки Х0 < Х1 < х2 < ... < хп + 1 с чередующимися положительными минимумами и отрицательными максимумами, равными по абсолютной величине максимуму модуля исследуемой функции, называются альтернансом Чебышева (см. [7, 8, 10-14]). Графики полиномов Чебышева первого рода [3], демонстрирующие выполнение данного условия на интервале [-1, +1] с весом ю(х) = 1, приводятся на рис. 1.
Смысл условий (3), (4) состоит в том, что для полинома р(х) имеется набор из п + 2 пробных точек хк, в которых невязка ю(х) (Дх) - р(х)) имеет локальные чередующиеся отрицательные минимумы и положительные максимумы, равные ±|в|, причем значения этих минимумов и максимумов являются глобальными на интервале [а, Ь]. Этот критерий является частным случаем теории Че-бышева минимаксной аппроксимации функций с помощью дробно-рациональных функций [7, 8], однако в случае аппроксимации с помощью полиномов доказательство соответствующих утверждений упрощается.
Утверждение 1 (теорема Валле - Пуссена [7, 8, 15-17]). Если рассматриваемая функция ю(х) (Дх) -- р(х)) для некоторого многочлена р(х) степени п принимает в N последовательных точках х0 < х1 < ... < xN _ 1 интервала [а, Ь] отличные от нуля значения V ..., ^ - 1 с чередующимися
Рис. 1. Полиномы Чебышева Тп(х), которые наименее отклоняются от нуля на отрезке [-1, +1]. а — п = 5, б — п = 8, в — п = 16
знаками, где N > n + 2, то для любого другого многочлена r(x) степени n справедливо
max |o>(x) (fix) - r(x))| > min ..., - i|}.
Доказательство. Пусть имеется многочлен r(x), для которого выполняется условие max |ш(х) (fix) - r(x))| < min (|^о|, ..., \kN - i|). Будем считать, что в последовательности xk для четных точек значения ®(xk) (fixk) - r(xk)) строго положительные, а для нечетных точек строго отрицательные. (Когда для нечетных точек значения положительные, а для четных точек отрицательные, рассуждения аналогичны.) Это означает, что в четных точках xk выполнено условие ®(xk) (fixk) -- r(xk)) < |Ak| = ®(xk) (fixk) - p(xk)), а в нечетных точках xk выполнено условие ®(xk) (f(xk) - r(xk)) > > -|Xk| = ®(xk) (fixk) - r(xk)). В результате величина ®(x) (p(x) - r(x)) строго отрицательна в четных точках xk и строго положительна в нечетных точках xk. Значения весовой функции ®(xk) в точках xk, для которых функция ®(x) (fix) - r(x)) принимает строго положительные или строго отрицательные значения, не обращаются в ноль и, следовательно, строго положительны. Следовательно, не равный тождественно нулю многочлен p(x) - r(x) степени n попеременно принимает строго положительные и строго отрицательные значения как минимум в n + 2 разных точках и в силу этого имеет не менее n + 1 нулей кратности 1. Значит, многочлен r(x), для которого выполняется условие max |®(xk) (fixk) - r(xk))| < min (|^o|, ..., |V - i|), не существует.
Замечание. Теорема Валле - Пуссена позволяет получить для решения оптимизационной задачи (2) оценку снизу. Если Х0, Х1, ..., V _ 1 — это локальные максимумы и минимумы функции ®(x) (fix) - r(x)) с чередующимися знаками, то такая конструкция называется альтернансом Валле -Пуссена, и для решения оптимизационной задачи (2) она дает не только оценку снизу, но и оценку сверху, равную max ..., IV1I}.
Утверждение 2. Если для полинома p(x) выполнены условия (3), (4) критерия Чебышева, то оптимальное значение для правой части оптимизационной задачи (2) равно |s|, а полином p(x) является ее решением (возможно, одним из многих).
Доказательство. Альтернанс Чебышева (3), (4) является частным случаем альтернанса Валле -Пуссена, поэтому, согласно Утверждению 1, правая часть оптимизационной задачи (2) не меньше |s|. C другой стороны, max|®(x) (fx) - p(x))| = |s|. Следовательно, минимум оптимизационной задачи (2) равен |s|, а полином p(x) является одним из ее решений.
Утверждение 3. Существует полином p(x), который обеспечивает решение оптимизационной задачи (2). (Это исключает вариант, когда есть полиномы pk(x) со все более и более уменьшающимися значениями для Pk = max |®(x) (fix) - pk(x))|, в то время как минимум этой величины никогда не достигается на множестве полиномов заданной степени.)
Доказательство. Значения величин Pk = = max |®(x) (fx) - pk(x))| ограничены снизу нулем, поэтому для значений Pk на множестве полиномов существует точная нижняя грань P. В соответствии с определением точной нижней грани существует последовательность полиномов pk(x) степени n:
Pk(x) = a0, k + a1, kx + a2, kx + . + an - 1, k x" 1+an, k x",
для которой lim Pk = P при k ^ да. Следовательно, начиная с некоторого номера k значения многочленов pk(x) на интервале [a + 5, b - 5], где 5 достаточно мало, будут ограничены сверху и снизу (от концов интервала пришлось сделать отступ, чтобы учесть случай, когда весовая функция на концах интервала может быть равна нулю). Из формулы Лагранжа [1, 2] для полинома степени n, который в n + 1 заданных точках xk принимает значения yk
f
p (x )= I yj п
j=1,n+1
/=1, n+1,/ Ф j
(x - x)
(xj - x)
л
следует, что когда в n + 1 фиксированных точках xk значения полинома yk ограничены сверху и снизу, то каждый из коэффициентов полинома также ограничен сверху и снизу. Следовательно, в соответствии с теоремой Больцано - Вейерштрасса (иначе называемой леммой Больцано - Вейерштрасса) о предельной точке, согласно которой из всякой бесконечной ограниченной последовательности точек пространства Rn можно выделить сходящуюся подпоследовательность, в последовательности полиномов pk(x) можно выбрать подпоследовательность полиномов, в которой для каждого из коэффициентов полинома имеется предел. То есть имеется такая последовательность полиномов pk(x), для которой lim Pk = P при k ^ да, а у коэффициентов полиномов имеются предельные значения bj = lim a;-, k при k ^ да.
Рассмотрим полином r(x) с коэффициентами bf.
r(x) = b0 + b1x + b2x2 + ... + bn _ 1xn _ 1 + bnxn.
Поскольку bj = lim ajk, то для любого фиксированного значения x выполнится условие lim pk(x) = r(x) при k ^ да, причем стремление к пределу является равномерным на рассматри-
ваемом конечном интервале значений x. Из неравенства
max |w(x) (fix) - r(x))| - max |w(x) (fix) - pk(x))| + + max |w(x)| max |r(x) - pk(x)|
следует, что max |w(x) (f(x) - r(x))| = P, т.к. меньше, чем P, эта величина быть не может, в то время как при k ^ да первое слагаемое в правой части стремится к P, а второе слагаемое стремится к нулю. Это означает, что для полинома r(x) величина (2) достигает своей нижней границы.
Утверждение 4 (теорема Чебышева). Полином p(x), обеспечивающий решение оптимизационной задачи (2), удовлетворяет критерию Чебышева.
Доказательство. Рассмотрим поведение функции F(x) = w(x) (fx) - p(x)) на отрезке [a, b]. Для упрощения рассуждений будем считать, что число интервалов, на которых функция F(x) меняет знак, конечно. Это не означает, что число нулей функции F(x) также конечно — допустим вариант, когда функция F(x) равна нулю в любой точке некоторого интервала внутри отрезка [a, b]. Пусть точки y1, y2, ..., ym разбивают отрезок [a, b] на m + 1 интервалов, на каждом из которых функция F(x) принимает попеременно то положительное, то отрицательное значение, причем допустим вариант, когда на всем отрезке [a, b] функция F(x) принимает только положительные или только отрицательные значения. Если имеется непрерывный ин-
тервал, целиком состоящий из нулей функции F(x), этот интервал присоединяется к соседнему интервалу с положительными либо отрицательными значениями функции.
Для интервалов с положительными значениями F(x) выбираем точку с максимальным значением функции на этом интервале (возможно, не единственную на этом интервале), а для интервалов с отрицательными значениями F(x) выбираем точку с минимальным значением функции на этом интервале. Получаем набор точек y1, y2, ..., ym, — нулей функции F(x), разбивающих отрезок [a, b] на m + 1 интервалов, для которых функция F(x) сохраняет свой знак, изменяя его при переходе через границу интервала, и набор принадлежащих этим интервалам точек x0, x1, x2, ..., xm, xm + 1, не совпадающих с концами интервалов, который состоит из чередующихся положительных максимумов и отрицательных минимумов Xo, Xb X2, ..., Xm + 1 функции F(x) (альтернанс Валле - Пуссена и прилагающееся к нему разбиение отрезка [a, b] на интервалы, на которых функция F(x) сохраняет знак).
Пусть X = max |Xk|. Если на каком-то интервале максимум равен X, а на соседнем (последующем или предшествующем) интервале модуль минимума меньше X, то этот интервал вместе с последующим интервалом с положительным максимумом, независимо от его значения, присоединяется к текущему интервалу (см. рис. 2, а).
Рис. 2. Объединение отрезков с чередующимися максимумами и минимумами.
а — случай одиночного минимума или нескольких последовательно расположенных минимумов, которые лежат выше минимального уровня функции F(x); б — случай одиночного максимума или нескольких последовательно расположенных максимумов, которые лежат ниже максимального уровня функции F(x). Вертикальными стрелками отмечены минимумы и максимумы, которые соответствуют максимальному отклонению. Наклонными стрелками отмечены минимумы и максимумы, которые по модулю меньше максимального отклонения, так что график функции можно сдвинуть вниз или соответственно вверх на ненулевое расстояние, не выходя за границы максимального отклонения. Жирными точками отмечены нули функции F(x), являющиеся границей комбинированного интервала со строго положительными либо со строго отрицательными максимальными отклонениями функции F(x)
Новый интервал обладает тем свойством, что из функции F(x) можно вычесть такую положительную константу, чтобы положительный максимум функции на интервале уменьшился, но при этом абсолютное значение отрицательного минимума функции на интервале не настолько увеличилось, чтобы превысить новый максимум функции.
Точно так же если на каком-то интервале минимум равен -X, а на соседнем интервале максимум меньше X, то этот интервал вместе с последующим интервалом, содержащим отрицательный минимум, присоединяется к текущему интервалу (см. рис. 2, б). В этом случае на полученном интервале к функции F(x) можно прибавить такую положительную константу, чтобы уменьшилось абсолютное значение отрицательного минимума на интервале, но при этом не настолько увеличился положительный максимум функции на интервале, чтобы он стал больше абсолютного значения нового минимума.
В конечном итоге на интервале [a, b] для заданного полинома p(x) будет получен набор из N интервалов и принадлежащих им точек x0, x1, x2, ..., xN_ 1 с чередующимися положительными максимумами и отрицательными минимумами функции F(x), равными ±А, где А = max |®(x) (fx) _ _p(x))|. Если N больше или равно n + 2, то такой полином p(x) удовлетворяет критерию Чебышева и, в соответствии с Утверждением 1, будет решением оптимизационной задачи (2).
Пусть N меньше или равно n + 1. На этапе объединения интервалов был получен набор точек y1, y2, ..., yN_ 1 (границы интервалов), в которых функция F(x) обращается в ноль и меняет свой знак. Эти точки разбивают интервал [a, b] на интервалы меньшего размера, на каждом из которых имеется точка с чебышевским минимумом или чебышевским максимумом. Также имеются константы 5k > 0, которые можно вычитать из функции F(x) на интервалах с положительными максимумами или прибавлять к функции F(x) на интервалах с отрицательными минимумами, уменьшая при этом минимаксную норму на соответствующем интервале.
Без ограничения общности можно считать, что на первом интервале функция F(x) принимает положительные значения, — когда это не так, можно изменить знак одновременно у функции fx) и у полинома p(x), не меняя оптимизационную задачу (2).
Рассмотрим функцию ^(x), которая представляет собой произведение полинома степени не выше n и весовой функции ®(x):
Q(x) = a>(x) (у1 _ x) (у2 _ x)...(yN_ 1 _ x). Функция Q(x) строго положительна на тех интер-
валах, на которых имеется чебышевский максимум функции F(x), и строго отрицательна на тех интервалах, на которых имеется чебышевский минимум функции F(x). Если вычесть из функции F(x) функцию sQ(x) с достаточно малым положительным множителем s (например, можно выбрать s = (5 / R), где R = max |fi(x)|, а 5 = min 5k > 0), это безопасным образом уменьшит положительные максимумы функции F(x) и уменьшит отрицательные минимумы функции F(x), тем самым уменьшив величину max |F(x)|. Следовательно, такой полином p(x) не может быть решением для оптимизационной задачи (2), поскольку можно получить еще меньшее значение для max |F(x)| при замене текущего полинома p(x) степени n на новый полином r(x) степени n:
G(x) = F(x) - (5 / R) Q(x) =
= a>(x) (fx) -p(x)) - (5 / R) a>(x) (y - x)...(yN- 1 - x) =
= ®(x) (fx) -p(x) - (5 / R) (y - x)...(yN- 1 - x)) =
= ф) (f(x) - r(x)), max |G(x)| < max |F(x)|.
Вывод: оптимальный полином степени n, рассмотренный в Утверждении 2, обязан удовлетворять критерию Чебышева.
Утверждение 5. Полином p(x), удовлетворяющий критерию Чебышева и обеспечивающий решение оптимизационной задачи (2), является единственным.
Доказательство. Идея доказательства заимствована из [7, 8]. Пусть имеются два полинома p(x) и r(x) степени n, которые являются решением оптимизационной задачи (2). Для вырожденного случая, для которого max |®(x) (fx) - p(x))| = 0 и max |®(x) (f(x) - r(x))| = 0, должны выполняться условия fx) - p(x) = 0 и fx) - r(x) = 0, — что означает, что p(x) = r(x). Поэтому будем считать, что
max |®(x) (f(x) - p(x))| = max |®(x) (fx) - r(x))| > 0,
так что у выражений ®(x) (fx) - p(x)) и ®(x) fx) -- r(x)) на интервале [a, b] имеются не равные нулю максимумы и минимумы.
В соответствии с Утверждением 3, каждый из полиномов p(x) и r(x) удовлетворяет критерию Чебышева, причем значение s ^ 0 для них одно и то же. Пусть x1, x2, ..., xn + 2 — это чебышевские точки уклонения с чередующимися максимумами и минимумами, равными ±s, для функции F(x) = ®(x) (f(x) - p(x)). Если для полинома r(x) функция G(x) = ®(x) (f(x) - r(x)) принимает в точ-
ках x1, x2, ..., xn + 2 те же самые значения ±s, кроме, может быть, одной точки из этого списка, то полиномы p(x) и r(x) тождественно равны друг другу. Действительно, в этом случае разность
d(x) = Fx) - G(x) =
= œ(x) (fx) - p(x)) - œ(x) fx) - r(x)) =
= œ(x) (r(x) -p(x)), (5)
которая представляет собой произведение положительной весовой функции œ(x) и полинома степени не выше n, обращается в ноль в n + 2 точках или, в крайнем случае, в n + 1 точке, — и поэтому полином r(x) - p(x) обязан быть тождественным нулем, поскольку иначе не объяснить наличия у него такого количества нулей (весовая функция œ(x), даже если она обращается в ноль на концах интервала [a, b], не равна нулю в точках минимумов и максимумов функции F(x)).
Пусть среди чебышевских точек уклонения x1, x2, xn + 2 для функции F(x) имеются хотя бы две точки xk и xk + m, в которых значение функции G(x) не совпадает со значением функции F(x) в этих точках (минимумом либо максимумом функции F(x)). Требуется доказать, что в таком случае внутри интервала [xk, xk + m] найдется по крайней мере еще одна точка, не совпадающая с точками x1, x2, xn + 2, в которой значения функций F(x) и G(x) совпадают (хотя, возможно, в этой точке общее значение функций F(x) и G(x) и не равно ±s), — а следовательно, у функции d(x), заданной формулой (5), имеется дополнительный неучтенный ноль.
Если в какой-то точке xj выполнено условие d(xj) Ф 0, то знак разности d(xj), заданной формулой (5), совпадает со знаком F(x;). Действительно, если xj это точка отрицательного минимума функции F(x), то G(x;) > F(xj), а если xj это точка положительного максимума функции F(x), то G(x;) < F(x;). Если
d(xk) Ф 0, d(xk + i) = 0, d(xk + 2) = 0, d(xk + m-1) = 0,
d(xk + m) Ф 0, (6)
то выражения (-1)k F(xk) и (-1)k + m F(xk + m) имеют одинаковые знаки и, следовательно, одинаковые знаки имеют выражения d(xk) и (-1)m d(xk + m). Следовательно, четность числа перемен знака функции d(x) на интервале [xk, xk + m] совпадает с четностью числа m, а тем самым и четность числа корней (с учетом их кратности) для функции d(x) на интервале [xk, xk + m] совпадает с четностью числа m. Согласно (6) у функции d(x) на интервале [xk, xk + m] имеется как минимум m - 1 не равных друг другу корней и, следовательно, либо должен существовать еще один корень, либо один из корней xk + 1, xk + 2, ..., xk + m _ 1 имеет четную кратность, как
минимум равную двум. Проделав эту операцию для всех интервалов, ограниченных точками xk с ненулевыми значениями d(xk), получим, что число корней функции d(x) на интервале [a, b] с учетом кратности корней будет не меньше, чем n + 1. Следовательно, в формуле d(x) = œ(x) (r(x) - p(x)) полином r(x) - p(x), имеющий степень не выше n, обязан быть тождественным нулем.
АЛГОРИТМ КОНСТРУИРОВАНИЯ МИНИМАКСНОЙ ПОЛИНОМИАЛЬНОЙ АППРОКСИМАЦИИ ФУНКЦИИ
Из критерия Чебышева для минимаксной аппроксимации функции выводится алгоритм численного построения аппроксимирующих полиномов, который представляет собой измененный и оптимизированный алгоритм Е.Я. Ремеза [1822]:
1. Произвольным образом выбирается начальный набор точек x0 < x1 < x2 < ...< xn + 1, принадлежащих интервалу [a, b]. Если весовая функция œ(x) обращается в ноль в точке x = a, необходимо обеспечить выполнение условия a < x0. Если весовая функция œ(x) обращается в ноль в точке x = b, необходимо обеспечить выполнение условия x0 < b.
2. Для текущего набора точек a < x0 < x1 < < x2 < .. .< xn + 1 < b конструируется вспомогательная функция g(x) = a fx) - p(x), где a это некоторая константа, а p(x) — полином степени n с некоторыми коэффициентами. Свободные коэффициенты функции g(x), т.е. множитель a и коэффициенты полинома p(x), выбираются так, чтобы были выполнены условия œ(xk)g(xk) = (-1)k при k = 0, 1, ..., n + 1. Такая функция существует и определяется единственным образом (см. далее).
3. Теперь надо определить корни алгебраического уравнения g(x) = 0, которые находятся на интервале [a, b]. Поскольку на краях интервалов [xk, xk + 1] значения функции g(x) имеют разные знаки, то внутри каждого такого интервала существует по крайней мере один корень уравнения g(x) = 0. Если функция f(x) меняется достаточно медленно, то на каждом из интервалов [xk, xk + 1] имеется ровно один корень кратности единицы, а численное определение этих корней не представляет трудности. Если же это условие нарушается и на рассматриваемом интервале имеется несколько корней, все эти корни добавляются к списку (кратные корни с четной кратностью, для которых значения функции g(x) по обе стороны от корня имеют один и тот же знак, в список не включаются). По завершении работы по поиску корней к списку значений a < y1 < y2 < ...< ym < b, являющихся корнями нечетной кратности функ-
ции g(x) на интервале [a, b], добавляются начало и конец интервала y0 = a и ym + 1 = b.
4. Теперь интервал [a, b] представляет собой совокупность не менее чем n + 2 интервалов [yk, yk + 1], которые соприкасаются своими конечными точками, причем на каждом интервале функция g(x) сохраняет знак, а при пересечении границы между интервалами знак функции меняется. Рассмотрим функцию G(x) = ю(х) g(x), которая на интервале [a, b] принимает значения того же знака, что и функция g(x). Найдем для функции G(x) на каждом интервале x е [yk, yk+1] глобальный (для рассматриваемого интервала) экстремум Gk и точку zk е [yk, yk + 1], в которой расположен этот экстремум. Для интервалов с положительными значениями функции G(x) надо искать положительный максимум, а для интервалов [yk, yk + 1] с отрицательными значениями функции — отрицательный минимум.
5. Если число найденных точек Gk = G(zk), zk е [yk, yk + 1] с чередующимися знаками больше n + 2, необходимо исключить лишние точки. Для этого выбирается наименьшее по модулю значение Gk. Если интервал с этим значением первый либо последний в списке, значение Gk и соответствующая этому значению точка zk отбрасываются. Если интервал с этим значением находится в середине списка, надо проверить предшествующий и последующий интервалы и отбросить не только рассматриваемое значение Gk, но также и то из двух соседних значений Gk _ 1 и Gk + 1, которое будет меньше по модулю. Действительно, знаки значений Gk _ 1 и Gk + 1 противоположны по знаку значению Gk, а т.к. из-за удаления из списка значения Gk два значения одного знака окажутся рядом, то из двух значений Gk + 1 и Gk_ 1 одно надо отбросить. Важно: поскольку на каждом шаге удаляется одна либо две точки, то есть риск из n + 3 точек zk получить сразу n + 1 точек, вместо требуемых n + 2 точек, — в этом случае следует оставить точку Gk в списке, а вместо нее удалить ту из начальной и конечной точек, которая меньше по модулю. Процесс продолжается, пока в списке не окажется n + 2 экстремумов Gk = G(zk), которым соответствуют n + 2 точек zk с чередующимися положительными максимумами и отрицательными минимумами функции G(x) = ю(х) g(x).
6. Для имеющегося списка точек a < z0 < z1 < < z2 < .. .< zn + 1 < b рассмотрим более внимательно значения экстремумов Gk = ®(zk) g(zk) = (-1)ksk. Если в пределах заданной точности sk « 1, то полином, удовлетворяющий требованиям критерия Чебышева, найден. При этом из теоремы Валле -Пуссена о наилучшем приближении функции полиномами (Утверждение 1) следует, что точный оптимум для задачи (2) лежит в диапазоне между min |sk / a| и max |sk / a| (где a это множитель в вы-
ражении g(x) = a fx) - p(x)). Если же условие sk « 1 не выполнено, заменяем пробные точки a < x0 < < xi < x2 < .. .< xn + i < b на точки a < z0 < zi < z2 < .. .< zn + 1 < b и возвращаемся к шагу 2.
7. Остается нормировать функцию g(x) = = a fx) - p(x) так, чтобы коэффициент при функции fx) стал равен единице. Такая нормировка всегда осуществима, поскольку коэффициент a не может быть нулем, — иначе окажется, что имеется отличный от тождественного нуля полином p(x) степени n, у которого есть n +1 корней, находящихся между соседними максимумами и минимумами разных знаков в количестве n + 2 штук.
Очевидно, что если этот алгоритм сходится, то он сходится к полиному, удовлетворяющему критерию Чебышева, и тем самым результат работы алгоритма является решением оптимизационной задачи (2), причем единственным таким решением. Теоретически алгоритм сходится всегда, причем со скоростью геометрической прогрессии: для любой непрерывной весовой функции и для любой степени n найдутся такие числа C > 0 и 0 < X < 1, для которых
max |rn(x) (afx) - p,-(x))| < 1 + СИ,
где i это номер итерации. Доказательство аналогично доказательству сходимости алгоритма Е.Я. Ремеза в монографии [20]. К сожалению, из-за ошибок округления алгоритм может расходиться для "плохих" функций fx) и ®(x), — а именно для функций, которые не меняются медленно на интервале [a, b], а вместо этого демонстрируют многочисленные осцилляции, так что для них имеется более одного корня для интервалов [xk, xk + 1] или более одного локального максимума или минимума для интервалов [yk, yk + 1]. В таких случаях алгоритм может выполнять бесконечный цикл, фактически переключаясь между одним и тем же списком тестируемых полиномов. Такое поведение алгоритма следует контролировать вручную, задавая максимальное количество циклов, которое разрешено выполнить. Другим полезным критерием для прекращения работы алгоритма, попавшего в бесконечный цикл, может быть проверка возрастания на очередном шаге разброса между максимумами и минимумами ±sk, найденными на интервалах [yk, yk + 1], по сравнению с этой величиной для предыдущего шага. Также подозрительным признаком, сигнализирующим о зацикливании алгоритма, может быть регулярное возникновение ситуации, когда на шаге 5 при "чистке" точек Gk на каком-то шаге в списке оказывается меньше чем n + 2 точки и возникает необходимость сделать шаг назад, удалив первую или последнюю точку, — ведь даже если оптимальным многочленом p(x) является тождест-
венный ноль (например, такая ситуация возникает при попытке аппроксимации на интервале [0, 1] многочленом конечной степени функции sin(1 / х)), то, согласно Утверждению 4, все равно должен существовать набор из п + 2 точек zk, которые формируют альтернанс Чебышева.
Так как весовая функция ю(х) на каком-то из концов интервала [а, Ь] может обращаться в ноль, то для выполнения шага 2 потребуется, чтобы пробные точки, выбираемые на шаге 1, лежали строго внутри интервала [а, Ь]: а < х0 < XI <
< х2 < ...< х„+1 < Ь. Можно убедиться, что если ю(а) = 0 либо ю(Ь) = 0 и при этом а <х0 < х1 < х2 < .< хп+1 < Ь, то условие а < х0 < х1 < х2 < .< хп+1 <
< Ь сохранится для всех последующих итераций, т.к. точки zk на шаге 6 также должны будут удовлетворять условию а < z0 < z! < z2 < ...< zn + 1 < Ь.
Для вычисления коэффициентов на шаге 2 нельзя применять явное решение системы линейных уравнений для нахождения коэффициентов функции £(х) = аДх) - р(х), поскольку обусловленность матрицы линейных уравнений (за исключением первого столбца, совпадающая с матрицей Вандермонда) быстро растет с увеличением степени полинома — в результате малые ошибки округления чисел с плавающей точкой приводят к совершенно неудовлетворительным значениям для неизвестных коэффициентов. Один из работоспособных способов нахождения искомой функции g(x) = аДх) - р(х) на шаге 2 состоит в следующем. Из условий р(хк) = а-Дхк) - (-1)к / ю(хк), которые равносильны условиям ю(хк) (аДхк) --р(хк)) = (-1)к, при к = 0, 1, ..., п по формуле Ла-гранжа2 [1, 2] находятся коэффициенты полинома р(х), который в заданных точках хк принимает заданные значения. Эти коэффициенты будут представлять собой линейные выражения со свободным коэффициентом а. После этого для уравнения ю(хк) (а-Дхк) - р(хк)) = (-1)к при к = п + 1, которое представляет собой линейное уравнение относительно переменной а, находится неизвестный коэффициент а.
Другой, более эффективный способ определения коэффициента а и коэффициентов полинома р(х), для которых на шаге 2 будут выполнены условия
ю(хк) (а-Дхк) -р(хк)) = (-1)к при к = 0, 1, .., п+1,
состоит в следующем. Для выполнения равенств р(хк) = а-Дхк) - (-1)к / ю(хк) по схеме конечных раз-
2 Вместо явной формулы Лагранжа, особенно при больших степенях полиномов, удобно применять итерационную схему Эйткена [2]. Хотя этот алгоритм требует больше памяти, он обеспечивает существенно большую точность коэффициентов и требует меньше компьютерного времени.
ностей Ньютона [1, 2] определяются такой полином г(х), для которого выполнены условия г(хк) = = (-1) / ю(хк) при к = 0, 1, ., п, и полином 5(х), для которого выполнены условия 5(хк) = Дхк) при к = 0, 1, ..., п. Затем из условия
ю(хп + 1) (а-Д(хп + 1) - а^(хп + 1) + г(хп + 1)) = (-1) п + 1
определяется константа а. После этого искомый полином р(х) определяется в соответствии с формулой р(х) = a•s(x) - г(х). Этот способ можно использовать также в следующем виде: вычисляются полином s(x) степени п + 1, для которого в точках хк при к = 0, 1, ., п + 1 выполнены условия s(xk) = Дхк), и полином г(х) степени п + 1, для которого в точках хк при к = 0, 1, ., п + 1 выполнены условия г(хк) = (-1)к / ю(хк). Остается подобрать коэффициент а так, чтобы полином a•s(x) - г(х) имел степень п, — что всегда возможно, если только функция Д(хк) не является полиномом степени не выше п.
Отметим, что на промежуточных этапах нет необходимости определять коэффициенты полинома р(х) в явном виде, достаточно уметь вычислять его значение в произвольной точке. Поэтому при компьютерной реализации данного алгоритма вплоть до получения итогового результата можно использовать вычисление значений полинома р(х) по формуле р(х) = a•s(x) - г(х), где для вспомогательных полиномов s(x) и г(х) используется форма конечных разностей Ньютона.
Для улучшения работы алгоритма на начальном этапе на шаге 1 рекомендуется выбирать для начальных точек а < х0 < х1 < х2 < .. .< хп + 1 < Ь нули полинома Чебышева первого рода степени п + 2, смещенные и масштабированные с интервала [-1, +1] к интервалу [а, Ь]. На шаге 2 вместо условия ю(хк) р(хк) = (-1) рекомендуется использовать условие ю(хк) р(хк) = (-1)к(4 / (Ь - а))-п - 1, чтобы избежать неоправданно больших коэффициентов у полинома р(х) и функции Дх) (такой выбор масштаба соответствует осцилляциям полинома Чебышева степени п + 1, пересчитанного от интервала [-1, +1] к интервалу [а, Ь] и масштабированного так, чтобы старший коэффициент был равен единице).
Пример. Пусть Д(х) = ехр(х), а = 0, Ь = 1, ю(х) = 1. Для п = 1, 2, ., 8 в результате работы описанного алгоритма получаем полиномы, которые обеспечивают оптимальную (по минимаксной норме) аппроксимацию функции ехр(х) на интервале [0, 1]. Ошибка аппроксимации в зависимости от степени полинома указана в табл. Зависимость ошибки аппроксимации от степени полинома в графической форме показана на рис. 3.
Табл. Ошибка аппроксимации функции ехр(х) на интервале [0, 1] многочленами степени п, оптимальными по минимаксной норме с весом д(х) = 1
п 1 2 3 4 5 6 7 8
г 0.1 810-3 5-10-4 3-10-5 10-6 4-10-8 10-9 3-1011
Рис. 3. Аппроксимация функции ехр(х) на отрезке [0, 1] минимаксными многочленами степени п. а — экспоненциальная функция и ее аппроксимации, б — разность между экспоненциальной функцией и ее аппроксимациями, в — минимаксная ошибка в логарифмическом масштабе
б
а
в
Рис. 4. Аппроксимация функции ехр(х) на интервале [0, 1] усеченными рядами Тейлора степени п. а — экспоненциальная функция и ее аппроксимации, б — разность между экспоненциальной функцией и ее аппроксимациями, в — минимаксная ошибка в логарифмическом масштабе
в
Для сравнения те же данные для усеченных рядов Тейлора показаны на рис. 4. Следует отметить, что ряд Тейлора для функции exp(x) сходится очень быстро (в отличие от ряда Тейлора для функции log(x), например), однако даже в этом случае минимаксное приближение дает существенное преимущество.
Работа выполнена в Институте аналитического приборостроения Российской академии наук (Санкт-Петербург) в рамках темы FFZM-2022-0009 (номер гос. регистрации 122040600002-3) государственного задания Министерства науки и высшего образования Российской Федерации № 075-01157-23-00 от 29.12.2022.
При проведении вычислений использовалась лицен-зионно чистая программа Wolfram Mathematica версии 11 (Home Edition) [23], лицензированная одному из авторов.
СПИСОК ЛИТЕРАТУРЫ
1. Ланцош К. Практические методы прикладного анализа. Справочное руководство. М.: ГИФМЛ, 1961. 524 с.
2. Березин И.С., Жидков Н.П. Методы вычислений, Т. 1. М.: ГИФМЛ, 1962. 464 с.
3. Данилов Ю.А. Многочлены Чебышева. Минск: Вы-шэйшая школа, 1984. 157 с.
4. Пашковский С. Вычислительные применения многочленов и рядов Чебышева. М.: Наука, 1983. 384 с.
5. Грибкова В.П. Эффективные методы равномерных приближений, основанные на полиномах Чебышева. М.: Изд-во "Спутник", 2017. 193 с.
6. Гончаров В.Л. Теория интерполирования и приближения функций. М.-Л.: ГТТИ, 1934. 316 с.
7. Ахиезер Н.И. Лекции по теории аппроксимации. М.-Л.: ОГИЗ, 1947. 324 с.
8. Ахиезер Н.И. Лекции по теории аппроксимации. 2-е изд. М.: Наука, 1965. 407 с.
9. Berdnikov A., Solovyev K., Krasnova N., Golovitski A., Syasko M. Algorithm for Constructing the Chebyshev-Type Polynomials and the Chebyshev-Type Approximations with a Given Weight // Proceedings of 2022 Int. Conference on Electrical Engineering and Photonics (EExPolytech). P. 143-145.
DOI: 10.1109/EExPolytech56308.2022
10. Чебышев П.Л. Теория механизмов, известных под именем параллелограммов // П.Л. Чебышев. Избранные труды. М.: Изд-во Академии Наук СССР, 1955. С. 611-648.
11. Чебышев П.Л. Вопросы о наименьших приближениях, связанных с приближенным представлением функций // П.Л. Чебышев. Избранные труды. М.: Изд-во Академии Наук СССР, 1955. С. 462-578.
12. Чебышев П.Л. О функциях, наименее уклоняющихся от нуля // П.Л. Чебышев. Избранные труды. М.: Изд-во Академии Наук СССР, 1955. С. 579-608.
13. Бернштейн С.Н. О наилучшем приближении непрерывных функций посредством многочленов данной степени // Сообщения Харьковского математического общества. Вторая сер. 1912. Т. 13, вып. 2-3. С. 49-144.
14. Бернштейн С.Н. О наилучшем приближении функций нескольких переменных посредством многочленов или тригонометрических сумм // Труды МИАН СССР. 1951. Т. 38. С. 24-29.
15. De la Vallée-Poussin Ch.J. Sur les polynomes d'approximation et la représentation approchée d'un angle // Bulletin de la classe des sciences, Académie royale de Belgique, 1910. No. 12. P. 808-845.
16. De la Vallée-Poussin Ch.J. Leçons sur l'approximation des fonctions d'une variable réelle: professées à la Sorbonne. Paris: Gauthier-Villars, 1919.
17. Butzer P.L., Nessel R.J. Aspects of de La Vallée Poussin's work in approximation and its influence // Archive for History of Exact Sciences. 1993. Vol. 46. P. 67-95. DOI: 10.1007/BF00387727
18. Remez E. Sur un procédé convergent d'approximations successives pour determiner les polynomes d'appproxima-tion // Compt. Rend. Acad. Sci. Paris, 1934. Vol. 198. P. 2063-2065.
19. Remez E.Ya. Sur le calcul effectif des polynomes d'approximation de Tschebyscheff // Compt. Rend. Acad. Sci. Paris, 1934. Vol. 199. P. 337-340.
20. Дзядык В.К. Введение в теорию равномерного приближения функций полиномами. М.: Наука, Главная редакция физико-математической литературы, 1977. 508 с.
21. Ремез Е.Я. Основы численных методов чебышевского приближения. Киев: Наукова думка, 1969. 624 с.
22. Лоран П.-Ж. Аппроксимация и оптимизация. М.: Мир, 1975. 496 с.
23. Wolfram Mathematica: наиболее полная система для математических и технических вычислений. URL: https://www.wolfram.com/mathematica/ (обращение 21.03.2023).
Институт аналитического приборостроения РАН, Санкт-Петербург
Контакты: Бердников Александр Сергеевич, asberd@yandex. ru
Материал поступил в редакцию 28.03.2023
ISSN 0868-5886
NAUCHNOE PRIBOROSTROENIE, 2023, Vol. 33, No. 3, pp. 74-91
NUMERICAL ALGORITHM FOR MINIMAX POLYNOMIAL APPROXIMATION OF FUNCTIONS WITH A GIVEN WEIGHT
A. S. Berdnikov, S. V. Masyukevich
Institute for Analytical Instrumentation of RAS, Saint Petersburg, Russia
The article discusses a rapidly converging numerical algorithm for determining the coefficients of polynomials. The algorithm provides the optimal approximation of a given function in the minimax norm with a given weight on a given interval. The approximation is made under the condition that the weight function does not turn to zero on the considered interval, except, perhaps, for the initial and/or final points of the interval.
Keywords: minimax norm, Chebyshev polynomials, optimal approximation
INTRODUCTION
The use of function approximations that are optimal in terms of the minimax (uniform) norm has significant advantages over the simpler approximation of functions using the least squares method [1, 2]. As a practical tool in constructing such approximations, the Chebyshev polynomials of the first kind are used, deviating least from zero and having a weight function equal to one [3, 4]. The decomposition of a function into a truncated series consisting of polynomials that deviate least from zero (the Chebyshev polynomials in the case of a unit weight function) is one of the methods for the close construction of polynomial approximations that are optimal in the minimax norm [4, 5]. Also, to construct polynomial approximations that are close to the optimal minimax approximation, one uses an approximation over a discrete set of collocation points corresponding to the zeros of the Chebyshev polynomial with a degree one greater than the degree of the approximating polynomial [1, 2, 4, 5].
Such methods are characterized by convenience and ease of design, as well as having certain empirical prerequisites to obtain quite good accuracy. However, the simplified methods of approximation differ from the true approximation, which is optimal in terms of the minimax norm, and, therefore, the former provide an uncontrollably larger approximation error than the exact solution of the corresponding minimax optimization problem. The use of approximations that provide an exact solution to the minimax optimization problem for a given function (when possible) is definitely preferable.
The set of functions for which there is an explicit algebraic (analytical) representation for a polynomial approximation that is optimal in the minimax norm, is not too large (practically all such cases are given in [6-8]). In this paper, we study an improved version of a rapidly converging numerical algorithm for determining the coefficients of polynomials that deviate
least from a given function on a given interval with a given weight, which was previously considered in [9].
MINIMAX FUNCTION APPROXIMATION AND THE THEORY OF CHEBYSHEV
Formulation of the problem
Let there be a continuous function fx) and a finite interval [a, b] for the latter a continuous weight function o>(x) is given, which is strictly positive on this interval, except, perhaps, for the ends of the interval, for which it can turn into zero.1 It is necessary to find a polynomial p(x) of degree at most n with previously undefined coefficients a0, a1, a2, ..., an:
p(x) = a0 + a1x + a2x2 + ... + an-1xn-1 + anxn, (1)
which is the solution to the optimization problem
max |©(x) (fx) - p(x))| ^ min, (2)
where the variable x e [a, b] is maximized, and the coefficients a0, a1, a2, ..., an. are minimized.
Without loss of generality, we can assume that the function f(x) is not a polynomial of degree n or lower, otherwise the optimal solution for problem (2) is known in advance and coincides with the function fx).
1 For a weight function with zeros inside the interval [a, b], it will be required that the Chebyshev maxima and minima (and trial points for the numerical algorithm) do not coincide with the zeros of the weight function, and the signs of the alternating maxima and minima change in accordance with the sign of the weight function. Also, the interval [a, b] can be infinite on the right and/or left, but both the approximated function fx) and the weight function q(x) must tend to zero at infinity no slower than the power functions 1/xk, and with the same rate [6-8].
The Chebyshev criterion for the minimax polynomial approximation
For a polynomial p(x) of degree n to be a solution to problem (2), it is necessary and sufficient that there exist such a set of n + 2 points x0 < xj< x2 < ... < xn + ь belonging to the interval [a, b], and such a number s (positive or negative) that the following conditions are met:
—|s| < ®(x) (fix) -p(x)) < |s| при x e [a, b], (3) m(xt) (fixk) -p(xk)) = (-1)k s при k = 0, 1, .., n + 1. (4)
The points x0 < xb < x2 < ... < xn + ь with alternating positive minima and negative maxima equal in absolute value to the maximum modulus of the function under study are called the Chebyshev alternance (see [7, 8, 10-14]). Graphs of Chebyshev polynomials of the first kind [3], demonstrating the fulfillment of this condition on the interval [-1, +1] with the weight ®(x) = 1, are shown in Fig. 1.
Fig. 1. Chebyshev polynomials Tn(x), which deviate least from zero on the interval [-1, +1]. a — n = 5, 6 — n = 8, b — n = 16
The meaning of conditions (3), (4) is that for the polynomial p(x) there is a set of n + 2 test points xk at which the discrepancy œ(x) (fx) - p(x)) has local alternating negative minima and positive maxima equal to ±|e|, and the values of these minima and maxima are global on the interval [a, b]. This criterion is a special case of Chebyshev's theory of minimax approximation of functions using fractional rational functions [7, 8], but in the case of approximation using polynomials, the proof of the appropriate statements is simplified.
Statement 1 (The de la Vallée - Poussin theorem [7, 8, 15-17]). If the considered function œ(x) (fx) --p(x)) for some polynomial p(x) of degree n takes non-zero values X0, X1, ..., ÎN_ 1 with alternating signs, where N > n + 2, at N successive points x0 < x1 < ... < xN _ 1 in interval [a, b], then for any other polynomial r(x) of degree n,
max |œ(x) (fx) - r(x))| > min {|Xo|, |^|, ..., |V - i|}.
Proof. Let there be a polynomial r(x) satisfying the condition max |œ(x) (fx) - r(x))| < min (|X0|, |X1|, ..., |V - 1|). We assume that in the sequence xk the values ®(xk) (f(xk) - r(xk)) are strictly positive for even points, and strictly negative for odd points. (When the values are positive for odd points and negative for even points, the reasoning is similar.) This means that the condition œ(xk) (/(xk) -- r(xk)) < |Xt| = œ(xk) (/(xk) -- p(xk)) is satisfied at even points, and the condition œ(xk) (fx) - r(xk)) > -|Xk| = œ(xk) (fx) - r(xk)) is sa-
tisfied at odd points. As a result, the value œ(x) (p(x) -
- r(x)) is strictly negative at even points xk and strictly positive at odd points xk. The values of the weight function œ(xk) at the points xk,, for which the function œ(x) (fx) - r(x)) takes strictly positive or strictly negative values, do not vanish and, therefore, are strictly positive. Therefore, the polynomial p(x) - r(x) of degree n, which is not identically zero, alternately takes strictly positive and strictly negative values at least at n + 2 different points and therefore has at least n + 1 zeros of multiplicity 1. Hence, the polynomial n + 2 satisfying the condition max |œ(xk) (fxk) -
- r(xk))| < min (|Ao|, |Xi|, ., |V - i|), does not exist.
Comment. The de la Vallée - Poussin theorem allows one to obtain a lower estimate for solving the optimization problem (2). If X0, X1, ..., XN - 1 are local maxima and minima of the œ(x) (fx) - r(x)) with alternating signs, then this construction is called the de la Vallée - Poussin alternance, and for solving the optimization problem (2) it gives not only a lower bound, but also an upper bound equal to max {|X0|,
|H ..., |V-1|}.
Statement 2. If the conditions (3), (4) of the Chebyshev criterion are satisfied for the polynomial p(x), then the optimal value for the right side of the optimization problem (2) is equal to |e|, and the polynomial p(x) is its solution (possibly, one of many).
Proof. The Chebyshev alternance (3), (4) is a special case of the de la Vallée - Poussin alternance, therefore, according to Statement 1, the right side of the optimization problem (2) is not less than |e|. On the other hand, max|œ(x) (fx) - p(x))| = |e|. Therefore, the minimum of the optimization problem (2) is equal to |e|, and the polynomial p(x) is one of its solutions.
Statement 3. There is a polynomial p(x) that provides a solution to the optimization problem (2). (This excludes the case when there are polynomials pk(x) with constantly decreasing values for Pk = max |œ(x) (fix) -
- pk(x))|, while the minimum of this quantity is never reached on the set of polynomials of a given degree.)
Proof. The values of Pk = max |œ(x) (fx) - pk(x))| are bounded from below by zero; therefore, for the values Pk on the set of polynomials, there exists an infimum P. In accordance with the definition of the infimum, there exists a sequence of polynomials pk(x) of degree n:
Pk(x) = ao, k + a1, k x + a2, kx2 + ... + an - 1, kxn - 1+a„, kxn,
for which lim Pk = P as k ^ <x>. Therefore, starting from a certain number k, the values of the polynomials pk(x) on the interval [a + S, b - S], where S is small enough, will be bounded from above and below (we had to make an indent from the ends of the interval to take into account the case when the weight function
may be zero at the ends of the interval). From the Lagrange formula [1, 2] for a polynomial of degree n, which at n + 1 given points xk takes the values yk
f
P ( ^ )= I
j=l,n+l
y
j n
ï=1,n+1,ï ^ j
(x - X ) ( xj - X)
A
it follows that when at n + 1 fixed points xk the values of the polynomial yk are bounded above and below, then each of the coefficients of the polynomial is also bounded above and below. Therefore, in accordance with the Bolzano - Weierstrass theorem (also called the Bolzano - Weierstrass lemma) on the limit point, according to which from any infinite bounded sequence of points in the space Rn, one can select a convergent subsequence, in the sequence of polynomials pk(x) one can choose a subsequence of polynomials in which there is a limit for each of the coefficients of the polynomial. That is, there is a sequence of polynomials pk(x) for which lim Pk = P as k ^ w and the polynomial coefficients have limit values bj = lim a;-, k as k ^ w.
Consider a polynomial r(x) with coefficients bf r(x) = b0 + b1x + b2x2 + ... + bn _ x- 1 + bxn-
Since bj = lim aJk, then for any fixed value x the condition lim pk(x) = r(x) is satisfied as k ^ w, and the tendency to the limit is uniform on the considered finite interval of values x. From inequality
max |©(x) (fx) - r(x))| < max |o>(x) (fx) - pk(x))| + + max |q(x)| max |r(x) - pk(x)|
it follows that max |©(x) (fx) - r(x))| = P, because this value cannot be less than P, while as k ^ w the first term on the right side tends to P, and the second term tends to zero. This means that for the polynomial r(x) the value (2) reaches its lower bound.
Statement 4 (The Chebyshev theorem). The polynomial p(x), which provides a solution to the optimization problem (2), satisfies the Chebyshev criterion.
Proof. Consider the behavior of the function F(x) = = ®(x) (fx) - p(x)) on the segment [a, b]. To simplify the reasoning, we assume that the number of intervals on which the function F(x) changes sign is finite. This does not mean that the number of zeros of the function F(x) is also finite: let's assume the variant is when the function F(x) is equal to zero at any point of some interval inside the segment [a, b]. Let the points y1, y2, ., ym divide the segment [a, b] into m + 1 intervals, on each of which the function F(x) takes alternately positive and negative values, and let us assume the option when on the entire segment [a, b] the function F(x) takes only positive or only negative values. If there is continuous interval consisting entirely of zeros of the function F(x), this interval is attached to the
neighboring interval with positive or negative values of the function.
For intervals with positive values of F(x), we choose a point with the maximum value of the function on this interval (perhaps not the only one on this interval), and for intervals with negative values of F(x), we select the point with the minimum value of the function on this interval. We get a set of points y1, y2, ..., ym, — function F(x) zeros dividing the segment [a, b] into m + 1 intervals, for which the function F(x) retains its sign, changing it when crossing the boundary of the interval. Also we get a set of points x0, x1, x2, ..., xm, xm + 1 belonging these intervals and not coinciding with the ends of the intervals. This set consists of alternating positive maxima and negative minima Xo, Xb X^ ., Xm+1 of the function F(x) (Vallée - Poussin alternance) and the accompanying partition of the segment [a, b] into intervals on which the function F(x) retains its sign).
Let X = max |Xk|. If the maximum is equal to X on some interval, and the modulus of the minimum is less than X on the neighboring (next or previous) interval, then this interval, together with the subsequent interval with a positive maximum, regardless of its value, merges the current interval (see Fig. 2, a).
Fig. 2. Combining segments with alternating maxima and minima.
a — the case of a single minimum or several successively located minima that lie above the minimum level of the function F(x); 6 — the case of a single maximum or several successively located maxima that lie below the maximum level of the function F(x). The vertical arrows mark the minima and maxima, which correspond to the maximum deviation. Slanted arrows mark minima and maxima, which are less than the maximum deviation in modulus, so that the graph of the function can be shifted down or up, respectively, by a non-zero distance without going beyond the limits of the maximum deviation. Bold dots mark the zeros of the function F(x), which are the boundary of the combined interval with strictly positive or strictly negative maximum deviations of the function F(x)
The new interval has the property that such a positive constant can be subtracted from the function F(x) so that the positive maximum of the function on the interval decreases, but the absolute value of the negative minimum of the function on the interval does not increase so much as to exceed the new maximum of the function.
In the same way, if on some interval the minimum is equal to -X, and on the neighboring interval the maximum is less than X, then this interval, together with the subsequent interval containing a negative minimum, merges the current interval (see Fig. 2, 6). In this case, on the obtained interval, a positive con-
stant can be added to the function F(x) so that the absolute value of the negative minimum on the interval decreases, but the positive maximum of the function on the interval does not increase so much that it becomes greater than the absolute value of the new minimum.
Finally, on the interval [a, b] for a given polynomial p(x), a set of N intervals with points x0, xi, x2, ..., xN _ i on them, and the function F(x) alternating positive maxima and negative minima equal to ±A, where A = max |©(x) (fx) _ p(x))|. If N is greater than or equal to n + 2, then such a polynomial p(x) satisfies the Chebyshev criterion and, in accordance with Statement 1, will be a solution to the optimization problem (2).
Let N be less than or equal to n + 1. At the stage of combining the intervals, a set of points y1, y2, ..., yN _ i (interval boundaries), at which the function F(x) vanishes and changes its sign, was obtained. These points divide the interval [a, b] into smaller intervals, each of which has a point with a Chebyshev minimum or a Chebyshev maximum. There are also constants Sk > 0, which can be subtracted from the function F(x) on intervals with positive maxima or added to the function F(x) on intervals with negative minima, while reducing the minimax norm on the appropriate interval.
Without loss of generality, we can assume that the function F(x) takes positive values on the first interval; otherwise, you can change the sign of the function fx) and the polynomial p(x) simultaneously without changing the optimization problem (2).
Consider the function Q(x), which is the product of a polynomial of degree at most n and a weight function ©(x):
Q(x) = ©(x) (y1 _ x) (y2 _ x).(yN _ 1 _ x).
The function Q(x) is strictly positive on those intervals where there is a Chebyshev maximum of the function F(x), and strictly negative on those intervals where there is a Chebyshev minimum of the function F(x). If we subtract the function sQ(x) with a sufficiently small positive factor 5 (for example, we can choose 5 = (S / R), where R = max |Q(x)| and S = = min Sk > 0) from the function F(x), this will safely reduce the positive maxima of F(x) and decrease the negative minima of F(x), thereby reducing max |F(x)|. Therefore, such a polynomial p(x) cannot be a solution for the optimization problem (2), since it is possible to obtain an even smaller value for max |F(x)| when replacing the current polynomial p(x) of degree n with a new polynomial r(x) of degree n:
G(x) = F(x) _ (S / R) Q(x) =
= ©(x) (fx) _ p(x)) _ (S / R) ©(x) (y _ x). (yN _ 1 _ x) = = ©(x) (fx) _p(x) _ (S / R) (y _ x)...(yN_ 1 _ x)) = = ©(x) (fx) _ r(x)), max |G(x)| < max |F(x)|.
Conclusion: the optimal polynomial of degree n, considered in Statement 2, must satisfy the Chebyshev criterion.
Statement 5. The polynomial p(x) that satisfies the Chebyshev criterion and provides a solution to the optimization problem (2) is the only one.
Proof. The idea of the proof is borrowed from [7, 8]. Let there be two polynomials p(x) and r(x) of degree n, which are the solution to the optimization problem (2). For the degenerate case, for which max |©(x) (fx) _ p(x))| = 0 and max |©(x) (fx) _ r(x))| = = 0, the conditions fx) _ p(x) = 0 h fx) _ r(x) = 0 must be satisfied, which means that p(x) = r(x). Therefore, we will assume that max |©(x) (fx) _ p(x))| = = max |©(x) (fx) _ r(x))| > 0 so that the expressions ©(x) (fx) _ p(x)) and ©(x) (fx) _ r(x)) on the interval [a, b] have nonzero maxima and minima.
According to Statement 3, each of the polynomials p(x) and r(x) satisfies the Chebyshev criterion, and the value s ^ 0 for them is the same. Let x1, x2, ..., xn + 2 be Chebyshev deviation points with alternating maxima and minima equal to ±s, for the function F(x) = ©(x) (fx) _ p(x)). If for the polynomial r(x) the function G(x) = ©(x) (fx) _ r(x)) takes the same values ±s at the points x1, x2, ., xn + 2, except, maybe one point from this list, then the polynomials p(x) and r(x) are identically equal to each other. Indeed, in this case, the difference
d(x) = F(x) _ G(x) =
= ©(x) (fx) _ p(x)) _ ©(x) (fx) _ r(x)) = = ©(x) (r(x) _ p(x)), (5)
which is the product of a positive weight function ©(x) and a polynomial of degree at most n, vanishes at n + 2 points or, in extremis, at n + 1 points — and therefore the polynomial r(x) _ p(x) must be identically zero, since otherwise it cannot be explained that it has so many zeros (the weight function ©(x), even if it vanishes at the ends of the interval [a, b], is not equal to zero at the points of function F(x) minima and maxima).
Let among the Chebyshev deviation points x1, x2, ..., xn + 2 for the function F(x) there be at least two points xk and xk + m, at which the value of the function G(x) does not coincide with the value of the function F(x) in the point (minimum or maximum of the function F(x)). It is required to prove that in this case inside the interval [xk, xk+m] there is at least one more point that does not coincide with the points x1, x2, ., xn + 2, and at which the values of the functions F(x) and G(x) coincide (although, perhaps, at this point the common value of the functions F(x) and G(x) is not equal to ±s), and, consequently, the function d(x) given by formula (5) has an additional unaccounted zero.
If at some point xj the condition d(xj) ^ 0 is satis-
fied, then the sign of the difference d(xf) given by formula (5) coincides with the sign of F(xf). Indeed, if xj is the point of the negative minimum of the function F(x), then G(xj) > F(x;), and if xj is the point of the positive maximum of the function F(x), then G(xf) <
< F(x). If
d(xk) * 0, d(xk + 1) = 0, d(xk + 2) = 0, ., d(xk + m-1) = 0, d(xk + m) * 0, (6)
then the expressions (-1)k F(xk) and (-1)k + m F(xk+m) have the same signs and, therefore, the expressions d(xk) and (-1)m d(xk + m) have the same signs. Therefore, the parity of the number of sign changes of the function d(x) on the interval [xk, xk+m] coincides with the parity of the number m, and thus the parity of the number of roots (taking into account their multiplicity) for the function d(x) on the interval [xk, xk + m] coincides with the parity of the number m. According to (6), the function d(x) on the interval [xk, xk+m] has at least m - 1 unequal roots and, therefore, either one more root must exist, or one of the roots xk + b xk + 2, ., xk + m _ 1 has an even multiplicity of at least 2. After performing this operation for all intervals bounded by points xk with non-zero values d(xk), we obtain that the number of roots of the function d(x) on the interval [a, b], taking into account the multiplicity of the roots, will be no less than n + 1. Therefore, the polynomial r(x) - p(x) of degree at most n must be identically zero in the formula d(x) = ©(x) (r(x) - p(x)).
ALGORITHM FOR CONSTRUCTING A MINIMAX POLYNOMIAL
APPROXIMATION OF A FUNCTION
An algorithm for the numerical construction of approximating polynomials is derived from the Cheby-shev criterion for the minimax approximation of a function. This algorithm is a modified and optimized E.Ya. Remez algorithm [18-22]:
1. An initial set of points x0 < xi < x2 < ...< xn + 1, on the interval [a, b] is chosen arbitrarily. If the weight function ©(x) vanishes at the point x = a, it is necessary to satisfy the condition a < x0. If the weight function ©(x) vanishes at the point x = b, then the condition x0 < b must be met.
2. For the current set of points a < x0 < x1 <
< x2 < ...< xn + 1 < b, an auxiliary function g(x) = = a-f(x) - p(x) is constructed, where a is some constant, and p(x) is a polynomial of degree n with some coefficients. The free coefficients of the function g(x), i.e. multiplier a and coefficients of polynomial p(x) are chosen so that the conditions ©(xk)g(xk) = (-1) are satisfied for k = 0, 1, ..., n + 1. Such a function exists and is uniquely defined (see further).
3. Now we need to determine the algebraic equation g(x) = 0 roots, which are on the interval [a, b].
Since the values of the function g(x) have different signs at the edges of the intervals [xk, xk + 1], then inside each such interval there is at least one root of the equation g(x) = 0. If the function f(x) changes slowly enough, then on each of the intervals [xk, xk + 1] there is exactly one root of multiplicity 1, and the numerical determination of these roots is not difficult. If this condition is violated and there are several roots on the interval under consideration, all these roots are added to the list (multiple roots with even multiplicity, for which the values of the function g(x) on both sides of the root have the same sign, are not included in the list). Upon completion of the search for roots, the beginning and the end of the interval y0 = a and ym + 1 = b are added to the list of values a < y1 < y2 < .< ym < b, which are roots of odd multiplicity of the function g(x) on the interval [a, b].
4. Now the interval [a, b] is a set of at least n + 2 intervals [yk, yk + 1], which are in contact with their end points, and on each interval the function g(x) retains its sign, and the sign of the function changes when the boundary between the intervals is crossed. Consider the function G(x) = ©(x) g(x), which on the interval [a, b] takes values of the same sign as the function g(x). Let us find the global (for the considered interval) ex-tremum G(x) and the point zk e [yk, yk + 1], where this extremum is located, for the function G(x) on each interval x e [yk, yk+1]. For intervals with positive values of the function G(x), one should look for a positive maximum, and for intervals [yk, yk + 1] with negative values of the function — a negative minimum.
5. It is necessary to exclude extra points if the number of determined points Gk = G(zk), zk e [yk, yk + 1] with alternating signs is greater than n + 2. To do this, the least modulo value of Gk is chosen. If the interval with this value is the first or last in the list, the value of Gk and the corresponding point zk are discarded. If the interval with this value is in the middle of the list, check the previous and subsequent intervals and discard not only the considered value of Gk, but also that of the two neighboring values Gk_ 1 and Gk + 1, which will be less in absolute value. Indeed, the signs of the values Gk _ 1 and Gk + 1 are opposite in sign to the value of Gk, and since due to the removal of the value Gk from the list, two values of the same sign will be nearby, then one of the two values Gk + 1 and Gk _ 1 must be discarded. It's important: there is a risk of getting n + 1 points from n + 3 zk points, instead of the required n + 2 points at once since one or two points are removed at each step, — in this case, you should keep the point Gk in the list, and delete instead one of the start or end points, which is less in absolute value. The process continues until the list contains n + 2 extrema Gk = G(zk), corresponding to n + 2 points zk with alternating positive maxima and negative minima of the function G(x) = ©(x) g(x).
6. For the existing list of points a < z0 < z1 < < z2 < ...< zn + 1 < b, let's take a closer look at the values of extrema Gk = ©(zk) g(zk) = (-1)kek. If ek « 1 within the given accuracy, then a polynomial that satisfies the requirements of the Chebyshev criterion is determined. Moreover, it follows from the de la Vallée -Poussin theorem on the best approximation of a function by polynomials (Statement 1) that the exact optimum for problem (2) lies in the range between min |ek / a| | and max |ek / a|, (where a is the factor in the expression g(x) = a fx) - p(x)). If the condition ek « 1 is not satisfied, we replace the trial points a < x0 < < x1 < x2 < .< xn + 1 < b with points a < z0 < z1 < z2 < .< zn + 1 < b and return to step 2.
7. It remains to normalize the function g(x) = = a fx) - p(x) so that the coefficient of the function fx) becomes equal to 1. Such a normalization is always feasible, since the coefficient a cannot be zero, otherwise it will turn out that there is a non-identical zero polynomial p(x) of degree n, which has n +1 roots located between neighboring maxima and minima of different signs in the amount of n + 2 pieces.
Obviously, if this algorithm converges, then it converges to a polynomial that satisfies the Cheby-shev criterion, and thus the result of the algorithm is a solution to the optimization problem (2), and the only such solution. Theoretically, the algorithm always converges at the rate of a geometric progression: for any continuous weight function and for any degree of n, there are numbers C > 0 and 0 < \ < 1 for which
max |©(x) (ai-f(x) -pi(x))| < 1 + Ck\
where i is the iteration number. The proof is similar to the proof of the convergence of the E.Ya. Remez algorithm in the monograph [20]. Unfortunately, due to rounding errors, the algorithm may diverge for "bad" functions fx) and ©(x), namely, for functions that do not change slowly on the interval [a, b], but instead show multiple oscillations, so that they have more than one root for the intervals [xk, xk + 1] or more than one local maximum or minimum for the intervals [yk, yk + 1]. In such cases, the algorithm can perform an infinite loop, in fact switching between the same list of tested polynomials. This behavior of the algorithm should be controlled manually by setting the maximum number of cycles that are allowed to run. Another useful criterion for terminating the operation of an algorithm that has fallen into an infinite loop can be, at the next step, to check the increase of the spread between the maxima and minima ±ek determined on the intervals [yk, yk + 1] compared to this value of the previous step. Also, a suspicious attribute, signaling the looping of the algorithm, may be the regular occurrence of a situation, at Step 5, when "cleaning" points Gk, at some step there are less than n + 2 points in the list and it becomes necessary to take a step
back, removing the first or the last point, because even if the optimal polynomial p(x) is the identical zero (for example, such a situation arises when trying to approximate with a polynomial of finite degree of the function sin(1 / x) on the interval [0, 1]), then, according to Statement 4, there must still be a set of n + 2 points zk that form the Chebyshev alternance.
Since the weight function œ(x) can vanish at an end of the interval [a, b], then to perform step 2 it will be required that the test points chosen at step 1 lie strictly inside the interval [a, b]: a < x0 < x1 < < x2 < ...< xn+1 < b. It can be seen that if œ(a) = 0 or œ(b) = 0 and, a <x0 < xi < x2 < .< x„+i < b, then the condition a < x0 < x1 < x2 < .< xn+1 < < b will be preserved for all subsequent iterations, because the points zk, at step 6, will also have to satisfy the condition a < z0 < z1 < z2 < ...< zn + 1 < b.
To calculate the coefficients at step 2, it is impossible to use an explicit solution of the system of linear equations to determine the coefficients of the function g(x) = afx) - p(x), since the conditionality of the matrix of linear equations (except for the first column, coinciding with the Vandermonde matrix) grows rapidly with the degree of the polynomial: as a result, small errors in rounding floating point numbers lead to completely unsatisfactory values for unknown coefficients. One of the efficient ways to determine the desired g function g(x) = a-fx) - p(x) at step 2 is as follows. From the conditions p(xk) = a-fxk) - (-1)k / œ(xk), which are equivalent to the conditions œ(xk) (a-fxk) -- p(xk)) = (-1)k, using the Lagrange formula [1, 2], we determine the coefficients of the polynomial p(x), which at given points xk takes the given values. These coefficients will be linear expressions with a free coefficient a. After that, for the equation œ(xk) (a-fxtk) -p(xk)) = (-1)k, which is a linear equation with respect to the variable a, the unknown coefficient a is determined.
Another, more efficient way to determine the coefficient a and the coefficients of the polynomial p(x), for which the conditions will be satisfied at step 2
œ(xk) (a-fxk) -p(xk)) = (-1)k при k = 0, 1, ., n+1,
consists of the following. To fulfill the equalities p(xk) = a-fx^ - (-1)k / œ(xk) according to Newton finite difference method [1, 2], we determine such a polynomial r(x) for which the conditions r(xk) = = (-1)k / a>(xk) are satisfied if k = 0, 1, ., n, and a polynomial s(x) for which the conditions s(xk) = f(xk) are satisfied if k = 0, 1, ..., n. Then from the condition
2 Instead of the explicit Lagrange formula, especially for large degrees of polynomials, it is convenient to use the Aitken iterative method [2]. Although this algorithm requires more memory, it provides significantly better coefficient accuracy and requires less computing time.
©(xn + 1) (a-f(xn + 1) - as(xn + 1) + r(xn + 1)) = (-1) n + 1
the constant a is defined. After that, the desired polynomial p(x) is determined in accordance with the formula p(x) = as(x) - r(x). This method can also be used in the following form: we calculate a polynomial s(x) of degree n + 1 for which the conditions s(xk) = fxk) are met at the points xk if k = 0, 1, ..., n + 1, and a polynomial r(x) of degree n + 1 for which the conditions r(xk) = (-1) / ©(xk) are met at the points xk if k = 0, 1, ..., n + 1. It remains to choose the coefficient a so that the polynomial as(x) - r(x) has degree n, which is always possible, unless the function fxk) is a polynomial of degree at most n.
Note that at intermediate stages, there is no need to determine the coefficients of the polynomial p(x) explicitly, it is enough to be able to calculate its value at an arbitrary point. Therefore, in the case of computer calculation of this algorithm, until the final result is obtained, one can use the calculation of the values of the polynomial p(x) according to the formula p(x) = = as(x) - r(x), where for the auxiliary polynomials s(x) and r(x) the Newton method of finite differences is used.
To improve the operation of the algorithm at the initial stage, at step 1, it is recommended to choose the zeros of the Chebyshev polynomial of the first kind of degree n + 2, shifted and scaled from the interval [-1, +1] to the interval [a, b] for the initial points a < x0 < < x1 < x2 < ...< xn + 1 < b. At step 2, instead of the condition ©(xk) p(xk) = (-1)k, it is recommended to use the condition ©(xk) p(xk) = (-1)k(4 / (b - a))-n - 1 to avoid unreasonably large coefficients of the polynomial p(x) and the function f(x) (this choice of scale corresponds to the oscillations of the Chebyshev polynomial of degree n + 1, recalculated from the interval [-1, +1] to the interval [a, b] and scaled so that the leading coefficient is equal to 1).
Example. Let fx) = exp(x), a = 0, b = 1, ©(x) = 1. For n = 1, 2, ., 8 as a result of the described algorithm, we obtain polynomials that provide the optimal (in terms of minimax norm) approximation of the function exp(x) on the interval [0, 1]. The approximation error depending on the polynomial degree is indicated in Tab. The dependence of the approximation error on the degree of the polynomial is shown in graphical form in Fig. 3.
Tab. Error of the function exp(x) approximation on the interval [0, 1] with polynomials of degree n, optimal in minimax norm with weight q(x) = 1
Fig. 3. Approximation of the function exp(x) on the interval [0, 1] with minimax polynomials of degree n. a — exponential function and its approximations, 6 — difference between the exponential function and its approximations, b — minimax error on a logarithmic scale
For comparison, the same data for the truncated Taylor series are shown in Fig. 4. It should be noted that the Taylor series for the function exp(x) converges very quickly (in contrast to the Taylor series for the function log(x), for example), but even in this case, the minimax approximation gives a significant advantage.
Fig. 4. Approximation of the function exp(x) on the interval [0, 1] with a truncated Taylor series of degree n.
a — exponential function and its approximations, 6 — difference between the exponential function and its approximations, b — minimax error on a logarithmic scale
REFERENCES
1. Lanczos C. Applied Analysis. New York, Prentice Hall Inc., 1956. (Russ. ed.: Lantsosh K. Prakticheskie metody prikladnogo analiza. Spravochnoe rukovodstvo. Moscow: GIFML, 1961. 524 p.).
2. Berezin I.S., Zhidkov N.P. Metody vychislenii, T. 1 [Methods of computing. Vol. 1]. Moscow, GIFML, 1962. 464 p. (In Russ.).
3. Danilov Yu.A. Mnogochleny Chebysheva [Chebyshev polynomials]. Minsk, Vyshehishaya shkola, 1984. 157 p. (In Russ.).
4. Paszkowski S. Numerical Applications of Chebyshev Polynomials. Warszawa: Panstwowe Wydawnictwo Nau-kowe, 1975. (Russ. ed.: Pashkovskii S. Vychislitel'nye primeneniya mnogochlenov i ryadov Chebysheva. Moscow, Nauka, 1983. 384 p.).
5. Gribkova V.P. Ehffektivnye metody ravnomernykh prib-lizhenii, osnovannye na polinomakh Chebysheva [Effective methods of uniform approximation based on the Po-lynomas of Chebyshev]. Moscow, Sputnik Publ., 2017. 193 p. (In Russ.).
6. Goncharov V.L. Teoriya interpolirovaniya i priblizheniya funktsii [Theory of interpolation and approximation of functions]. Moscow, Leningrad, GTTI Publ., 1934. 316 p. (In Russ.).
7. Achiezer N.I. Lektsiipo teorii approksimatsii [Lectures on approximation theory]. Moscow, Leningrad, OGIZ Publ., 1947. 324 p. (In Russ.).
8. Achiezer N.I. Lektsii po teorii approksimatsii. 2-e izd [Lectures on approximation theory, 2nd ed.]. Moscow, Nauka Publ., 1965. 407 p. (In Russ.).
9. Berdnikov A., Solovyev K., Krasnova N., Golovitski A., Syasko M. Algorithm for Constructing the Chebyshev-Type Polynomials and the Chebyshev-Type Approximations with a Given Weight. Proceedings of2022 Int. Conference on Electrical Engineering and Photonics (EExPo-lytech), 2022 Oct 21-22, pp. 143-145. DOI: 10.1109/EExPolytech56308.2022
10. Chebyshev P.L. [Theory of mechanisms known by the name of parallelograms]. P.L. Chebyshev. Izbrannye trudy [See P.L. Chebyshev. Selected works], Moscow, Izd. Akademii Nauk SSSR, 1955, pp. 611-648. (In Russ.).
11. Chebyshev P.L. [Questions about the best approximations related to approximate representations of functions]. P.L. Chebyshev. Izbrannye trudy [See P.L. Chebyshev. Selected works], Moscow, Izd. Akademii Nauk SSSR, 1955, pp. 462-578. (In Russ.).
12. Chebyshev P.L. [On the functions that are the least deviating from zero]. P.L. Chebyshev. Izbrannye trudy [See P.L. Chebyshev. Selected works], Moscow. Izd. Akademii Nauk SSSR, 1955, pp. 579-608. (In Russ.).
13. Bernstein S.N. [On the best approximation of continuous functions through polynomials of a given degree]. Soobshcheniya Khar'kovskogo matematicheskogo obsh-chestva. Vtoraya ser. [Transactions of the Kharkov mathematical society. 2 series], 1912, vol. 13, no. 2-3, pp. 49-144. (In Russ.).
14. Bernstein S.N. [On the best approximation of the functions of several variables through polynomials or trigonometric sums]. Trudy MIAN SSSR [Proceedings of the Steklov Institute of Mathematics], 1951, vol. 38, pp. 2429. (In Russ.).
15. De la Vallée-Poussin Ch.J. Sur les polynomes d'approximation et la représentation approchée d'un angle. Bulletin de la classe des sciences, Académie royale de Belgique, 1910, no. 12, pp. 808-845.
16. De la Vallée-Poussin Ch.J. Leçons sur l'approximation des fonctions d'une variable réelle: professées à la Sorbonne. Paris, Gauthier-Villars, 1919.
17. Butzer P.L., Nessel R.J. Aspects of de La Vallée Poussin's work in approximation and its influence. Archive for History of Exact Sciences, 1993, vol. 46, pp. 67-95. DOI: 10.1007/BF00387727
18. Remez E. Sur un procédé convergent d'approximations successives pour determiner les polynomes d'appproxima-tion. Compt. Rend. Acad. Sci., Paris, 1934, vol. 198, pp. 2063-2065.
19. Remez E.Ya. Sur le calcul effectif des polynomes d'approximation de Tschebyscheff. Compt. Rend. Acad. Sci., Paris, 1934, vol. 199, pp. 337-340.
20. Dzyadyk V.K. Vvedenie v teoriyu ravnomernogo pribliz-heniya funktsii polinomami [Introduction to the theory of uniform approximation of functions by polynomas]. Moscow, Nauka Publ., Glavnaya redaktsiya fiziko-matematicheskoi literatury, 1977. 508 p. (In Russ.).
21. Remez E.Ya. Osnovy chislennykh metodov chebyshevsko-go priblizheniya [Fundamentals of the numerical methods of the Chebyshev approximation]. Kyiv, Naukova Dumka, 1969. 624 p. (In Russ.).
22. Laurent P.-G. Approximation et optimization. Paris, Herrmann, 1972 (Russ. ed.: Loran P.-Zh. Approksimatsiya i op-timizatsiya. Moscow: Mir Publ., 1975. 496 p.).
23. Wolfram Mathematica: naibolee polnaya sistema dlya matematicheskikh i tekhnicheskikh vychislenii [The world's definitive system for modern technical computing]. URL: https://www.wolfram.com/mathematica/ (accessed 21.03.2023). (In Russ.).
Contacts: Berdnikov Aleksandr Sergeevich, Article received ЬУ the editorial office on 28 03 2023
