Сравнение экстраполяционных методовобработки результатов измерений Текст научной статьи по специальности «Математика»

УДК 621.317: 519.67

Л. Н. Бондаренко, Д. И. Нефедьев

СРАВНЕНИЕ ЭКСТРАПОЛЯЦИОННЫХ МЕТОДОВ ОБРАБОТКИ РЕЗУЛЬТАТОВ ИЗМЕРЕНИЙ

L. N. Bondarenko, D. I. Nefed'ev

COMPARISON OF EXTRAPOLATION METHODS PROCESSING OF RESULTS OF MEASUREMENTS

Аннотация. Рассмотрены некоторые методы интерполяции и экстраполяции, применяемые при обработке результатов измерений. С помощью тестовых задач проведено сравнение методов на базе многочленов, сплайнов и рациональных функций. Результаты сравнения показывают, что в ряде случаев для экстраполяции наиболее выгодно использовать алгоритм Кронекера-Чебышева.

Abstract. Some methods of interpolation and the extrapolations applied at processing of results of measurements are considered. By means of test tasks comparison of methods on the basis of polynomials, splines and rational functions is carried out. Results of comparison show that in some cases for extrapolation it is most favorable to use of Kronecker-Chebyshev algorithm.

Ключевые слова: интерполяция, экстраполяция, многочлен, сплайн, рациональная функция, алгоритм Кронекера-Чебышева.

Key words: interpolation, extrapolation, polynomial, spline, rational function, Kronecker-Chebyshev algorithm.

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

Если в фиксированных точках (узлах) дискретизации {tt}”, где t1 < t2 <... < tn и все

tt е [a, b], найдены значения {f (tt )}П измеряемой величины f (t), то для решения задачи интерполяции часто используются многочлены, значения которых совпадают со значениями f (t) в узлах интерполяции [1], или полиномиальные сплайны, построенные на сетке {t;-}П [2].

Для решения задачи экстраполяции также возможно применение многочленов и сплайнов, используемых при интерполяции, но особенности функции f (t), расположенные в комплексной плоскости вблизи промежутка [a, b], могут сильно ухудшить это решение.

Нередко значения, принимаемые функцией f (t) вне отрезка [a, b], представляют наибольший интерес. Это характерно для изучения процессов, имеющих фазовые переходы,

.........4 Измерение. Мониторинг. Управление. Контроль

я

резонансные значения и т.п. В этих случаях важные для исследователя значения измеряемой величины могут находиться вне рассматриваемого промежутка [a, b].

■

; Следует отметить, что даже при решении задачи интерполяции многочлены и сплайны

могут не дать приемлемого результата, причем при построении интерполяционных многочле-\ нов существенную роль играет правильный выбор узлов интерполяции {tt}” .

Так, например, часто используемые на практике равноотстоящие узлы интерполяции являются наихудшими для рассматриваемой задачи, т.е. уже при достаточно невысокой степени интерполяционного многочлена, равной n -1, он может в ряде точек сильно отличаться от интерполируемой функции [1]. Положение значительно улучшается при выборе асимптотически оптимальных чебышевских узлов интерполяции. Если отрезок [a, b] преобразовать в

■

; стандартный отрезок [-1,1] с помощью линейной замены

■

я

I b - a . л. 2t - a - b r ,

t = a + —— (t +1), t = —----, те [-1,1], (1)

■ 2 b - a

я

■

l то чебышевские узлы интерполяции описываются следующими соотношениями:

■

(2n - 2i + 1)я b - a

T = cos----- ------, ti = a + —— (Ti +1), i = 1,..., n. (2)

; 2n 2

■

При преобразовании (1) функция f (t) заменяется функцией g (t) , а интерполяционный многочлен Hn-х(т) для функции g (t) = f (t) на узлах (2) определяется выражением

я

■

Hn-1(T) = I; * hjTj (т), (3)

j=о

■

■

в котором звездочка у знака суммы означает, что при j = 0 соответствующее слагаемое берет-\ ся с коэффициентом 1/2, Tj (т) = cos j arccos т, j = 0,1,..., n -1, - многочлены Чебышева, орто-

; гональные на отрезке [-1,1] с весом (1 -т2)-1/2, а коэффициенты {И,}П-1 в выражении (3) лег-

ко вычисляются по формуле [1]

■

hj = -|g(Ti) Tj (Ti), j = 0,1,., n -1.

