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

УДК 621.396.01

Е.Г. ЖИЛЯКОВ, д-р. техн. наук, БелГУ (г. Белгород),

Т.Н. СОЗОНОВА, БелГУ (г. Белгород)

О ВЫЧИСЛЕНИИ ВТОРЫХ ПРОИЗВОДНЫХ ПО

ЭМПИРИЧЕСКИМ ДАННЫМ НА ОСНОВЕ ЧАСТОТНЫХ

ПРЕДСТАВЛЕНИЙ

У статті розглядається новий метод оцінювання першої та другої похідних та інтерполяції за дискретними емпіричними даними. Інтерполяція та оцінювання похідних виповнюється у класі цілих функцій. При цьому для вибору конкретної функції, яка інтерполює, використовується варіаційний принцип мінімізації евклідової норми похідної, трансформанта Фур’є якої має фінітну область визначення.

In the article is considered a new assessment of the first and the second derivative and the interpolation of signals according the quantified empirical data. Interpolation and appraisement derivative are carried out in the class of entire functions. Whereat for choice of a concrete interpolating functions is used variatsionnyi principle of the minimization of the Evklid norm of the derivative, the coefficient of the Fure row which has the limited area of the determining.

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

Проблема заключается в том, что для анализа доступны только значения сигнала в некотором конечном наборе точек интервала изменений аргумента, то есть дальнейшей обработке доступен только некоторый набор чисел u(tk), где tk, k = 0,...,N - набор значений аргумента, причем для определенности полагаем, что

10 = 0; tN = T , (1)

а большему индексу соответствует большее значение времени. Для простоты рассмотрим случай эквидистантной дискретизации с шагом Дt, то есть имеют место равенства tk = k • Д. В дальнейшем зарегистрированное значение сигнала будем именовать отсчетом.

Анализ литературы. Проблема оценивания производных по эмпирическим данным исследовалась в работах многих авторов на протяжении сотен лет. В результате этих исследований предложено достаточно много методов оценивания производных [1 - 3].

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

ик = и (кА) = ик, к = 0,1,...,N. (2)

Далее приближенно полагают, что

