CC BY
94
УДК 517.925.54:517.962.27/.8
О ДИСКРЕТИЗАЦИИ ЛИНЕЙНЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ
А.О. Егоршин
Рассмотрены некоторые вопросы получения дискретного описания дифференциальной системы (ДС) на равномерной сетке. Рассматриваются ДС в виде системы п линейных обыкновенных дифференциальных уравнений первого порядка с постоянными коэффициентами или одно уравнение п-го порядка для наблюдаемого функционала состояния ДС. Изучаемые вопросы дискретизации важны для задач вариационной идентификации и аппроксимации динамических процессов моделями этого типа в конечном интервале. Дано сравнение аналитического равномерного (на основе теоремы Гамильтона-Кэли) и локальных методов дискретизации: на основе разделенных разностей и с помощью интерполяции выборок из п +1 отсчетов многочленами Тейлора степени п. Получена общая формула локальной дискретизации, прозволяющая сравнивать ее разностный и интерполяционные методы. Показано с использованием свойств обратных матриц Вандермонда, что в полученной общей формуле локальной дискретизации ее интерполяционному методу соответствуют (п + 1)-матрицы Тейлора (из коэффициентов многочленов Тейлора), а разностному — (п + 1)-матрицы Паскаля (из чисел треугольников Паскаля).
Показано, что невырожденность матрицы наблюдаемости ДС на сетке есть необходимое и достаточное условие как для аналитической дискретизируемости, так и для приведения дискретной системы (описания ДС сетке) к каноническому фробениусов-скому виду. Он эквивалентен одному обыкновенному разностному уравнению для наблюдаемой переменной с постоянными коэффициентами. Это уравнение есть основа известного вариационного метода идентификации. Показано, что интерполяционный метод локальной дискретизации есть первое (линейное) приближение формулы равномерной аналитической дискретизации. Показано, что нулевое приближение ее не зависит от коэфффициентов ДС и есть вектор коэффициентов п-й разности. Показано также, что нулевое приближение матрицы наблюдаемости ДС н и матрицы наблюда-мости полиномиальной системы у(п) =0 на сетке есть п-матрица Тейлора.
Ключевые слова: вариационная аппроксимация и идентификация, дискретизация дифферециальных уравнений, аналитическая дискретизация, линейное приближение, теорема Гамильтона-Кэли, локальная дискретизация, многочлены Тейлора, матрица Вандермонда, треугольник Паскаля.
Введение
Дискретизацией дифференциального уравнения называем получение разностного уравнения, описывающее отсчеты решения дифференциального уравнения на заданной Л,-сетке. Мы говорим об однородных уравнениях.
Вопрос о дискретизации становится актуальным при решении некоторых задач идентификации (оценки коэффициентов) дифференциальных уравнений, для которых удовлетворительные методы решений существуют лишь для их разностных аналогов. К таким задачам относятся, в частности, вариационные задачи идентификации [1].
Отметим два подхода к дискретизации дифференциальных уравнений: локальной — с помощью конечных разностей (он в некоторых случаях применим и для уравнений с переменными коэффициентами) и равномерной. Один из них — аналитический метод дискретизации на основе теоремы Гамильтона-Кэли [2]. Он дает разностное уравнение для отсчетов на Л,-сетке любого решения ДУ на любом интервале.
И дифференциальные уравнения (ДУ), дискретизации которых посвящена статья (обыкновенные линейные уравнения с постоянными коэффициентами), и получаемые разностные уравнения (РУ) есть, с точки зрения идентификации и математического моделирования, некоторые (наиболее простые в нашем случае) модели реальных динамических процессов. Поэтому для обозначения исходных данных (идентифицируемых реальных процессов) используются обычные латинские буквы (например, y(t)), а их описания с помощью моделей отмечаем надсимвольным знаком, в данном случае, y(t).
Условимся и о других обозначениях. Конструкции вида {xi }m будут обозначать вектор-столбец, | Xj |П или |xi, ...,xn| — вектор-строку, а {xj }m]n — матрицу (i — номера строк) множества элементов вида x. Они могут быть и скалярами, и строками, и столбцами, и матрицами. Если элементы xi в {x^^ — (n — l + 1)-векторы-строки, а xj в |xjЩ или в |xi,..., xn| — (m — k+1)-векторы-столбцы, то эти конструкции определяют матрицу {xj}m]n-Верхний индекс * у скаляров обозначает комплексное сопряжение. У векторов и матриц, в том числе и действительных — он обозначает еще и транспонирование.
Нумерация L + 1 узлов ti, i = 0, L, h-сетки с L ячейками интервале наблюдения Ту, T = Lh ведется с нуля. Поэтому нумерация и отсчетов y*, i = 0, L решетчатой функции y* — строки y* = |y*|o, и компонент (L + 1)-вектора y = {yi}° € E = E0+1 ведется также с нуля. Последовательность длины L + 1 из l2 [L] компонент вектора y называется далее реализацией (отсчетов yi, i = 0, L).
С нуля ведется и нумерация других элементов в используемых евклидовых пространствах E = E0+1, En+1, En. Первая позиция этих и других векторов считается нулевой. Это согласуется с нумерацией компонент в векторах производных — векторах состояний дифференциальных уравнений.
С нуля нумеруются и орты (столбцы) в евклидовых пространствах. Канонические орты обозначаются так: в E — это ej = e.j/L+1 = №,j}о, j = 0, L, в En+1 — это ej/n+1, j =
0,n ив En — это e j/n, j = 0,n — 1. В частности, eo/n = |1, 0^.-11* и en_1/n = |0^-1, 1|* — есть соответственно первый и последний n-й орты в En. Здесь и далее 0^ — нулевой к-вектор-столбец, а 0^ — нулевой k-вектор-строка. Аналогичные орты ej = (•, ej ), ej/n+1 =
( ', ej/n+\) , ej/n = ( ■, ej/n) в сопряженных пространствах — E1 = E/0+1, E/n+1, E/n — будут обозначаться как e*, e* . ,, e* соответственно.
j ]/n+V j/n
1. Аналитическая равномерная дискретизация
1.1. Модели вариационной идентификации
Для получения аппроксимации функций y(t) € L2 минимизируется выпуклый функционал метода наименьших квадратов (МНК) —
Jc = ||y(t) — y(t) ||0,2, где L2 = L2[It ], It = [0,T ]. (1)
Здесь y(t) € C(n)[Iy] - переходные процессы линейных обыкновенных ДУ с постоянными коэффициентами -
n
^y{i)(t)a* = 0 = (Wt,a)En+1, где Wt = (y)t = {^(t)}n, t € It. (2)
i=0
Задача аппроксимации (1), (2) при неизвестных коэффициентах a уравнения (2) является задачей их вариационного оценивания (идентификации) [1].
Такая постановка задачи аппроксимации (функционал (1) и ДУ (2) в качестве условий его минимизации) есть обобщение задачи МНК-аппроксимации функций полиномом степени n — 1 на конечном интервале. Это легко увидеть, если указанную известную задачу сформулировать аналогичным образом: минимизировать Jc из (1) при условии, что y(n) = 0.
Для аппроксимации решетчатых функций у = {у}о £ Е^1 = Е на Л,-сетке
/ь = /ь(Ь) = {¿г = гЛ, : г = 0, Ь, Т = ЬЛ,, Л, > 0} С /т = [0, ЬЛ,] (3)
мы используем функционал МНК в Е (или, что то же, в 12[Ь])
о
^ = ||у - у||| = ^ ^ - ^|2’ (4)
о
а в качестве условий его минимизации — РУ вида, аналогичного (2) (линейное, однородное с постоянными коэффициентами) для реализации у £ Е.
Пусть к = 0,Ж, N = Ь — п, у = у(^), ^ £ 4, г = 0,Ь,
п (5)
тогда ^=0 ук+г«* =0 = (у ,а)Еп+1 , где у = [у]; = {у+г}о и а„ = 0.
Векторы вида V подряд следующих отсчетов из их последовательности у £ Е = Ео+1
называются далее выборками. Векторы V;, к = 0, N из (4) являются (п + 1)-выборками — (п + 1)-векторами из Ео+1.
Для вычисления и прогнозирования решения РУ (5) важна его рекуррентная форма:
о—1
у;+п = - ^ у;+га* = —(у;, а>Еп, где ^ = [У]; = {у;+г}о—1 (6)
г=0
— вектор состояния модели динамической системы, описываемой уравнением (5). Вектор а * назван прогнозирующим вектором коэффициентов РУ (5), (6).
Уточним эти обозначения для коэффициентов РУ (5) и ДУ (2):
а * = | а0,...,а**—1,а0|, если ао = 1, то
>к >к >к >к 1 а = 1 а0,..., ап—1, ап|, если ао = 1, то
Запишем ДУ (2) в эквивалентной канонической форме линейной системы п уравнений
первого порядка с фробениусовской матрицей перехода — еп— 1/оа*. Здесь /^ =
{¿¿+1а-}0— 1 — матрица сдвига, частично изометрический оператор в Ео.
Вектором состояния динамической системы с матрицей перехода является «укороченный> п-вектор го производных наблюдаемого функционала у = (го*, в0/га)еп: го* = {У(г)}о— 1 • Пусть в ДУ (2) ао = 1. Тогда, если
а = | а0,..., ап _ 1,11 = | а , 11,
(7)
^ 1 —^ 1
a = I a0,..., an _ 1, 11 = I a , 11.
Af
—an
1 n_1
—a1, ■ ■ ■ , —аП_1
AF (a) In era_1/raa ,
(8)
то на It : wt = Afwt, yt = e0/nwt = (гоt,eo/n)E„
Вследствие теоремы Гамильтона-Кэли, для переходной матрицы B = exp (Af h) дискретной системы го^+1 = Bw& имеем равенство ^П a*e0/nB = 0, из которого и вычисляется вектор а * = |а0,..., а**_1,1| = |а*, 1| коэффициентов РУ (5), (6), если положить an = 1.
0
1.2. Матрица наблюдаемости на равномерной сетке
Пусть имеется система п дифференциальных уравнений с наблюдаемым значением у(£) линейного функционала в * = в! = (•,^)Ёп € Е/га на п-векторе состояния ж(£):
ж(і) = Аж(і), у(£) = й*ж(£) = (ж(і), й)Еп , і Є /т = [0, ЬЛ,].
(9)
Дискретный аналог системы (9) на сетке /^ (3) для Ск = зс(кН), к = 0, Ь:
Ок+1 — ФхСк, ук+п — й ОСк+п — (З'к+Ъ Ф
* п—1
Еп
, Ф = ехр(ЛЛ).
(10)
Определим операторы наблюдаемости ^к = {в *Фк+г}П—о1, к = 0, N. Эти операторы определяют связь между состояниями систем (10) и (6). Образуем линейную функцию , к = 0, N значений у(^), г = 0, £ наблюдаемого функционала (•, в) в (9) на сетке :
Еф‘
г+к а*
,г+к\
С , а> Со.
J г=0 / ЕП+1
(11)
Здесь
Ск
^к
Ук+п
п— 1
)к+Л п о
= \ й* Фк'^ ОСо
^к й * Фк+п
Хо, где
(12)
^к = {й*Фк+г}0 '= КФ^П 1 Фк = ^оФк, а ^о = {й*Ф^-1 = ^
— матрицы наблюдаемости дискретной системы (10) или, точнее, матрица наблюдаемости дифференциальной системы (9) на Л,-сетке интервала /у.
Выборки Ук есть п-векторы состояния эквивалентной уравнению (6) следующей дискретной системы с канонической матрицей Фу фробениусовского типа (см. (8)).
Пусть
Фу = Ф^ (а*) =
0
ао
/п—1
-а1, ■ ■ ■ , -аП—1
— /1 — в і / а *
/п вп—1/па
тогда ^к+1 = ФуСк, Ук+п = еП—1/п^к+1 = ^Ск+1, вп—1/п
Еп
(13)
1.3. Равномерная аналитическая дискретизация
Рассмотрим две задачи анализа системы (9) на сетке /^ с помощью матрицы ^ наблюдаемости дифференциальной системы (9) на этой сетке. Обе эти задачи приводят к формуле равномерной аналитической равномерной дискретизации любого решения системы (9). Первая — на основе теоремы Гамильтона-Кэли, вторая — с помощью преобразования системы
(10) к каноническому виду (13).
Теорема 1. Обратимость матриц ^к, к = —п, N наблюдаемости системы (9) на сетке /^ есть необходимое и достаточное условие для решения следующих задач.
я. Равномерная аналитическая дискретизация: получение коэффициентов РУ (6) для вычисления отсчетов любого решения ДУ (1) на сетке /^ интервала /у.
Ь. Получение любого решения системы (9) на сетке /^ для наблюдаемой переменной с помощью РУ (6) или дискретной системы (13) с фробениусовской матрицей Фу.
Доказательство. Во-первых, отметим, что невырожденности матриц наблюдаемости ^к, к = —п, N необходима невырожденность матрицы Ф дискретной системы (10).
а. Пусть а * — характеристические числа матрицы Ф системы (10). Тогда, по теореме Гамильтона-Кэли, Фг+ка* = 0 для всех к = 0, N. Это тождество, во-первых, дает из
п
(11) разностное уравнение (5). Во-вторых, поскольку вектор ж° произволен, то, положив а« = 1, из тождества Гамильтона-Кэли получаем, для прогнозирующего вектора а* алгебраические уравнения при любом конечном к. Для их решения необходима и достаточна невырожденность матрицы V:
*Л «+к — _ 7^7*
^*Ф
—№
^*Ф
-а*^ Ф
к
а* = - ^*Фга^-1 = - ^_«.
(14)
Ь. От системы (10) с переходной матрицей Ф общего вида не трудно заменой переменных V = ^°ж перейти к каноническому описанию (13) системы (10). Эти описания неразличимы с точки зрения наблюдаемой переменной — значений функционала у = ^*Ж. Определим:
В * = кфч«-1, ь* = ^*ф° = ¿*, ь« = ^*Фга.
(15)
Делая теперь в (10) указанную замену переменных ж = V, получим систему с вектором состояния 5к = [V] к = (Уг+к}о-1, к = 0, Ж, с матрицей перехода ФЕ и с наблюдателем ^:
Ф^ = V Ф^-1, Ук+п = 41 ук+1 = ¿*^Г-1гаг;к+1 = ¿*^Т«г;к = -а*^- (16)
Последнее равенство показывает, что рекуррентное уравнение (6) может рассматриваться и как наблюдатель, и как предсказатель.
Покажем, что Ф^ в (16) есть матрица Фробениуса, а наблюдатель = е«_1/га из (13). С учетом определений (15) можем записать:
ф^ = (V ф)^°-1 =
в*
ь*
V-1 = г0 =
в*
Ь*
*
ь0
в*
_1
в*
Ь*
ь(-1),в(-1)
(17)
Из последнего равенства в (17) очевидны тождества: В*В( 1) = /«_1 — единичная п — 1-матрица, а В*Ь° 1) = 0. Отсюда видим, что матрица Ф^ в (17) — матрица Фробениуса из (13), если обозначить Ь«^0-1 = —а*. Учет (15) приводит к формуле —а* = ^*^_« из (14).
О каноническом наблюдателе. Пусть В* = (^*Фг_п+1}П_2. Тогда из (16) имеем:
В*
_1
в(-1),^(-1)
= |0«-1,1|,
(18)
что и требовалось. Необходимость и достаточность обратимости матрицы V очевидна. □
2. Локальная дискретизация и конечные разности
2.1. Дискретизация уравнения для полиномов
Широко используемая альтернатива равномерной аналитической (точной) дискретизации — локальная приближенная дискретизация. Она осуществляется при помощи операторов = Дг/Лг, г = 0,п, конечных разделенных разностей ук, к = 0, N в пределах г + 1 ближайших к текущей точке к отсчетов ук реализации у. Такая локальная дискретизация основана на приближении решения полиномами порядка г = 1,п (п — порядок уравнения) на интервалах длины гЛ.
В п. 1.1 было отмечено, что модель (2) для аппроксимации функций (сигналов) есть обобщение модели полиномиальной аппроксимации у](”) = 0, где £ € /у. Модель (5) есть аналогичное обобщение модели:
ДгаУк = 0 = [ДТ?к, где [Д«]* = |(—1)"_' С«|«, С« = п/ !(п — ;)! (19)
для полиномиальной аппроксимации отсчетов на сетке /^ для к = 0, N. Здесь через Д« обозначен оператор разности порядка п, а через [Д«]* — (п + 1)-вектор-строка его коэффициентов С«. Это вектор-строка коэффициентов РУ (19).
Прогнозирующую часть вектора коэффициентов рекуррентной формы РУ (19) отмечаем, как и ранее (см. (3),(5)) надсимвольной чертой: [Дга] = |(—1)”_гС«|« . Это значит, что
отсчеты ук, к = 0, N полинома порядка п определяются рекуррентной формулой, аналогичной формуле (6):
ук+1 = —[Дга] "«к = —(^к, [Дга])еп . (20)
Далее мы опишем способ локальной дискретизации на основе интерполяции выборок отсчетов заданной функции на интервалах длины пЛ полиномами порядка п. Этот метод локальной дискретизации основан на интерполяции многочленами Тейлора и поэтому назван тейлоровским. В п. 2.4 мы покажем его связь с упомянутым выше методом и операторами разделеных разностей. Будет также показано, что локальная тейлоровская дискретизация есть линейное приближение к равномерной аналитической на основе теоремы Гамильтона-Кэли.
2.2. Матрицы Тейлора и Вандермонда
Обозначим через /т(^) минимальный интервал (длины Лт), содержащий выборку -и(т) из т + 1 подряд следующих отсчетов реализации у. Если индекс и аргумент т опущены, то т = п, то есть длина выборки есть п + 1 или п как определено ранее для стандартных выборок V и V в уравнениях (5) и (6).
Запишем для отсчетов выборки Vk(т) = {ук+г}™, где к = 0, Ь — т представление с помощью их интерполяции многочленами Тейлора порядка т. Пусть ¿к+т = кЛ+г°Л — некоторая точка в окрестности точки ¿к выборки 5к В этой точке значение т производных полинома, интерполирующего отсчеты выборки, и аппроксимируемой по этим отсчетам функции совпадают. Представление для отсчетов ук+г, г = 0, т с помощью такого интерполирующего полинома может быть записано так:
т Л
^у0)(£к + т) ■ (и(г — !1 —►
«=о )
Vk(т) = Т« (г°) ■ ад(т) (¿к + т), где Т™ (г°) = (Т*«}т = { (к (г — г°))/{;!}^. (21)
В формуле (21) через -ш(т)(£к + т) = (у(«) (¿к + т)}т° обозначен (т + 1)-вектор производных функции у(£) € С(га) в точке ¿к + т, а Тт(г°) — квадратная (т + 1)-матрица коэффициентов т-степенной интерполяции т + 1 отсчетов выборки Vk(т).
Матрицы Т вида (21) мы называем матрицами Тейлора или тейлоровскими матрицами. Если индекс т = п, то он опускается. Если т = п — 1, то используется надсимвольная черта: Тга_1 = Т. Опускается и аргумент г°, если его значение не существенно. Основные его значения будут 0, п и матрицы Т(0), Т(п) для граничных точек интервала /«(V).
Пусть для простоты т = п и пусть д = { дг}«. Матрица Тейлора Т из (21) может быть представлена в виде:
Т(г°) = W(д) ■ ^{Л/,;!}« = ЖД где Ж(д) = {Ж«}« = {д«}« = {(г — г°)«}«. (22)
Так определенную матрицу Ж = Ж(д) назовем матрицей Вандермонда. Ее компоненты: п + 1 степеней ; (строки) п + 1 чисел д = {дг}« (столбцы).
^(т) = {ук+^Г=° =
Лемма 1. Последняя строка т” 1)+ обратной матрицы Тейлора T 1 есть [Д”]*/h” — коэффициенты оператора D” = Дга/Л° разделенной разности порядка n, где [Д”]* = |(—1)0-jCj|”_0 (19), а Cj = n!/(n — j)!j! — коэффициенты бинома (s — 1)” =
5^0 Cj( —1)”-j sJ в порядке степеней s — оператора сдвига: syi = Уг+1-
Доказательство. Этот факт следует из свойств матрицы Вандермонда и ее обращения, поскольку TT-1 = WDD-1W-1 = WW-1 = I. Из последнего следует, что столбец W-обратной матрицы W-1 Вандермонда W есть коэффициенты полинома
n / n / n p(j)(A) = Y1 p(j)rAV= П(A — q^/^ , где ^ = П (qj— q^ j = 0>™
r_0 0,i_j 0, i_j
в порядке степеней Л. Поскольку коэффициент P(j)” = 1, то нормирующие множители
1/^j есть компоненты последней строки г,! 1)+ обратной матрицы Тейлора T-1, так как в соответствии с равенством WW-1 = I должно выполняться для всех i,j равенство ki |0_0 • {P(j)r}”_0 • Я-1 = Æij, где Æÿ — символ Кронекера.
Далее. Из определений (21), (22) и выборок Vk, k = О, N следует, что числа qi + ¿0 = 0, n, то есть они есть отрезок натурального ряда. Следовательно, сомножители qj — qi в формуле для множителя ^j — две части указанного отрезка, разделенные числом j = 0, n. Поэтому
П
& = П(j — ¿) =(n—j)!j !(—1)o-j, j = м.
0, i_j
Учтем снова представление (22), а именно, что T-1 = diag {j!/hj}” • W-1(q) (22). Последнее дает в последнню строку числа n!/(^jh”). Это и требовалось доказать. □
Следствие 1. Последняя (n + 1)-строка т(-1)* матрицы Тейлора T-1(i0) не зависит от точки ¿0. Число т(-1)*v есть оценка n-й производной в разложении Тейлора функции y(t) относительно средней точки интервала In. Тогда в (21): т = т0 = hî0n = hn/2 ^ ¿0 = n/2.
Доказательство. Первое утверждение вытекает из леммы 1. Второе следует из того, что у n-полинома n-я производная — число т(-1)*v — постоянна. Это число —
т(-1)*Vk = (vfc, [Д°])е„+1 /h” = (|(—1)”-jC”|0 • Vk)/h”
в силу симметричности коэффициентов вектора |( —1)”-jC”|” естественно считать оценкой производной y(n) (tk + т0) функции y(t) в средней точке т0 = hn/2 интервала In(vk) по ее разложению Тейлора на этом интервале относительно средней его точки т0. □
2.3. Локальная тейлоровская дискретизация
Лемма 2. Общая формула локальной дискретизации имеет вид
а* = а* • T-1 • с, (23)
где с — множитель, зависящий от вида матрицы T-1 (матрицы разностей) и нормировки вектора коэффициентов РУ (например, к единичной длине [3, 4] или к значению а0 = 1).
Доказательство. Из (23) следует, что матрица Тейлора невырождена, если не вырождена соответствующая матрица Вандермонда. Как известно, определитель последней отличен от нуля, если все числа q различны.
Пусть в (21) т = п. Тогда из полученной в (21) невырожденной системы алгебраических уравнений, получим выражения вида ад = Т-1^ для (п + 1)-векторов производных ад, определяющих ДУ (2). Подставляя эти выражения в уравнение(2) для значений £ на сетке , получим формулу локальной тейлоровской дискретизации:
= Т ^, к = 0, N —►
(ад^, а)£п + 1 = (Т % , а^> Еп+1 = ,Т а) Еп + ! = (Ук, а)йп + 1
а* = а*Т-1,
(24)
что и требовалось. □
Представим матрицу Т в окаймленном виде с учетом ее общего представления (21):
Т
Т=
* =
т* = |т*, $| , где $ = (Л(п — ¿0))7п! ^ $(п) = 0,
(25)
а * = {(Л (г — го))П/п!}”=0, т* = КЛ(п — ¿о^/і! |"=0•
і и\ 1”_1 і=о •
Следствие 2. Если в лемме 2 матрица Т есть матрица Тейлора из (21), а нормировка вектора коэфффициентов а РУ (5) осуществляется к значению ап = 1 (как в уравнении (6)), то нормирующий множитель в формуле (23) леммы 2 равен:
_1_;
—1_;
с = є ■ є 1, где є = $ — т*Т % є = 1 — а*Т "*• Доказательство. Используем формулу Фробениуса [2, 5]
(26)
Т_1____
Т 1 0
0 0
+ ¿є
-1г*,
где
г = — _ 1—Т 1* , )* = — _1 —т*Т 1,1
1
(27)
и є = (т *) = ()* *) = $ — т *Т 1*. Умножая это слева на а*, получим требуемый результат.
Лемма 3.
____і ___________.
т *Т 1 = — [Д”]'
є = $ — т *Т 1* = $ + [Д”]~ * =
Доказательство. Эти равенства следуют из представления последней строки матрицы Т 1 в формулах (26), (27) и из леммы 1. В соответствии с этими результатами
є-1)* = є-11 — т*Т \ 11 = [Д”]*/Лга = |[Д™] , •
Следствие 3. В формуле (23) тейлоровской дискретизации множитель с = Л”/є.
□
Теорема 2. При тейлоровской локальной дискретизации формула для прогнозирующих коэффициентов имеет вид:
а
= а * (єТ) 1 Л” + [Д”] *, где є = 1 — а*Т 1*, а * ={ (Л(г — ¿о))7/ п!}
\га/ і -і га—1 о
Доказательство. Этот факт следует из леммы 2, из ее следствия, а также из равенств, установленных в лемме 3. □
т *
2.4. Конечные разделенные разности
Наиболее известен метод дискретизации с помощью замены производных конечными разделенными разностями порядка от 0 до п.
Лемма 4. Разностная дискретизация эквивалентна использованию в общей формуле локальной дискретизации (23) нижнетреугольной матрицы разностей Т-1. Строка г, г = 0, п ее треугольника с точностью до множителя К есть строка г знакопеременного треугольника Паскаля. В этом случае ап = 1, если в (23) множитель с = с д = Л”.
Доказательство. По лемме 1 последние строки обратных матриц Тейлора Т* (21) размерности г + 1 есть (г + 1)-векторы коэффициентов разделенных разностей порядка г и эти векторы не зависят от точки т = Ыо разложения (21) выборок отсчетов V. Отсюда следует, во-первых, что метод конечных разностей означает использование в лемме 2 вместо матрицы Т-1 — обратной матрицы Тейлора вида (21) — нижнетреугольной матрицы Т- . Во-вторых, г-я строка ее треугольника есть вектор [Дг]*/К г-й разделенной разности. Матрица конечных разностей Т-1, следовательно, имеет вид:
(28)
Т-1 = {| [Af/h , 0n_i |}0 = DP, D = diag{h-i}П, P = P-{Pij}П, где
P_ij = 0 при j > i и P_ij = (—1)i-j при i ^ j.
Если an = 1, то формула (23) примет вид а* = а* ■ T—1 ■ с д, где с д = hn. □
Нижнетреугольные матрицы вида P можно назвать матрицами Паскаля. Они имеют интересную особенность, используемую в доказательстве теоремы 3. Эта теорема отвечает на вопрос: какому представлению (виду интерполяции) отсчетов выборок v через производные (вместо тейлоровского v = Tw из (21)) соответствует их выражение v = Тд-w, вытекающее из способа дискретизации в лемме 4? Каков вид нижнетреугольной матрицы Тд?
Теорема 3. Дискретизация ДУ разделенными разностями по формуле а* = а*Т_1hn эквивалентна, вместо их тейлоровского разложения (21), следующему представлению отсчетов выборок % с помощью нижнетреугольной матрицы Тд и чисел треугольника Паскаля:
{i Л п
¿0+ hj/2)j , = (i -ij.),j,. (29)
Доказательство. При тейлоровском приближении Wk = Т-1ik производных на основе разложения (21), все n + 1 отсчетов выборки % интерполируются одним полиномом порядка n, и именно его производные в заданной точке т = hio принимаются в качестве оценок производных аппроксимируемой им функции. Ими заменяются производные в ДУ (лемма 4).
Разностное приближение производных Wk = Т_1 % означает, что для оценки j-й производной используется «свой> полином j-го порядка. Как было показано в следствии леммы 1 (на примере полинома порядка n — n-полинома), j-я производная j-полинома, постоянная на минимальном интервале Р/(vk), содержащего эту (j + 1)-выборку, есть оценка j-й производной аппроксимируемой им функции y(t) в центре этого интервала, то есть в точке hj/2. Это и указано в формуле (29) теоремы 3.
Определим оператор dT = sTd = (s — 1)/h, 0 ^ т ^ h разностного дифференцирования функции y(t) в точке tk + т: dTyk = sTd ■ yk = y;(tk + т). Здесь d = d/dt —оператор дифференцирования, а sT — символ сдвига sTy(tk) = y(tk + т) на величину т ^ h. Он указывает какой точке tk + т на интервале длины h значению производной y;(tk + т) приписывается
значение разделенной разности ((в — 1)/Н)у^ . Наконец, в в операторе разности (в — 1) — «:обычный> оператор сдвига из леммы 3: ву^ = у^+ь Обратим оператор разностного дифференцирования йт = (в — 1)/Н, то есть выразим оператор сдвига в через оператор йт. Получим: в = йтН +1 и
вУк = Ук+1 = (1 + Н^т )ук = (1 + Нвт ^)уй = ук + Ну'(¿й + т).
Степени (йт)г = ((в — 1)/Н)г оператора йт (их коэффициенты) определяют строки матрицы Т-1 (28) и преобразование w = Тд1^, используемое для разностной дискретизации в лемме 4. Обратное нижнетреугольное преобразование (29) V = Т^, представлено в теореме 3. Его строки есть коэффициенты степеней вг = (йтН + 1)г оператора в. Представление (29) есть альтернатива тейлоровскому представлению (21). Докажем теперь формулу (29).
Можно показать, что Р+Р_ = I, то есть Р—1 = Р+. Здесь Р+ — это матрица, аналогичная матрице Р_ из (28), но со знакопостоянными (положительными) числами треугольника Паскаля (без множителя ( — 1)г_ в компонентах Р^).
В соответствии с леммой 4 матрица Тд = Р_1^_1 заменяет тейлоровскую матрицу Т в системе уравнений (21) для отсчетов: V = Tдw. Это означает использование для интерполяции отсчетов выборки % из РУ (5) вместо интерполяционного полинома (21) следующего формального разложения
% = = {(вт й + 1)^ }”=0, где т = Н/2.
Отсюда получаем результат (29) — альтернативу разложения Тейлора (21): і-й отсчет выборки V есть сумма производных у(^(’’Н/2), j = 0, і с весами бинома і-й степени. □
3. Линейное приближение равномерной дискретизации
3.1. Основной результат
В этом разделе будет доказан следующий результат.
Теорема 4. Линейным по а приближением к формуле (14) аналитической дискретизации
а* = -^Р-Щ = -^*ФгаР-1, где Р = Р(А) = {^Ф*}^"1 = {^* ехр(АЛг)}п_1, (30)
является формула локальной тейлоровской дискретизации (23) леммы 2.
Докажем вспомогательные и представляющие самостоятельный интерес утверждения.
Лемма 5. Результат аналитической дискретизации по формулам (30) инвариантен к преобразованию фазовых координат в исходной системе (9).
Доказательство. Пусть Ф = ехр(АуЛ) и Ф = ехр(АЛ) есть матрицы перехода дискретной системы (10) в случаях, когда исходная континуальная система (9) представлена в каноническом (8) и общем виде (9) соответственно. Матрицы наблюдаемости этих систем на сетке есть соответственно Р = Р(а) = Р(Ау) = Р(Ау(а)) = {е0/га фг}П и Р из (30). Заменой переменных ад = Рдх, где Рд = {^*Аг}Щ 1 в (9) придем к равенству Ау = РдАРД-1. Отсюда выводим: Ф = РД 1фРд ^ Ф* = Рд1фгРд. Из замены переменных в (9) следует также, что е0/п = ¿*Рд1 = й*|й, В| 1*, где В* = {^*Аг}Щ 1 и, кроме того,
e0/ra
F = Ігі*Р-1Ф^ Рд ^ F-1 = Р-1і?_1. Используем эти результаты в (30):
а* = -^*Рд1фгаР-1 = -e^F-1 где F = Р(а) = {e0/n exp (Af(a)hi) j^-1- (31)
Получили требуемое: результат дискретизации динамических систем, связанных преобразованием фазовых координат, одинаков. В частности, он одинаков для общего (9) и канонического (8) описаний. Этот результат (31) и необходим для доказательства теоремы 4. □
3.2. Полиномиальная система и нулевое приближение
Лемма 6. Первая строка т0 1) + (п) обратной матрицы Тейлора Т 1(п) на (п — 1)-выборке с точкой разложения п, следующей за правой граничной точкой /га_1^) интервала интерполяции выборки V (п — 1)-полиномом, есть — [Дп] — прогнозирующий вектор п-й конечной
разности с обратным знаком.
Доказательство. Формула ад = Т 1(п) ■ V выражает п-вектор ад оценок производных (от
0-й до (п — 1)-й) через п-вектор отсчетов — п-выборку V. Поэтому, первая строка т0 ^*(п)
___1 ____________________________________________________________________^
матрицы Т (п) должна быть — со знаком «минуо — прогнозирущим вектором [Дп]
коэффициентов разности Дп для прогноза значения нулевой производной (то есть отсчета
самой функции) на точку п. Эта точка находится вне интервала /га_1(^), но в интервале
1п^).
РУ вида (5) для (п — 1)-полинома есть [Дга]* ■ V = 0. Прогнозирующий вектор его коэффициентов есть [Дп] , а рекуррентная формула прогноза (п — 1)-полинома, подобная
----^
уравнению (6) для РУ общего вида, есть у^+1 = — [Дп] ^. Коэффициенты правой части этого уравнения и есть компоненты первой строки матрицы Т (п). □
___і __________^ _______і
Следствие 4. е0/пТ (п) = — [Дп] = —т*Т .
Доказательство. В результате леммы 6 нужно учесть первое равенство леммы 3. □
Полиномиальной системой мы называем систему п уравнений (8), если а = 0. Тогда в (8): Ау(а) = Ау(0) = /П — матрица сдвига. Эта система эквивалентна уравнению у(п) (¿) = 0, степенного полинома порядка п — 1 в заданном интервале /у.
Лемма 7. Пусть дискретная система (10) описывает отсчеты полинома (п — 1)-й степени, то есть является реализацией полиномиальной (при а = 0, то есть Ау = Ау(0)) системы (8) на сетке /^. Матрица наблюдаемости (12) такой системы (8) на Н-сетке —
ф = |e*/nфг j фk, где фг = фг(0) = exp (Ay(0)hi) — есть матрица Тейлора T(—k).
Доказательство. Вычислим матрицу наблюдаемости F = F(0) = Fo (12) для этого случая:
F = ie0/n Ф lo-1 = ie0/n eXP (JКак известно [2L exp (/lh) = En-1(h)j(11)У.Л-
Следовательно, ф(i+k)(0) — верхнетреугольная теплицева матрица с числами (h(i + k))j/j! на j-й наддиагонали (на главной — j = 0). Ее первая строка e0/o Фг(0) = e0/o exp (/ 1h^ — есть
|(hi)YjllO^o1- Это не что иное как i-я строка матрицы Тейлора T(—k) (21). Из выписанной
выше формулы для матрицы наблюдаемости F следует, что ее i-я строка и есть i-я строка этой матрицы Тейлора. Лемма доказана. □
Следствие 5. Нулевое по а приближение F = F(A(0)) = { exp (A(a)h(i+fc)) Щ^о1 матрицы наблюдаемости F есть матрица Тейлора T(—k).
Доказательство. Для получения нулевого приближения надо положить а = 0 в формуле (31) для матрицы наблюдаемости F системы (8) на сетке □
Теорема 5. Нулевое по a приближение к аналитической дискретизации (30) есть оэффи-циенты уравнения для полинома (n — 1)-го порядка: а* = — [Дп] .
Доказательство. Учтем, во-первых, лемму 5. Во-вторых, учтем в формуле (31) следствие 5, а именно, что F—n(0) = T (n). Результат леммы 6 и формула (31) для а* приводят к доказываемому утверждению. □
3.3. Линейное приближение
Доказательство теоремы 4. Для доказательства необходимо убедиться, что линейные приближения по a компонент формулы дискретизации (30) дают результат леммы 2 и теоремы 2. С учетом результата леммы 5 достаточно провести доказательство для дискретизации канонического описания (8). Будем исходить из первой формулы в (31): а* = —ео/п 1П ' F1-1 (Af(a)). Она доказывается несколько сложнее, но соответствует более общему случаю, когда угловая компонента $ = 0 в представлении (25) матрицы T.
Получим сначала разложение Тейлора по степеням h до степени n включительно векторной функции — строки /* = /n(â) = e0/n Фп = e*/n exp (Af(â)hn), где Af(a) — матрица фробениуса из (8). Воспользуемся результатом леммы 7, где это разложение было получено до степени n — 1. Эти n членов ряда Тейлора функции e*/n Фn от вектора a не зависят, поскольку, как можно увидеть из формулы для матрицы Af в (8), e0/nAF
= eo/) ,('1)' = e*/n,
если j < n. Если j = n, то e0/nAF = —a*.
Используем эти факты, а также обозначения в (25) представления матрицы Тейлора T = Tn в окаймленном виде. Выпишем n + 1 членов ряда Тейлора функции f*(a). Последний член содержит линейную часть по a:
n n—1
fn = e0/nфп = Y e0/n(Ahn)ÿj! = Y e0/n(hn)j/j! — $a* = T*— $a*, $ = (hn)7n! • (32)
oo
Теперь выписываем аналогичное (до линейного по a члена) разложение матричной функции F = F (a) = {/*}n— 1. Используем формулы (21), (32) и результаты леммы 7:
n— 1
Р = ШП-1 = е**/га(Ьг )Ш - * ■ а* = Т - * ■ а*, где Т = Т(0).
I 0 J г=0
Для обращения матрицы Р применим известную формулу для обращения матриц с аддитивными приращениями [5]. Используем обозначения из представления (26):
і — — 1 — — 1- f ——1-\ — 1 і.— — 1 — — 1 — — 1- і і.— — 1
— 1 ___ ГТ1 1 І ГТ1 1 1 І-j + \ 1 __ 'Т'1 1 I 'Г'1 + >- — 1 1
р-1 = т " + т ^ 1 -ат ^ 1 ат "= т "+ т ч■ е_1■ ат
Подставляем результат (32) и последнее выражение для Р-1 в первую формулу в (31): а* = - (т* - $а*) (т-1 + Т-1£ ■ е-1 ■ а*Т-1) =
а*Т “V + а*Т-1£ ■ е-1 ■ а*Т-1$ - т*Т-1 - т*Т-1£ ■ е-1 ■ а*Т-1 =
—*т^—1^ ,-1—*т^-Ча I гдп* ^*гг-^ ,-1 а*Т-1 =
(е + â*T 1^Є—1â*T V + [Дга]* - т*T 1t ■ є
- т*Т — 1^ е—1â*T—1 + [ДП]*.
Это совпадает с результатом лемм 2, 3 и теоремы 2, с учетом формул (23) и (26). Теорема 4 доказана. □
Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (проект № 10-01-00035) и Сибирского отделения Российской академии наук (междисциплинарный проект № 80).
Литература
1. Егоршин, А.О. Идентификация стационарных моделей в унитарном пространстве / А.О. Егоршин // Автоматика и телемеханика. - 2004. - Т. 65(12). - С. 29-48.
2. Гантмахер, Ф.Р. Теория матриц / Ф.Р. Гантмахер. - М.: Наука, 1966.
3. Егоршин, А.О. Об отслеживании параметров экстремума в вариационной задаче идентификации / А.О. Егоршин // Вестник НГУ. Серия: Математика, механика, информатика. - 2011. - Т. 11, вып. 3. - С. 95-114.
4. Егоршин, А.О. Об одном способе оценки коэффициентов моделирующих уравнений для последовательностей / А.О. Егоршин // Сиб. журн. индустр. матем. - 2000. - Т. 3, № 2.
- С. 78-96.
5. Эльясберг, П.Е. Определение движения по результатам измерений / П.Е. Эльясберг. -М.: Либроком, 2011.
Алексей Олегович Егоршин, кандидат физико-математических наук, старший научный сотрудник, Институт математики им. С.Л. Соболева СО РАН, egorshin@math.nsc.ru.
MSC 65F25; 15A03
On Linear Differential Equation Discretization
A.O. Egorshin, Sobolev Institute of Mathematics, Siberian Branch of the Russian Academy of Sciences (Novosibirsk, Russian Federation)
Some problems of obtaining the discrete description of the first order differential system (DS) on the uniform lattice have been considered. These DS are regarded in the form of system n of the first order ordinary linear differential equations with constant coefficients or as one n-order equation for the observed functional of the DS state. The problems under consideration are of some importance for the problems of the variational identification and approximation of the dynamic processes by means of that type models on the finite interval. There are compared the analytic uniform method of discretization (based on Cayley
- Hamilton theorem) and that of the local one on the basis of the interpolation of the samples of n +1 counting by Taylor polynomials to the power n. There have been obtained the general formula of the local discretization that makes it possible to compare its difference and interpolarization methods. It has been shown by using Vandermond inverse matrices that in the obtained general formula of the local discretization n +1 Taylor matrices (from Taylor polynomial coefficients) correspond to its interpolational method while n+1 Pascal matrices (from Pascal triangle numbers) correspond to the difference method.
It has been shown that matrix nondegeneracy of the DS observability on the lattice is a necessary and sufficient condition both for analytic discretizability and for reducing the discete system (of the DS description of the lattice) to Frobenius canonical form. It is equivalent to one ordinary difference equation for the observed variable with constant coefficients. This equation is a basis of the well-known variational method of identification. It has been shown that interpolation method of the local discretization is the first order linear approximation of the uniform analytic discretization formula. It has been demonstrated that its zero order approximation does not depend on the DS coefficients and is a vector of the coefficients of the n-th difference. We conclude that zero order approximation of the observability matrix of DS and of the observability matrix of the polynomial system y(n) =0 on the lattice is Taylor n-matrix.
Keywords: variational approximation and identification, discretization of differential equation, analytical discretization, linear approximation, Cayley - Hamilton theorem, local discretization, Teylor polynomial, Vandermond matrices, Pascal triangle.
References
1. Egorshin A.O. Parameters Optimisation of the Stationary Models in Unitary Space, Automatica and Remote Control, 2004, vol. 65, no. 12, pp. 1734-1756.
2. Gantmacher F.R. The Theory of Matrices, New York, Chelsea Publishing Company, N.Y. 1959.
3. Egorshin, A.O. On Extremum Parameters Tracing in Variational Identification Problem, [Ob otslezhivanii parametrov ekstremuma v variatsionnoy zadache identifikatsii]. Vestnik NGU, Seriya: Matematika, Mechanika, Informatika, 2011, vol. 11, no. 3, pp. 95-114.
4. Egorshin A.O. On One Estimation Method of Modelling Equation Coefficients for Sequences [Ob Odnom sposobe otsenki koeffitsientov modeliruyushchikh uravneniy dlya posledovatel‘nostey]. Sibirskiy Zhurnal Industrialnoy Matematiki [Sib. Zh. Ind. Math.], 2000, vol. III, no. 2, pp. 78-96.
5. El‘yasberg P.E. Determination of Moving via Result of the Measurement [Opredelenie dvizheniya po rezul‘tatam izmereniy]. Moscow, URSS, 2011.
Поступила в редакцию 17 июля 2012 г.