! n i=1

■

Если особенностями функции f (t) являются полюсы, то рациональные функции в ка-■ u

честве аппроксимирующего агрегата выступают как главный конкурент многочленов и поли-

I номиальных сплайнов. Так, например, для практически важного случая, когда f (t) представ-

■

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

Так как рациональная функция инвариантна относительно преобразования (1), рассмотрим для простоты рациональную аппроксимацию функции g(т) на отрезке [-1,1], получаемой из f (t) с помощью этого преобразования. В этом случае функция g (т) аппроксимируется дробью вида Rk m (т) = Pk (т) / Qm (т), где многочлены Pk (т) и Qm (т) имеют соответствен-

■

! но степени: deg Pk (т) < к и deg Qm (т) < m . По сравнению с полиномиальным приближением

■

рациональная аппроксимация в ряде случаев более эффективна, причем она часто позволяет получить информацию об особенностях аналитической функции.

Если значения функции получены путем измерений на фиксированных узлах, то естественным аппроксимирующим агрегатом является рациональная функция, значения которой в узлах совпадают со значениями данной функции. Такая рациональная функция называется многоточечной аппроксимацией Паде, или рациональной интерполяцией, а соответствующая задача интерполяции рациональными функциями называется задачей Коши-Якоби [3].

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

вычислений [4] и реализованный, например, в виде процедуры в пакете Мар1е, использует только точную арифметику, и поэтому исходные данные для этой процедуры должны задаваться только в алгебраическом виде (их нельзя записывать в системе с плавающей точкой).

Рассмотрим алгоритм рациональной интерполяции Кронекера [3], который из-за ряда существенных недостатков был модифицирован для использования в задаче структурнопараметрической идентификации в частотной области и назван алгоритмом Кронекера-Чебышева, причем некоторые недостатки в результате такого преобразования превратились в определенные достоинства [5-7].

В задаче рациональной интерполяции по заданному набору точек {(тг-, gl )}П , gi = g(тг-)

требуется построить рациональную функцию Як т (т) = Рк (т)/ Qm (т), которая в узлах {тг-}П, тг- е [-1,1], совпадает с функцией g(т), т.е. в этих узлах должны выполняться равенства Як т (тг-) = gi, причем число неизвестных параметров этой задачи т + к + 1 = п .

Алгоритм Кронекера получения диагональной аппроксимации Як т (т), для которой к = [(п -1) / 2], т = [п / 2], где [•] - целая часть числа, заключается в построении интерполяционного полинома Нп-1(т), удовлетворяющего равенствам Нп-1(тi) = gi, i = 1,..., п, и применении рекуррентной процедуры, состоящей в вычислении последовательностей многочленов:

- {Р п - у (т)}т=:() , где Рп-1(т) = Нп-1(т) и Рп (т) = Г1 (т-тг );

i=1

- {Q^■ (т)}^ где Q-l(т) = 0 и Qо(т) = ь

а также нахождении параметров {а у, р у }т-1 с помощью следующих простых соотношений:

Рп-у -2 (т) = (ау т + Ру ) Рп-] -1(т) - Рп-] (т) ; (4)

Q у+1 (т) = (а у т + Р у) Q у (т) - Q у-1 (т), (5)

которые применяются при значениях итерационного шага у = 0,1, ..., т -1.

Если степени всех многочленов {Рп_у(т)}т=+01 совпадают с их индексами, то эта рекуррентная процедура не вырождена и на ее -м шаге из соотношения (4) находятся параметры

ау-, Ру с использованием частного от деления многочлена Рп ■ (т) на многочлен Рп-у-1 (т), а затем по формулам (4) и (5) вычисляются многочлены Рп-у-2(т) и Qу+l(т) .

На последнем (т -1) -м шаге алгоритма Кронекера, кроме параметров ат-1, Рт-1, определяются многочлены Рк (т) , Qm (т) и записывается искомая функция Як,т (т) = Рк (т)/ Qm (т) . Соотношения (4) и (5) также позволяют легко получить равенство детерминантов

1? 7 Г Г Рп—у—1(т) Рп—у—2(т)

б—1(т) Я](т) б](т) Я у+1(Т)

из которого следует, что при невырожденности алгоритма, т.е. при выполнении условий degР -(т) = п — у, у = 0,1,..., т — 1, детерминанты не зависят от итерационного шага у, выполняются соотношения