и(г)(() И й(г)((), 0 < г < N . (3)

где г - порядок производной.

В качестве примера можно привести полином п-й степени Рп ((), например, в форме Лагранжа или в форме Ньютона и тогда оценки производных будут иметь вид:

и(г)(()и рПг)(() 0 < г < N. (4)

В ряде случаев наряду с приближенным равенством удается получить точное равенство (например, используя формулу Тейлора), содержащее остаточный член Я, характеризующий порядок погрешности численного дифференцирования:

,(г )(Л= р(г)

- п

Наиболее простыми и широко используемыми для приближенных вычислений производных являются конечно-разностные формулы:

- для первой производной

и'(0» й('+ АЬй('), (6)

А

где вычисления имеют первый порядок точности;

- для второй производной

и"(()И и(( ~А()_ 2и(() + и(( + А() (7)

~ (А()2 ’

где вычисления имеют второй порядок точности.

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

нулю при А ^ 0. В результате общая погрешность, возникающая при

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

иг(() = Р(г)(() + Я(А) 0 < г < N. (5)

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

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

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

вариационный принцип минимизации евклидовой нормы первой и второй производных.

Оценивание первой производной. С учетом изложенного

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

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

-да < ( < да, (8)

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

В основе дальнейших построений используется представление

(

и(() = и0 +| / (т)й?т, (9)

0

которое позволяет по производной

/(х) = ^ (10)

ах

вычислить интерполирующую функцию.

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

определения (8) и ограниченной евклидовой нормой, представимых в виде

/(х) = |F(ю^^ю, (11)

2% юеО

где 0 = [-02,-^1)^[^1,П2); <да; <да; Р(ю) - трансформанта

Фурье,

+да

^ (ю) =| / (х)е -“Ух. (12)

да

Выбор области определения О трансформанты Фурье может быть продиктован априорными сведениями о свойствах сигнала.

Подстановка представления (11) в правую часть (9) позволяет получить соотношение для вычисления интерполирующей функции на основе трансформанты Фурье производной

и(() = и0 +— [ Р(ю)“—(—/—^е 2 Ою, (13)

2п юеО Ю/2

так что интерполяционным равенствам можно придать вид

/юА(

1 Г ж Л51—(юА(//2) —>

Т" I Р(ю)-------гт;— е 2 Ою= V, / А(, (14)

2люеО юА( /2

где

V, = (и, -и0), 7 = 1,...,^ (15)

Можно привести достаточно много аргументов использования вариационного принципа

да

| / 2(х)ОХ=— | |Р (ю)|2 Ою = шт, (16)

-да 2л юеО

основанного на использовании равенства Парсеваля [2].

Один из аргументов заключается в целесообразности построения

функции с наименьшей в смысле евклидовой нормы производной скорости

изменения значений.

Другим важным соображением может служить необходимость повышения устойчивости вычислений к воздействиям случайных ошибок измерений (регуляризация).

Нетрудно показать [5], что искомое решение вариационной задачи (14), (16) представимо в виде

^ “ш(юМ /2)

Р(ю) = -----777^—е 2 ^ (17)

7=1 юА(/2

когда юеО и нулю в противном случае.

—* гр

Для вычисления вектора множителей Лагранжа р = (Р: ,...>Р N ) следует воспользоваться подстановкой представления (17) в левые части равенств (14). В результате нетрудно получить систему линейных алгебраических уравнений (СЛАУ)

Ар = V = (V!, ...,VN)т. (18)

Элементы матрицы

А = {% } (19)

определяются следующим выражением

К }=-1

1т -А-

1 °-2 “ш(хк / 2) “ш(х7 / 2)

- I -----------------------С0“

п0 (X /2)2

Iх (к - 7)

Ох, к, 7 = 1,..., N, (20)

где [01,О 2 ] = [А(О1, А(0 2 ] и [-О2 ,-01) ^ [01,0 2) - частотный интервал, в котором должна быть сосредоточена максимальная часть энергии интерполирующей функции.

Границы частотного интервала 01 и 02 задаются следующим образом:

— — (г -1) г

О = П , 02 = #2П , #1 =-М , #2 =— М , (21)

Л Л

где г = 1, ...,Л, Л - количество частотных интервалов.

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

N(^-^2 ) <п, (22)

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

О2 - О, “ш(хск / 2) “ш(хс7 / 2)

а 1к =-------------------------------------С0“

п (хс / 2)2

хс (к - 7) 2

1

< N' (23)

Здесь хс - некоторая средняя точка интервала интегрирования, которая в рассматриваемом смысле будет слабо зависеть от сочетания индексов к, 7. Поэтому значения элементов матрицы будут близки и ее определитель будет близок к нулю. Следовательно, решение СЛАУ вида (18) будет неустойчивым.

В этом случае для вычисления искомого р вектора следует использовать выражение

Р = А+у, (24)

где А + - псевдообратная матрица.

А+ = , (25)

если

А = дЬ£т. (26)

Q - матрица, состоящая из собственных векторов, так что

АО = LQ, (27)

где L = d7ag(X1,...,XM,0,...,0) - диагональная матрица, составленная из собственных чисел матрицы А, и - матрица собственных векторов, соответствующих ненулевым собственным числам,

01 = (Ч\,...Яз), А = d7ag('kl,...^J), ^1 >Х2 >...>^ >0. (28)

Так как определитель матрицы равен произведению ее собственных

чисел, то

з

аеА = Пі • (29)

і=1

Для определения ранга 3 используется неравенство:

з

Пі ^, (3°)

і=1

где є - задаваемое число, конкретно є > 1° -7 •

Используя тригонометрические тождества, подынтегральную функцию в (20) можно представить в следующем виде

ът(хк / 2) 8Іп(хі / 2) 008|(х / 2)(к - і)] _ 1 + со8[х(к - і)] - 008(хі) - 008(хк)

(х /2) 2 х2

-•(31)

Отсюда и из (22) следует, что если выполняется неравенство

^2 — = NАt(^2 — ) — 2 л, (32)

то хотя бы одна из тригонометрических функций числителя (20) при

изменениях г,к = 1,...,N, на интервале интегрирования претерпит изменения на полном периоде. При этом гарантировано, что наименьшее из собственных чисел, а следовательно и определитель матрицы А будет заметно больше нуля, то есть матрица А будет неособенной. В этом случае получаем точное решение системы (18)

Р = А~^ . (33)

После решения СЛАУ (18) для вычислений оценки производной и интерполирующей функции можно воспользоваться соотношениями (11) и (13), куда следует подставить представление (17). Предполагая, что М -количество интерполируемых значений внутри одного интервала

интерполяции, нетрудно получить вычислительные формулы [5].

Формула для вычисления интерполирующего вектора имеет вид

и = и0е + СА+V, (34)

где V = (у1, v2,...,vN) - вектор разностей; матрица

С = К-}, (35)

элементы которой определяются следующим образом:

а

К !=11

іт А.

8Іп(хк / 2М) бш(х/ / 2)

>сй/=т I---------------------008

(х /2)2

(36)

■4

і = 1,..., N; к = 1,..., ЫМ.

Для вычисления оценки производной на заданном наборе точек используется формула

1 = ВА+V , (37)

где матрица

задается выражением

П- г

х ( 2к

К }= - /

тг -А

01

1 г 8Ш(х/ / 2)

008

— г

2 I М

л х/2

О}

г = 1,...,N +1; к = 1,...,NM.

йх,

(39)

Оценивание второй производной. Для оценки второй производной по эмпирическим данным предлагаются два метода:

1. Второй метод предполагает вычисление второй производной на основании первой. Вычисление первой производной осуществляется по выражению (35). Затем при оценке второй производной в качестве исходных данных принимаются вычисленные значения первой производной. Оценка второй производной осуществляется по выражению (35).

2. Вычисление второй производной без предварительного вычисления первой.

Рассмотрим теоретические основы второго метода.

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

I

и(7) = и0 +| / (т)йт, (40)

о

где /(х) - первая производная, в свою очередь представимая через вторую производную ф(9);

X

/ (х) = /о +1 ф(9)йх. (41)

о

Для определения второй производной предлагается использовать класс непрерывных вещественных дифференцируемых функций с областью определения (8), представимых в виде

ф(9) = — |ф(ю)е]ю9йю . (42)

2л юеО

Следовательно, выражение для первой производной при подстановке (34) в (33) будет иметь вид

х 1 1 е ^'юх — 1

/(х) = /о +1 | ф(ю)е-/ю 9йюй9 = /о + — | Ф(ю)------:---йю . (43)

0 2л юеО 2л юеО }Ю

Подставляя выражение (43) в правую часть (40) и допуская изменение порядка интегрирования, получим соотношение для вычисления второй производной

(Ґ) = и0 + /0Ґ + — І Ф(ю)І(ю, Ґ^ю ,

2л юєО

где І (ю, Ґ) =

Ґ

Г х^х = —, ю = 0; 1 2

0

Ґ /юх і Г ^ -1

0

/Ю

dх = — /ю

?/ЮҐ -і

:-Ґ

/Ю

(45)

ю Ф 0.

Таким образом

и(ґ) - и0 - /0Ґ = — | Ф(ю)І (ю,Ґ)dю

2л юєО

Потребуем выполнение условия

1 ||ф(ю)|2 dЮ = ті

шт

при ограничениях

’г = и(іД) - и0 - / (іД) = — | Ф(ю)І (ю,іД^ю .

2л юєО

(46)

(47)

(48)

Поставленная задача сводится к поиску решения вариационной изопериметрической задачи, определяемой условием (47) и ограничениями вида (48).

Нетрудно показать, что искомое решение представимо в виде

N

ф* (ю) = ХРі7 (ю, іД);

і=1

N

Ф(ю) =£р*/* (ю, кД);

к=1

N 1 г N

^$к— 11*(ю,кД)І(ю,іД)аю> = ^Рк(сс)к ;

к=1 2% юєО к=1

(асі)ік = — 11(ю,іД)І* (ю, кД)dю .

(49)

(50)

(51)

(52)

юєО

Подстановка сюда представления (45) дает:

(сс) ік =

(ік)2 (Дґ)3 Іл(^2 - 41)

2л4

4N

N

м

т=1

СОБ

[т5(і - к)]5Ш(т5і) 8і"(т5к 1 - со5(т508іП(т5і)

(т5і) (т5к) (т5і)

- оо8(т5к) 8ІП(т5к 1 +1

/ т 2 }.

(т5к)

—► гр

Для вычисления вектора р = (Р1,...,Р N) следует воспользоваться

подстановкой представления (50) в равенства (48). В результате нетрудно получить систему линейных алгебраических уравнений

ААр = її,

(53)

где wi = иі - и0 - /0іД ;

/0 - значение первой производной, определенное ранее

АА = {(аа) й }. (54)

После решения СЛАУ (53) для вычислений оценки второй производной можно воспользоваться соотношениями (42) и (44), куда следует подставить представление (51). В результате получаем вычислительные формулы

1 М 1

ф(0) = — Гф(ю)е'ю0<*о = Урк — [I*(ю,МОе^ю . (55)

2^ ^ з

' юєО к=1 2л юєО

Если теперь положить, что 9 меняется дискретно

0 = пД9; Д9 = — ; г =

ДҐ

=----; п = 0,1, ...,(Ы-1)М .

М п Д МД м где N - первоначальное количество отсчетов иі, то будем иметь

N ( Ґ - Л „;„2

ф(9 п)=Д £к 2р к 2 к=1

42-141 + „і ш5| 2тП51 +

2Л' ЯІ I М ! (т5к)2

біпі 2т — 5

М

(т5к)

008( т5к )8!П(т5к2 -1

(т5к)

Следовательно

где

ф = ввр,

ВВ = {(ЪЪ) к1};

(56)

(57)

= п

+

(bb) nk = — k

nk 2N

q 2- qi

2N

m=1 I

X] cosI 2m — 5

n ^ sin (шЪк)

M

(ш5k )2

sin| 2m — 5

M

ш5k

/ 4sm(m5k)

cos (даек )----------1

m5k

Таким образом, искомое решение представимо в виде

At

и(t) = uo + fot + — Xк2 xp 2 к=1

N

+ Z

m=1

_ s . sinl 2m — 5

I. n Л sin (m5k) ^ M

sl 2m — 5 I------------------

M

(m5k)

(m5k)

q 2 - q1 2N

sin(m5k)

cos<mck)-------------1

(m5k)

(5S)

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

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

Рунге | у =--1—- I, функции у = cos X , у = cos 10х.

^ 1 + 25х2)

2

Значения аргумента х = t -1, t е [0; 2], At = — .

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

В качестве начальных данных использовались вычисленные значения модельных функций и их первых и вторых производных, взятых с шагом 5At, то есть каждое пятое значение модельной функции (M = 5). Оценка первой и второй производных осуществлялась с шагом дискретизации At, то есть осуществлялось вычисление производных в четырех точках между соседними значениями начальных данных.

Вычисление первой и второй производных осуществлялось двумя описанными выше методами:

1. На основе выражения (58).

2. На основе выражения (35).

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

N

2

+

+

+

к

so

X (У - 2,)2

где Уг - значения производной модельной функции; г = 1,...,ИМ; Zi -значения производной аппроксимирующей функции.

Результаты эксперимента приведены в таблице.

Таблица

Функция Количество начальных данных, N Погрешность оценки первой производной, 8 Погрешность оценки второй производной, 8 (1 метод) Погрешность оценки второй производной, 8 (2 метод)

1 у 1 + 25х2 11 0.1136 0.2779 0.3834

21 0.0106 0.0425 0.0769

31 0.0029 0.0078 0.0275

41 0.0014 0.0079 0.0258

51 0.0016 0.0058 0.0296

у = 008 X 11 0.0354 0.0271 0.0961

21 0.0150 0.0229 0.0632

31 0.0087 0.0217 0.0504

у = 008 10 X 21 0.0083 0.0131 0.1523

31 7.3537-10'4 0.0224 0.0752

41 8.9917-10'4 0.0202 0.0531

Из приведенных данных видно, что оптимальным является метод, в котором для оценки второй производной в качестве исходных данных используется оценка первой производной, рассчитанная по выражению (35). Данный факт может быть объяснен следующим образом. При расчете второй производной по формуле (58) может возникать погрешность за счет неточного задания значения первой производной /0, так как наибольшее отклонение первой производной от ее "истинного" значения наблюдается вблизи концов интервала области определения.

Действительно, оценку первой производной в точке и0 можно

представить выражением /0 = /0 + у, где у - погрешность оценки первой производной, /0 - ее "истинное" значение. Таким образом

/0 = /0гА1 + у/А/, г = 1,..., N. При достаточно большом значении г увеличивается значение слагаемого у1Д, а, следовательно, погрешность оценки первой производной.

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

эмпирическим данным. Предлагаемый подход основывается на классе функций с финитными областями определения трансформанты Фурье (финитными спектрами) и в этом он близок к условиям применимости интерполяционной формулы Котельникова-Шенона.

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

Список литературы: 1. Ланцош К. Практические методы прикладного анализа: Справ. рук. - М.: Физматгиз, 1961. - 524 с. 2. Хургин Я.И., Яковлев В.П. Финитные функции в физике и технике. -М.: Наука, 1971. - 408 с. 3. Колесников А.П. Введение в численный анализ: Учеб. пособие для физ.-мат. фак. ун-тов. - М.: Из-во РУДН, 2002. - 218 с. 4. Де Бор К. Практическое руководство по сплайнам - М.: Радио и связь, 1985. - 304 с. 5. Смирнов В.И. Курс высшей математики: Учеб. пособие для мех.-мат. и физ.-мат. фак. ун-тов. В 5-ти т. - М.: Наука, 1974. - Т. 4. - Ч. 1. - 336 с. 6. Жиляков Е.Г., Мисливец И.Ю., Созонова Т.Н. Вариационный метод оценивания производных и интерполяции сигналов по эмпирическим данным // Вестник ВГУ, серия "Системный анализ и информационные технологии". - Воронеж, 2006. - Вып.2. - С. 70-73.

Поступила в редакцию 13.09.2007