Рп—у — 2(т) Рп—у — 1(т) Рп СО . 0 1 1

---------=---------------------------, ] = 0,1, ..., т — 1,

0у.+і(Т) бу.(Х) бу(т) Єу+1(Т) 7 ’

которые показывают, что все рациональные функции Рп—у—2 (т) / Q у+1 (т) , вычисляемые по алгоритму Кронекера на шагах у = 0,1, ..., т — 1, совпадают со значениями исходной функции g (т) в узлах интерполяции (т і Щ, ті є [—1,1].

Основными недостатками алгоритма Кронекера являются возможность его вырождения и высокая чувствительность к изменению исходных данных.

Первый недостаток можно легко устранить, если в соотношениях (4) и (5) вместо линейной функции (ttyT + Рj) использовать многочлен более высокой степени с неопределенными коэффициентами. Тогда вычисление частного от деления многочлена P -(t) на многочлен Pn_j_j(t) позволяет найти все эти коэффициенты.

Для устранения второго недостатка было предложено для записи многочленов {Pn_j(t)}™+o , j = 0,1,..., m _ 1, использовать чебышевский базис {Tj(т)}П, а в качестве узлов интерполяции рассматривать чебышевские узлы (2), являющиеся нулями многочлена Чебышева Tn(t) = cosnarccos t . Для записи многочленов {Qj (т)}т1, j = 0,1, ..., m _ 1, в модифицированном алгоритме Кронекера применяется степенной базис {т] }0.

Модифицированная процедура Кронекера была названа алгоритмом Кронекера-Чебышева. В этом алгоритме используются многочлены Pn (т) = Tn (т) и Pn_1 (т) = Hn_1 (т),

где Hn_j(т) находится по формуле (3), причем вычисления, связанные с соотношением (4), производятся в чебышевском базисе, а связанные с соотношением (5) - в степенном базисе.

Чувствительность алгоритма Кронекера-Чебышева к изменению исходных данных существенно ниже чувствительности алгоритма Кронекера, что достигнуто за счет построения исходного многочлена Hn_x( т) на асимптотически оптимальных чебышевских узлах и проведения вычислений, связанных с соотношением (4), в чебышевском базисе.

Так как многочлены Чебышева можно определить для комплексной переменной p, для

чего достаточно заменить вещественную переменную т на p = ja, j = >/—1, юе[_1,1], алгоритм Кронекера-Чебышева подходит также для частотной области, где он был успешно использован для решения соответствующих задач идентификации, в которых результат должен быть представлен в виде рациональной функции.

Также для частотной области было разработано обобщение алгоритма Кронекера - алгоритм Кронекера-Чебышева-Ахиезера, который позволяет рассматривать задачу идентификации в любом диапазоне частот [6, 7].

При проведении компьютерных вычислений вырождение отмеченных алгоритмов практически не встречается, но возможность вырождения можно использовать для остановки этих алгоритмов на некотором промежуточном шаге j < т _ 1, в результате чего знаменатель результирующей рациональной функции будет иметь степень s < т . Отметим, что число s является важным структурным параметром при идентификации в частотной области.

Для такой остановки алгоритмов Кронекера-Чебышева и Кронекера-Чебышева-Ахиезера в [6] был введен параметр tol, характеризующий задаваемый допуск.

Так, в случае невырожденности алгоритма Кронекера-Чебышева при задании tol = 0 в результате интерполяции функции g (т) на отрезке [_1,1] будет получена рациональная

функция Rk т (т) = Pk (т) / Qm (т) . При увеличении задаваемого параметра tol можно получить

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

Следует отметить, что использование только диагональных рациональных аппроксимаций Rk т (т) = Pk (т) / Qni (т), при которых степень знаменателя превышает степень числителя

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

Алгоритмы Кронекера-Чебышева и Кронекера-Чебышева-Ахиезера реализованы в виде соответствующих процедур в пакете аналитических вычислений Maple. При этом была учтена возможность их применения как в вещественной, так и комплексной областях с использованием задаваемого параметра tol, причем при выборе tol = 0 были исключены основные случаи вырождения этих алгоритмов.

В результате тестирования алгоритма Кронекера-Чебышева в вещественной области на ряде показательных, тригонометрических, гиперболических и бесселевых функций по-

f (t) = 1

1 + 2512 ’

для которой полиномиальная интерполяция, рассматриваемая на равноотстоящих узлах в промежутке [—1,1], расходится в промежутке у - j t j < 1, где значение у ~ 0,72б . Это происходит из-за близости полюсов функции Рунге, равных ±0, 2 ] , где ] =V—I , к отрезку [—1,1].

Для функции Рунге при t є [0,1] и n = б вычислялись три интерполяционных агрегата:

- интерполяционный многочлен P5(t), построенный на чебышевских узлах (2);

- кубический сплайн S3(t) дефекта 1, построенный на равномерной сетке;

- рациональная функция, построенная с помощью алгоритма Кронекера-Чебышева.

C помощью пакета Maple был найден многочлен P (t), для чего использовалась формула (3), S3(t) находился с помощью встроенной процедуры «Spline», а рациональная функция вычислялась с помощью разработанной процедуры при выбранном значении параметра tol = 0,1, причем все вычисления проводились с десятью значащими цифрами, а в дальнейших записях все нули, не влияющие на результат, опущены.

Значения абсолютных погрешностей на отрезке [0,1] оцениваются неравенствами

тах | f (0 - Р5(0| < 0,374 ; тах | f (0 - 53(01 < 0,0787 ; тах | f (?) - Я] 2(t)| < 0,739-10 ге[0,1] ге[0,1] /е[0,1]

а полученная рациональная функция имеет вид

1,000000176 - 0,22475924 • 10-6 t

—7

R12(t)= — 6 2 .

1,00000025 _ 0,35 -10_61 + 2512

Графики функции Рунге и ее полученных аппроксимаций на отрезке [0,1] показаны на рис. 1, на котором график функции Рунге сливается с графиком рациональной аппроксимации R\2(t), полученной с помощью алгоритма Кронекера-Чебышева.

Экстраполяция этих аппроксимаций на отрезок [1, 2] показывает, что агрегаты P5(t) и S3(t) дают результаты значительно хуже, чем рациональная функция R1 2(t) . На рис. 2 приведены графики логарифмических относительных погрешностей lg| f (t) _ P5(t)| _ lg| f (t)| и lg | f (t) _ S3 (t) | _ lg | f (t) | , построенные на отрезке [1,2], а максимальные значения абсолютных погрешностей на отрезке [1,2] оцениваются следующими неравенствами:

max | f (t) _ P5(t)| < 45,3; max | f (t) _ S3(t)| < 0,581; max | f (t) _ R, 2(t)| < 0,266-10_8,

te[1,2] йе[1,2] te[1,2]

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

Таким образом, можно предположить, что алгоритм Кронекера-Чебышева дает рациональную аппроксимацию, близкую к наилучшей. Именно это свойство данного алгоритма показывает целесообразность его использования в задаче экстраполяции аналитических функций. В [8] на основе теории наилучшей аппроксимации непрерывных функций многочленами рассмотрены вопросы экстраполяции аналитических и квазианалитических функций вещественной переменной. Наилучшая рациональная аппроксимация позволяет получать не только аналитическое продолжение, но и определять особенности функции в виде полюсов.

Для сравнения некоторых методов экстраполяции в качестве первого теста используем известную функцию Рунге

причем они достигаются при t = 2, а значение функции Рунге в этой точке / (2) = 1/101

Рис. 1. Графики функции Рунге и ее аппроксимаций на отрезке [0,1]

Рис. 2. Графики логарифмических относительных погрешностей на отрезке [1, 2]

При рациональной интерполяции функции Рунге для t е [0,1] и n = 6 методом работы [4] с помощью встроенной процедуры «Rationallnterpolation» в пакете Maple требуется задание исходных данных в виде рациональных дробей. При отсутствии возмущений в этих исходных данных на равномерной сетке в результате получается функция Рунге.

При наличии малых возмущений в исходных данных, которые также записываются в виде рациональных дробей, процедура «Rationallnterpolation» находит рациональную функцию в форме R2 3(t), а алгоритм Кронекера-Чебышева с любым заданием исходных данных при соответствующем подборе параметра tol находит рациональную функцию в форме R1 2(t).

В качестве второго теста рассмотрим функцию

Г / \ ^ t

f(t) = tg—,

которая при четных целых значениях переменной t имеет простые полюсы.

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

Для исследуемой функции на отрезке [0,1] и n = 7 при выборе tol = 0,1 с помощью алгоритма Кронекера-Чебышева была найдена следующая рациональная функция:

„ /ч 0,02560196575 + 3,07346861

R1 2(t) = 2,

4,122049938 - 0,0570256245t -12

которая имеет полюсы: -2,058996232 и 2,001970608.

При выборе tol = 0 для этого случая найденная рациональная функция R^^(t) имеет полюсы: -2,049635455, 2,001695883 и 19,87598399.

Рациональная интерполяция методом работы [4] на равномерной сетке для отрезка [0,1] и n = 7 при использовании точной арифметики дает для исследуемой функции рациональную функцию R 33(t), которая имеет полюсы: -2,049950908, 2,001732621 и 19,89498128.

Максимальные значения абсолютных погрешностей на отрезке [1,1,9] оцениваются следующими неравенствами:

max | f(t)-R12(t)I < 1,823; max | f(t)-R33(f)| <0,131; max | f(t)-R33(01 < 0,135,

te[1,1,9] te[1,1,9] te [1,1,9]

причем f (1,9) = 12,70620482 .

Хотя на практике измеряемые значения функции f (t) обычно имеют ограниченное число верных знаков, приведенные тестовые примеры позволяют правильно оценить возмож-

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

Чебышевские узлы (2), используемые в алгоритме рациональной аппроксимации Кроне-кера-Чебышева, в процессе измерения также не могут быть заданы точно, но и в этом случае получаемые погрешности можно отнести к возмущению исходных данных. При обработке же результатов измерений на компьютере чебышевские узлы (2) должны задаваться с наибольшей возможной точностью, что особенно легко выполнить при применении современных математических пакетов аналитических вычислений, например, Maple.

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

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

Список литературы

1. Бабенко, К. И. Основы численного анализа / К. И. Бабенко. - М. : Наука, 1986. - 744 с.

2. Завьялов, Ю. С. Методы сплайн-функций / Ю. С. Завьялов, Б. И. Квасов, В. Л. Мирошниченко. - М. : Наука, 1980. - 352 с.

3. Бейкер, Дж. Аппроксимации Паде / Дж. Бейкер, П. Грейвс-Моррис. - М. : Мир, 1986. -502 с.

4. Beckermann, B. Fraction-free computation of matrix rational interpolants and matrix GCDs /

B. Beckermann, G. Labahn // SIAM Journal on Matrix Analysis and Applications. - 2000. -Vol. 22, № 1. - P. 114-144.

5. Бондаренко, Л. Н. Определение параметров передаточной функции средств измерений по значениям амплитудно-частотной и фазочастотной характеристик / Л. Н. Бондаренко // Датчики и системы. - 2004. - № 7. - С. 18-20.

6. Бондаренко, Л. Н. Методы идентификации в частотной области при наличии шума /

Л. Н. Бондаренко // Известия высших учебных заведений. Поволжский регион. Технические науки. - 2009. - № 2. - С. 113-123.

7. Бондаренко, Л. Н. О точности идентификации систем на базе RC-схем замещения /

Л. Н. Бондаренко // Измерительная техника. - 2011. - № 4. - С. 27-32.

8. Бернштейн, С. Н. Аналитические функции вещественной переменной, их возникновение и пути обобщения / С. Н. Бернштейн // Собрание сочинений. - М. : Изд-во АН СССР, 1952. - Т. 1. - С. 285-320.

Бондаренко Леонид Николаевич

кандидат технических наук, доцент, кафедра дискретной математики, Пензенский государственный университет E-mail: leobond5@mail.ru, dm@pnzgu.ru

Нефедьев Дмитрий Иванович

доктор технических наук, заведующий кафедрой информационно-измерительной техники, Пензенский государственный университет E-mail:iit@pnzgu.ru

Bondarenko Leonid Nikolaevich

candidate of technical sciences, associate professor, sub-department of diskrete mathematics,

Penza State University

Nefed'ev Dmitriy Ivanovich

doctor of technical sciences, head of sub-department of information and measuring technique,

Penza State University

УДК 621.317 : 519.67 Бондаренко, Л. Н.

Сравнение экстраполяционных методов обработки результатов измерений / Л. Н. Бондаренко, Д. И. Нефедьев // Измерение. Мониторинг. Управление. Контроль. - 2013. - № 2(4). - С. 3-9.