close

Вход

Забыли?

вход по аккаунту

?

Mironovsky

код для вставкиСкачать
ФЕДЕРАЛЬНОЕ АГЕНТСТВО ПО ОБРАЗОВАНИЮ
Государственное образовательное учреждение
высшего профессионального образования
САНКТ!ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
АЭРОКОСМИЧЕСКОГО ПРИБОРОСТРОЕНИЯ
Л. А. Мироновский
МОДЕЛИРОВАНИЕ
РАЗНОСТНЫХ
УРАВНЕНИЙ
Учебное пособие
Рекомендовано УМО по образованию в области радиотехники,
электроники, биомедицинской техники и автоматизации
в качестве учебного пособия для студентов высших учебных
заведений, обучающихся по специальности 210100 –
"Управление и информатика в технических системах"
Санкт!Петербург
2004
УДК 681.3.06:51(075)
ББК 22.18я7
М63
Мироновский Л. А.
М63 Моделирование разностных уравнений: Учеб. пособие/
СПбГУАП. СПб., 2004. 108 с.: ил. ISBN 5!8088!0135!4
Изложены основы теории разностных уравнений, методы их анали!
тического решения и компьютерного моделирования. Описаны основ!
ные свойства линейных разностных уравнений и применение жорда!
новых цепочек векторов для их решения. Приводятся примеры моде!
лирования разностных уравнений с начальными и краевыми условия!
ми в пакете SIMULINK.
Учебное пособие предназначено студентам, обучающимся по спе!
циальности 210100 «Управление и информатика в технических систе!
мах», а также родственным специальностям и направлениям, таким
как 220100 «Информатика и вычислительная техника», 552800 «Вы!
числительные сети и системы» и др.
Оно будет также полезно магистрантам и аспирантам, которым при!
ходится сталкиваться с необходимостью аналитического и компьютер!
ного анализа дискретных систем и их математических моделей.
Рецензенты:
доктор технических наук профессор В. А. Слаев;
кандидат технических наук профессор Г. С. Бритов
Утверждено
редакционно!издательским советом университета
в качестве учебного пособия
ISBN 5!8088!0135!4
2
©
ГОУ ВПО “Санкт!Петербургский
государственный университет
аэрокосмического приборостроения”,
2004
ПРЕДИСЛОВИЕ
Среди различных математических моделей, применяемых для опи!
сания реальных систем, важное место занимают разностные уравнения,
причем их роль постоянно возрастает. Они широко используются в на!
уке и технике при описании самых разных процессов и систем – элект!
рических, механических, биологических, демографических, экономи!
ческих и др. В качестве примеров можно назвать анализ цепных (лест!
ничных) схем в теории цепей, модели длинных линий в электротехни!
ке, методы численного интегрирования в вычислительной математике,
методы сеток и конечных элементов в математической физике. К разно!
стным уравнениям приводят многие экологические задачи и модели
популяционной динамики, экономические задачи (расчет сложных про!
центов, управление банковскими депозитами и т. п.), демографические
модели (прогнозирование половозрастной структуры населения, пост!
роение демографических пирамид).
В то же время в учебных планах технических вузов разностным
уравнениям уделяется недостаточное внимание. Они не рассматри!
ваются в стандартных курсах высшей математики, и лишь в курсах
по теории автоматического управления студентов кратко знакомят с
применением z!преобразования для исследования импульсных сис!
тем. Учебная литература по разностным уравнениям также весьма
скудна. Настоящее учебное пособие призвано хотя бы немного ис!
править существующее положение. Оно написано на базе лекций по
курсу «Специальные разделы прикладной математики», читавше!
муся автором на протяжении ряда лет в ГУАП.
В пособии изложены основы теории разностных уравнений, мето!
ды их аналитического решения и компьютерного моделирования. Все
разделы имеют примерно одинаковый объем и снабжены большим
количеством примеров.
В первом разделе описаны основные свойства линейных разно!
стных уравнений и методы их решения.
Второй раздел посвящен применению жордановых цепочек векто!
ров для решения дифференциальных и разностных уравнений.
Третий раздел содержит материал о z!преобразовании и дискрет!
ных передаточных функциях.
3
В четвертом разделе приводятся примеры моделирования разно!
стных уравнений с начальными и краевыми условиями в пакете
SIMULINK.
Последний, пятый, раздел посвящен нелинейным разностным урав!
нениям.
4
1. ОСНОВНЫЕ СВОЙСТВА РАЗНОСТНЫХ УРАВНЕНИЙ
1.1. Примеры разностных уравнений
Многие научные и технические задачи приводят к необходимости
решать разностные уравнения. Примерами могут служить проекти!
рование импульсных систем в теории автоматического управления,
расчет цифровых фильтров в теории связи, анализ погрешностей ре!
шения дифференциальных уравнений на ЭВМ и др.
Разностные уравнения (другие названия: уравнения в конечных
разностях; возвратные последовательности) по своим свойствам и
области применения очень близки к дифференциальным уравнени!
ям. Отличие состоит в том, что дифференциальные уравнения связы!
вают значение функции и производных от нее в один и тот же момент
времени
f (x(n) (t), ..., x1 (t), x(t)) 1 0,
а разностные уравнения – значения функции в различные моменты
времени
f (x(t 1 n), ..., x(t 1 1), x(t)) 2 0.
(1.1)
Перечислим некоторые источники разностных уравнений: диск!
ретизация обыкновенных дифференциальных уравнений и уравне!
ний в частных производных; модели объектов с дискретным време!
нем (задача Фибоначчи о размножении кроликов); объекты с диск!
ретным пространством (R!2R!цепь); анализ математических рядов и
рекуррентных соотношений.
Пример 1. Непрерывную функцию x(t) 1 e 12t можно описать с по!
мощью дифференциального уравнения
x1 (t) 1 2x(t) 3 0, x(0) 3 1.
Если рассматривать дискретные моменты времени t 1 0, 1, 2... , то
эту же функцию можно описать с помощью разностного уравнения
x(t 1 1) 1 2x(t) 3 0, x(0) 3 1, 2 3 4e 12 .
(1.2)
5
Непосредственной подстановкой легко убедиться, что x 1 e 12t удов!
летворяет обоим уравнениям.
Так же как и дифференциальные, разностные уравнения делятся
на линейные и нелинейные, на однородные и неоднородные, на урав!
нения с постоянными и переменными коэффициентами. Например,
уравнение (1.2) является линейным однородным разностным урав!
нением первого порядка.
Линейное однородное разностное уравнение n!го порядка с посто!
янными коэффициентами имеет вид:
x(t 1 n) 1 an 11x(t 1 n 2 1) 1 ... 1 a1x(t 1 1) 1 a0 x(t) 3 0.
(1.3)
Решением этого уравнения называется решетчатая функция x(t),
t = 0, 1, 2, ..., обращающая уравнение в тождество. Можно говорить
и о непрерывном решении – функции x(t), удовлетворяющей уравне!
нию (1.3) при любом 0 1 t 1 2 .
Иногда, чтобы подчеркнуть дискретный характер изменения вре!
мени, уравнение (1.3) записывают в форме
xn 1 k 1 an 21xn 1 k21 1 ... 1 a1xk11 1 a0xk 2 0.
Рассмотрим несколько простых задач, приводящих к линейным
разностным уравнениям.
Пример 2. Пусть имеется геометрическая прогрессия
1, a, a2, a3, .... Ее можно описать однородным разностным уравнени!
ем первого порядка
x(t 1 1) 2 ax(t), x(0) 2 1.
(1.4)
Действительно, подставляя последовательно t 1 0,1, 2, 3, 1 , по!
лучаем
x(0) 1 1, x(1) 1 a, x(2) 1 a2,
x(3) 1 a3, ... .
Таким образом, геометрическая прогрессия является решением
простейшего разностного уравнения, подобно тому как экспонента
является решением простейшего дифференциального уравнения. Эта
аналогия распространяется и на уравнения более высоких порядков.
Как правило, общее решение уравнения (1.3) представляет собой
линейную комбинацию n геометрических прогрессий (для дифферен!
циальных уравнений мы имели линейную комбинацию экспонент).
Пример 3. Пусть имеется арифметическая прогрессия
a, a 1 d, a 1 2d, ... . Ее можно описать неоднородным разностным урав!
нением первого порядка
x(t 1 1) 2 x(t) 1 d, x(0) 2 a
6
(1.5)
или однородным уравнением второго порядка
x(t 1 2) 2 2x(t 1 1) 1 x(t) 3 0 .
(1.6)
(Убедитесь, что это действительно так).
Пример 4. Предположим, что необходимо вычислить напряже!
ние un на выходе электрической цепи (рис. 1.1), зная входное напря!
жение u0 и полагая R = r . Такие цепи встречаются, например, в циф!
роаналоговых преобразователях.
Попытка непосредственно решить эту задачу с помощью законов
Кирхгофа и Ома приводит к громоздким вычислениям. Значительно
более коротким оказывается путь решения, использующий состав!
ление разностного уравнения:
R
u0
R
u1
R
u2
R
un–1
un
r
r
r
r
Рис. 1.1
Запишем уравнение токов для первого узла цепи:
u0 1 u1 u1 u1 1 u2
2 3
,
R
r
R
R
u2 3 16 2 4 27 u1 4 u0 5 0.
r 9
8
Аналогично, для k + 1!го узла получаем
R
uk12 3 16 2 4 27 uk 11 4 uk 5 0,
r9
8
(1.7)
или, учитывая, что R = r:
uk12 1 3uk 11 2 uk 3 0, k 3 0, n 1 2.
Уравнение (1.7) представляет собой разностное уравнение второ!
го порядка, поскольку связывает значение функции u в точке k + 2
с ее значениями в двух предыдущих точках k + 1 и k.
Пример 5. Рассмотрим старинную задачу о размножении кроли!
ков (задача Фибоначчи или Леонардо Пизанского, итальянского ма!
тематика, жившего около 1200 г.). В этой задаче требуется опре!
делить число пар зрелых кроликов, образовавшихся от одной пары
7
в течение года, если каждая пара кроликов рождает ежемесячно но!
вую пару и новорожденные достигают зрелости через месяц. Таким
образом, получается некоторая последовательность; нас интересует,
каким уравнением она описывается.
В первый момент времени число пар равно x0 = 1. Через месяц
прибавится пара новорожденных, но число зрелых пар не изменится
x1 = 1.
Через два месяца крольчата достигнут зрелости и общее число пар
будет x2 = 2.
Через k месяцев число зрелых пар будет xk, а через k + 1 месяцев
xk+1 , но так как к этому времени появится xk пар приплода, то через
k + 2 месяцев общее число зрелых пар
xk 12 1 xk 11 2 xk .
(1.8)
Мы получили разностное уравнение второго порядка. Оно описы!
вает последовательность, где каждый следующий член равен сумме
двух предыдущих:
1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ... .
(1.9)
Числа, входящие в эту последовательность, называются числами
Фибоначчи и обладают многими интересными свойствами. Напри!
мер, в ботанике известно так называемое явление филлотаксиса,
связанное со спиральным расположением листьев на стеблях, се!
мян в плодах. Спиральные линии образуют чешуйки еловых и сосно!
вых шишек, семена в головке подсолнечника и т. д. Оказывается,
что число “правых” и “левых” спиралей при этом представляет собой
соседние числа Фибоначчи. Так, для сосновых шишек число “пра!
вых” и “левых” спиралей равно 5 и 8, а для еловых – 8 и 13.
В двух последних примерах мы свели исходные задачи к решению
разностных уравнений. Перейдем к описанию методов решения та!
ких уравнений.
1.2. Решение однородных разностных уравнений
Свойства разностных и дифференциальных уравнений и методы
их решения во многом совпадают. Так, разностное уравнение n!го
порядка (1.3) имеет ровно n линейно независимых частных решений
x1(t), ..., xn(t). Любая их линейная комбинация, например x1(t) +
+ 3x2(t), также является решением.
Общее решение имеет вид
x(t) 1 c1x1 (t) 2 ... 2 cnxn (t),
где c1, ..., cn – произвольные коэффициенты.
8
(1.10)
Известны два основных метода решения линейных разностных урав!
нений – с помощью характеристического полинома и с использованием
z!преобразования, аналогичного преобразованию Лапласа.
Рассмотрим первый из них. Будем искать решение в виде x(t) 1 zt ,
где z – некоторое число. Если t – дискретно, то x(t) – геометрическая
прогрессия, если же t – непрерывно, то x(t) – показательная функ!
ция. В обоих случаях ее можно записать в виде eat, где a = ln z.
Подставляя x(t) в (1.3) и сокращая на zt, получаем алгебраи!
ческое уравнение
zn 1 an 11zn 11 1 ... 1 a0 2 0,
(1.11)
которое называется характеристическим (сравните с аналогичной
процедурой получения характеристического полинома для дифферен!
циального уравнения).
Вид решения разностных уравнений зависит от типа корней ха!
рактеристического полинома, которые могут быть вещественными,
комплексными, простыми и кратными.
Вещественные корни
Если уравнение (1.11) имеет вещественные и различные корни
z1, ..., zn, то получаем n линейно независимых частных решений урав!
нения (1.3):
x1(t) 1 z1t, ..., xn (t) 1 znt .
Общее решение в соответствии с (1.10) имеет вид:
x(t) 1 c1z1t 2 ... 2 cnznt .
(1.12)
(Сравните с аналогичной формулой решения дифференциальных
уравнений).
Рассмотрим несколько примеров.
Пример 1. Решим разностное уравнение
x(t + 2) – 5x(t + 1) + 6x(t) = 0.
Характеристическое уравнение имеет вид z2 – 5z + 6 = 0.
Его корни вещественные и различные: z1 = 3, z2 = 2.
Общее решение: x(t) = c1 3t + c2 2t.
Пример 2. Найдем общее решение уравнения электрической
цепи (1.7).
Составляем характеристическое уравнение:
z2 – 3z + 1 = 0.
9
Оно имеет вещественные и различные корни
z1,2 1 0,5(3 2 5); z1 3 2,62; z2 3 0,38 .
Следовательно, общее решение имеет вид
(1.13)
uk = c1 2,62k + c2 0,38k.
Пример 3. Найдем общее решение уравнения Фибоначчи (1.8).
Составляем характеристическое уравнение
z2 – z – 1 = 0.
Корни этого уравнения также вещественны и различны:
z1 1 0,5(1 2 5) 3 1,62 z2 1 0,5(1 4 5) 3 40,62.
Общее решение уравнения Фибоначчи имеет вид
xk 1 c1 1,62k 2 c2 (30,62)k.
(1.14)
Пример 4. Решим разностное уравнение
x(t + 2) – 3x (t + 1) + 2x(t) = 0.
Составляем характеристическое уравнение: z2 – 3z + 2 = 0.
Находим его корни: z1 = 1, z2 = 2.
Общее решение: x = c1 + c2 2t.
Комплексные корни
Если характеристическое уравнение имеет пару комплексно!со!
пряженных корней z1,2 1 2 3 i4 , то соответствующее им реше!
ние x(t) 1 c1z1t 2 c2z2t может быть записано в вещественном виде
x(t) 1 2t (c1 sin 3t 4 c2 cos 3t),
1
(сравните с аналогичной формулой для
4
дифференциальных уравнений).
Из приведенной формулы вытекает, что если корни характерис!
тического уравнения лежат внутри единичного круга, то решение
разностного уравнения будет устойчивым. Это широко известный
критерий устойчивости дискретных систем. Напомним, что для ус!
тойчивости непрерывных систем корни характеристического уравне!
ния должны лежать в левой полуплоскости.
Пример 5. Решим разностное уравнение
x(t + 2) + 2x(t + 1) + 4x(t) = 0.
2
2
где 2 3 4 5 1 , 6 3 arctg
Характеристическое уравнение: z2 + 2z + 4 = 0.
10
Его корни: z1,2 1 21 3 3 i . Модуль и аргумент корней:
2 3 2, 4 3
51
.
3
Общее решение:
x(t) 2 2t (c1 cos
51
51
t 3 c2 sin t).
3
3
Кратные корни
Если один из корней характеристического уравнения zi имеет крат!
ность k, то соответствующее ему слагаемое в общем решении умно!
жается на полином степени k – 1
xi (t) 1 (c1 2 c2t 2 ... 2 cktk11)zit
(это полностью аналогично случаю дифференциальных уравнений).
Пример 6. Решим разностное уравнение с кратными веществен!
ными корнями
x(t + 3) + 7x(t + 2) + 15x(t + 1) + 9x(t) = 0.
Характеристическое уравнение: z3 + 7z2 + 15z + 9 = 0.
Его корни: z1 = – 1, z2 = –3, z3 = –3.
Общее решение:
x(t) = c1 (–1)t + (c2+c3 t) (–3)t.
Пример 7. Решим разностное уравнение с кратными комплексны!
ми корнями
x(t + 4) + 2x(t + 2) + x(t) = 0.
Характеристическое уравнение: z4 + 2z2 + 1 = 0.
Его корни: z1 = z2 = i; : z3 = z4 = –i.
Общее решение:
1
1
x(t) 2 (c1 3 c2t)cos t 3 (c3 3 c4t)sin t.
2
2
Учет начальных и краевых условий
Общее решение, в которое входят произвольные постоянные ci,
задает семейство решений разностного уравнения. Для выделения из
этого семейства нужного нам частного решения, необходимо найти
эти постоянные. В случае дифференциальных уравнений их нахо!
дят, задавая начальные или краевые условия, общее число которых
11
равно порядку уравнения. Аналогично для разностных уравнений
для нахождения постоянных c1, ..., cn надо задать n условий, связы!
вающих значения x(t) в различных точках.
Чаще всего встречаются два типа условий:
– когда задается начальная последовательность x(1), x(2), ... , x(n);
– когда часть значений x(t) задается в начале, а часть – в конце
интервала решения.
Условия второго типа называются к р а е в ы м и .
Пример 8. Требуется решить разностное уравнение
xn 12 1 3xn 11 2 2xn 3 0
для трех вариантов начальной последовательности: а) x0 = 1, x1 = 3;
б) x1 = 1, x2 = 3; в) x1 = 2, x2 = 0.
Р е ш е н и е . Составляем характеристическое уравнение и нахо!
дим его корни:
z2 1 3z 2 2 3 0,
z1 3 1, z2 3 2.
Общее решение имеет вид xn 1 c1 2 c2 3 2n.
Постоянные с1 и с2 находим из следующих начальных условий:
а) с1 + с2 = 1, с1 + 2с2 = 3, откуда с1 = – 1, с2 = 2, т. е. xn = 2n+1 – 1;
б) с1 + 2с2 = 1, с1 + 4с2 = 3, откуда с1 = – 1, с2 = 1, т. е. xn = 2n – 1;
в) с1 + 2с2 = 2, с1 + 4с2 = 0, откуда с1 = 4, с2 = – 1, т. е. xn = 4 – 2n.
Пример 9. Найти решение уравнения из предыдущего примера,
если заданы краевые условия x1 = 2, x10 = 1024.
Полагая в формуле для общего решения n = 1 и n = 10, получаем
систему уравнений для определения c1 и c2:
c1 + 2c2 = 2, c1 + 210 c2 = 1024,
откуда с1 = 0, с2 = 1.
Следовательно, искомое решение имеет вид: x(t) = 2n .
Пример 10. Выделить из общего решения (1.14) частное решение,
соответствующее последовательности Фибоначчи (1.9).
Общее решение: xk = c1 z1k + c2 z2k . Начальные условия: x0 = 1, x1 = 1.
Cистема уравнений для определения с1 и с2:
3c1 1 c2 2 1,
4
5c1z1 1 c2z2 2 1,
откуда
c1 2
1 1 z2
3 0,72,
z1 1 z2
c2 2 1 1 c1 3 0,28.
12
Решение задачи Фибоначчи (общий член последовательности
Фибоначчи):
xk 1 0,72 2 1,62k 3 0,28 2 (40,62)k.
Не выполняя округлений, то же решение можно записать в виде
k 11
k 11
1
316 5 4 2
1 73 1 5 5 4
xk 9
6
8.
5 7 2 2 8
Это известная формула Бине. Хотя в нее входят радикалы, она
дает при любом k целые числа. В частности, при k = 0 и k = 1 получаем
x0 = x1 = 1, а при k = 12 (решение задачи о кроликах) имеем x12 = 233.
Пример 11. Требуется вычислить определитель трехдиагональ!
ной n ´ n матрицы A вида
11 1 0 2
1 4
3
An 5 31 1 1 4 .
1
30 1 1 4
6
7
Для n = 1, 2, 3 находим
A1 1 1;
A2 1
1 1
1 0;
1 1
1 1 0
A3 1 1 1 1 1 21.
0 1 1
Чтобы получить общую формулу, раскроем определитель матри!
цы An+2 по элементам первой строки:
An 12 1 An 11 2 An ;
1
1
например, A4 1
0
0
1
1
1
0
0
1
1
1
0
1 0 0
0
1 A3 2 0 1 1 1 A3 2 A2 .
1
0 1 1
1
Обозначая xn 1 An , приходим к разностному уравнению
xn 12 1 xn 21 2 xn 3 0 .
Выписываем характеристическое уравнение z2 – z + 1 = 0 и нахо!
дим его корни
1
41
z1,2 1 2 3 i 3 /2, 2 3 1, tg4 3 5 3; 4 3 1200 3 21 и 2 3 2400 3 .
3
2
3
13
21
21
n 3 c2 sin n.
3
3
Для определения произвольных коэффициентов используем на!
чальные условия:
Общее решение имеет вид: xn 2 c1 cos
3
11
4 c2 5
1
,
c 1
2
25
76 1
3
3
1
4 c2 5
n 2 2 : 0 2 c1 cos240 3 c2 sin240 2 c1
c2 1 1.
2
2 58
n 2 1:
1 2 c1 cos120 3 c2 sin120 2 c1
1
21
21
cos n 3 sin n .
3
3
3
Результат выглядит неожиданно и вряд ли может быть «угадан»
путем рассмотрения значений определителей для n = 1, 2, 3, ... .
Пример 12. Вычислить определитель трехдиагональной n ´ n мат!
рицы вида
О т в е т . An 2
14 2 0 2
1 4
3
Bn 5 32 1 2 4 .
1
30 2 4 4
6
7
Для n = 1, 2 находим
B1 1 4;
B2 1
4 2
1 12.
2 4
Записываем рекуррентное соотношение, связывающее три после!
довательных определителя: Bn 12 1 4 Bn 11 2 4 Bn .
Ему соответствует разностное уравнение
xn 12 1 4xn 11 2 4xn 3 0,
z2 1 4z 2 4 3 0.
Корни характеристического уравнения z1 = z2 = 2.
Общее решение имеет вид: xn 1 c12n 2 c2n2n 1 (c1 2 c2n)2n .
Находим коэффициенты с1 и с2:
n 1 1:
( c1 2 c2 )2 1 4 3
5 4 c1 1 c2 1 1.
n 1 2 : (c1 2 2c2 )4 1 126
О т в е т . Bn 3 1 n 4 12 2n .
14
Пример 13. Решить разностное уравнение третьего порядка
yn+3 – (a + 4)yn+2 + 4(a + 1)yn+1 – 4yn = 0;
y0 = b, y1 = 2b + 2c, y2 = 4b + 8c.
Выписываем характеристическое уравнение
z3 – (a + 4)z2 + 4(a + 1)z – 4 = (z – a)(z – 2)2 = 0
и находим его корни z1 = z2 = 2, z3 = a,
Общее решение имеет вид: yn 1 (c1 2 c2n)2n 2 c3an .
Находим коэффициенты с1, с2 и с3:
с1+ с3 = x0, (с1+ с2)2 + aс3 = x1, (с1+ 2с2)4 + a2с3 = x2.
Следовательно, с1 = b; с2 = c; с3 = 0.
О т в е т . yn = (b + cn)2n . Удивительно, но решение не зависит от
параметра a.
Пример 14. Разностные уравне!
T i (x)
ния и многочлены Чебышева.
1
Классические многочлены Чебы!
шева Tn(x) определены на интерва!
0,5
ле [–1, 1]. Они обладают рядом за!
x
мечательных свойств. В частности,
0,5
1
–1
–0,5
0
среди всех полиномов с фиксиро!
ванным коэффициентом при стар!
–0,5
шей степени х они имеют минималь!
–1
ную амплитуду на интервале [–1, 1]
(полиномы, наименее уклоняющи!
Рис. 1.2
еся от нуля). Первые четыре много!
члена имеют вид (рис. 1.2)
T0 (x) 1 1, T1(x) 1 x,
T2 (x) 1 2x2 2 1,
T3 (x) 1 4x3 2 3x.
Многочлены могут быть заданы с помощью рекуррентного соот!
ношения
Tn 11(x) 1 2xTn (x) 2 Tn 21(x).
Оно представляет собой линейное однородное разностное урав!
нение второго порядка. Ему соответствует характеристическое
уравнение
z2 1 2xz 2 1 3 0,
z1,2 3 x 4 x2 1 1.
15
Общее решение разностного уравнения имеет вид Tn (x) 1 c1z1n 2 c2z2n.
1
T1(x) 1 x находим c1 1 c2 1 . Отсю!
2
да получаем общую формулу для многочлена Чебышева
Из начальных условий T0 (x) 1 1,
1
13
Tn (x) 5 8 x 6 x2 7 1
2
2 6 1x 7
n
x2 7 1
2
n
4
9.
Вхождение радикалов в эту формулу – кажущееся, при раскрытии
скобок все они взаимно сокращаются. Коэффициент при старшем
члене равен 2n–1. Все нули многочленов Чебышева и все их экстрему!
мы сосредоточены на интервале [–1, 1].
Например, при n = 2 получаем
T2 (x) 3
1
2
1 2
x 4 2 x2 5 1 4 x2 5 1 4 x2 5 2 x2 5 1 4 x2 5 1 3 2x2 5 1.
2
Этот многочлен имеет два нуля, расположенных в точках 1 2 /2,
и минимум в точке 0 (рис. 1.2). Значения многочлена в точке мини!
мума и на краях интервала по абсолютной величине одинаковы и
равны единице.
Решение неоднородных разностных уравнений
Рассмотрим линейное неоднородное разностное уравнение
xn 1 k 1 an 21xn 1 k21 1 ... 1a1xk11 1 a0xk 2 fk,
где fk – известная функция дискретного аргумента.
Как и в случае дифференциальных уравнений, общее решение это!
го уравнения представляет собой сумму общего решения соответству!
ющего однородного уравнения и любого частного решения неодно!
родного уравнения. Частное решение обычно находят в том же виде,
что и функцию fk. Для определения постоянных коэффициентов, вхо!
дящих в общее решение неоднородного уравнения, используют на!
чальные или краевые условия.
Пример 15. Неоднородное разностное уравнение второго порядка
с комплексными корнями
xk 12 1 xk 11 1 xk 2 1;
x0 2 1;
x1 2 2.
Находим частное решение и корни характеристического полинома
11 2 3 i
1
xчаст 3 ; z2 4 z 4 1 3 0; z1,2 3
.
3
2
16
Общее решение неоднородного уравнения получаем как сумму ча!
стного решения и общего решения соответствующего однородного
уравнения
xk 2
1
1
1
3 c1 sin k 3 c2 cos k.
3
3
3
Коэффициенты с1, с2 находим из начальных условий: с1 = 8 3 ,
9
с2 = 2/3.
1.3. Системы линейных разностных уравнений
Матричная запись системы линейных однородных разностных
уравнений имеет вид
X(t + 1) = AX(t), X(0) 1 X0,
(1.15)
где X 1 R – вектор переменных; А – квадратная матрица постоян!
ных коэффициентов; X0 – вектор начальных условий.
Решение этой системы может быть записано в компактной сте!
пенной форме:
X(t) = At X0,
(1.16)
однако на деле процедура аналитического определения элементов
матрицы At представляет значительные трудности.
Для матриц простой структуры существует более удобный способ
записи решения. В соответствии с ним сначала надо найти собствен!
ные числа и собственные векторы матрицы A. Собственные числа zi
получаем как корни характеристического уравнения det(A – zE) = 0,
собственные векторы Hi находим из соотношений АHi = ziHi. Форму!
ла для общего решения имеет вид
n
X(t) = с1H1 z1t + ... + сnHn znt .
(1.17)
Постоянные сi определяются из начальных условий путем решения
системы линейных алгебраических уравнений X(0) = с1H1 + ... + сnHn.
В случае комплексных и кратных корней формула заметно услож!
няется. При наличии пары комплексных корней нужно выделять
вещественную и мнимую части соответствующих слагаемых в общем
решении и рассматривать их линейную комбинацию. При наличии
кратных корней в формуле (1.17) появляются полиномиальные мно!
жители. В частности, корню zi кратности k в общем решении соот!
ветствует слагаемое вида
1H t 1
k 1
1
2
3 1 3 Hk 11t 3 Hk zik,
(1.18)
17
где Н1, ..., Нк корневые векторы матрицы А, отвечающее собствен!
ному числу zi, причем Н1 – ее собственный вектор, а Нk – цикличес!
кий вектор, порождающий инвариантные подпространства матрицы
для указанного собственного числа.
К сожалению, элементы векторов Н2, ..., Нk зависят от началь!
ных условий и не могут быть рассчитаны заранее. Подробнее соот!
ветствующая теория и методика отыскания решения системы (1.15)
с помощью корневых векторов будет изложена в следующем разделе.
Пример 1. Найдем решение системы разностных уравнений вто!
рого порядка
x(t + 1) = x(t) – 3y(t),
y(t + 1) = x(t) + y(t).
1 1 z 13
2 0 и на!
1 11 z
ходим его корни z1,2 1 1 2 3 i . Собственные векторы матрицы A:
Выписываем характеристическое уравнение
1i 3 2
1 3i 3 2
H1 = 3 1 4 , H2 = 4 1 5 .
5
6
6
7
Запишем первую компоненту решения:
1
1
x 2 i 3 2t (cos t 3 i sin t)
1 3 i(1 3 3 i)t 2
x
1
2
3
3
6 ; или
H1 z1t = 5 y 6 4 5
7 8 7 (1 3 3 i)t 8
1
1
y 2 2t (cos t 3 i sin t).
3
3
Выделяем вещественную и мнимую части
1
x1 2 3 2t 3 sin t,
3
1
t
y1 2 2 cos t,
3
1
x2 2 2t 3 cos t,
3
1
t
y2 2 2 sin t
3
и строим из них общее решение
4 6 c 3 sin 3 t 7 c 3 cos 3 t 5
1
2
4x1 5
4 x2 5
t8
3
3 9.
X 1 t 2 c1 8 9 7 c2 8 9 2 8
9
3
3
y1 y2 c1 cos t 7 c2 sin t
8
9
3
3
Пример 2. Задача о размножении бактерий. В колонию бактерий,
численность которых удваивается каждую секунду, попадает вирус.
18
Каждую секунду он съедает одну бактерию, после чего размножается
делением. Какова судьба колонии бактерий, если их начальная чис!
ленность равна N?
Таким образом, сначала было N бактерий и один вирус. На первой
секунде вирус съел одну бактерию и размножился, число бактерий так!
же удвоилось, в результате стало 2(N – 1) бактерий и 2 вируса и т. д.
Р е ш е н и е . Обозначим через x и y текущее число вирусов и бакте!
1x 2
рий и рассмотрим вектор X 3 4 5 . Матричное описание этой задачи
6y7
имеет вид
1 2 02
1x(0) 2 1 1 2
, X(0) 3 4
3
.
X(t + 1) = AX(t), A 3 4
5
6
2
2
7
8
6 y(0) 75 64 N 75
Для получения решения воспользуемся сначала степенной форму!
лой (1.16). Последовательно возводя матрицу А в степени 2, 3, ..., n,
находим
1 1 02
3 1 04 3 1 4 t
2.
An = 2n 3
, X 1 t 2 5 A t X0 5 2t 6
6 5n 1 74
9 8t 1 7
69 N 8 t 7
Переходя к скалярной форме записи, получаем
x = 2t, y = (N – t) 2t.
Переменная y обращается в нуль при t = N, следовательно, коло!
ния вымрет через N секунд.
Второй способ решения опирается на вычисление собственных век!
торов матрицы А. В данном случае мы имеем дело со случаем крат!
ных собственных чисел z1 = z2 = 2, поэтому решение находят в форме
(1.18):
X(t) 1 (c1H1t 2 H2 )2t,
(1.19)
10 2
где H1 3 4 5 – единственный собственный вектор матрицы А.
61 7
Для определения корневого вектора H2 подставляем t = 0 и ис!
пользуем начальные условия
112
H2 3 X(0) 3 4 5 .
6N 7
Коэффициент с1 = –1 находим из равенства (A 1 2E)H2 2 2c1H1, ко!
торое получается после подстановки формулы (1.19) в исходную си!
стему уравнений.
19
Следовательно, решение имеет вид
1 30 4 3 1 4 2
3 1 4 t
X(t) 5 8 6 t 7 9 2t 5 2,
N 6 t 1 N что совпадает с полученным первым способом.
Для численного решения этой задачи можно воспользоваться од!
ним из математических пакетов. График решения, полученный с по!
мощью пакета MATLAB для N = 10, показан на рис. 1.3. Видно, что
на девятой секунде численность вирусов и бактерий совпала (точка
пересечения кривых), в результате чего на десятой секунде количе!
ство бактерий стало равным нулю. Дальнейшая часть кривой y(t)
физического смысла не имеет.
VIRUSES (x) and BACTERIES (y)
1200
x
1000
800
600
y
400
200
t
0
Dead
–200
–400
0
2
4
6
8
10
12
Рис. 1.3
Если начальная численность бактерий равна миллиону, то коло!
ния погибнет через миллион секунд (это примерно 11 сут. 13 час. и
47 мин.).
Отметим еще два пути решения этой задачи. Один из них связан
с переходом от системы двух разностных уравнений первого порядка
к одному уравнению второго порядка относительно y
y(t + 2) – 4y(t + 1) + 4y(t) = 0.
Второй путь – составить разностное уравнение для переменной
z = y/x. Оно получается очень простым z(t + 1) – z(t) = – 1 и имеет
линейное решение z(t) = N – t, которое обращается в нуль при t = N.
20
Пример 3. Задача Нарайаны (Индия, XIV в.). Сколько коров и
телок произойдет от одной коровы в течение 20 лет, если корова в
начале каждого года приносит телку, а телка дает такое же потом!
ство в начале года, достигнув трех лет.
Р е ш е н и е . Пусть xn – число коров, которые принесли потом!
ство в начале n!го года (число новорожденных телок тоже равно xn);
yn – число телок, которые принесут потомство через год; zn – число
телок, которые принесут потомство через два года.
Тогда общее число голов в стаде равно vn = 2 xn + yn + zn.
В силу условий задачи имеем систему разностных уравнений:
xn 11 1 xn 2 yn ,
x1 1 1,
yn 11 1 zn ,
y1 1 0,
zn 11 1 xn ,
z1 1 0,
vn 1 2xn 2 yn 2 zn .
(1.20)
Перепишем ее в матричной форме:
Xn 11 7 AXn ;
n 7 CXn ;
3x 4
X 7 5y6 ,
5z6
8 9
C 7 12 1
31 1 0 4
A 7 50 0 1 6 ,
51 0 0 6
8
9
12.
31 4
X1 7 50 6 .
50 6
8 9
(1.21)
Рассмотрим три подхода к решению этой системы разностных урав!
нений – скалярный, векторный и матричный.
Скалярный способ. Исключая переменные yn, zn из системы (1.20),
перейдем от нее к одному скалярному разностному уравнению тре!
тьего порядка
xn 13 1 xn 12 2 xn
с начальными условиями x1 = x2 = x3 = 1.
Будем вычислять значения xn по шагам, подставляя n = 1, 2, ..., 16.
Получим последовательность чисел 1, 1, 1, 2, 3, 4, 6, 9, 13, 19, 28,
41, 60, 88, 129, 189, 277, 406, 595, 872.
Заметим, что численность стада определяется формулой
vn 1 2xn 2 xn 11 2 xn 12.
Значит, через 20 лет стадо будет состоять из n = 2 · 872 + 595 +
+ 406 = 2745 голов.
При ручном счете этот способ быстрее всего приводит к цели.
21
Векторный способ. Используем представление решения через соб!
ственные числа li и собственные векторы Нi матрицы А (1.21):
Xn 1 c1H121n 3 c2H222n 3 c3H323n.
(1.22)
Характеристическое уравнение l3 – l2 – 1 = 0 имеет один веще!
ственный корень и пару комплексных корней
11 2 1,4656;
12,3 2 30,2328 4 0,7926i.
Собственные числа и собственные векторы матрицы А могут быть
вычислены в пакете MATLAB с помощью команды [H, L] = eig(A):
50,7710 30,3623 3 0,2924i 30,3623 4 0,2924i 6
H 9 70,3589 0,6784 4 0,0733i 0,6784 3 0,0733i 8 ,
7 0,5261 30,2160 4 0,5206i 30,2160 3 0,5206i 8
L 9 11,4656
30,2328 4 0,7926i
30,2328 3 0,7926i 2.
После подстановки этих значений в формулу (1.22) коэффициен!
ты сi находятся из начальных условий.
Если необходимо иметь не приближенные, а точные значения кор!
ней характеристического уравнения, можно воспользоваться тулбок!
сом SYMBOLIC пакета MATLAB. Выполним с его помощью символь!
ные вычисления корней:
y=solve(‘x^3 x^2 1 = 0’);
y(1)=1/6*(116+12*93^(1/2))^(1/3)+2/3/(116+12*93^(1/2))^(1/3) + 1/3
y(2)=1/12*(116+12*93^(1/2))^(1/3)1/3/(116+12*93^(1/2))^(1/3)+1/3+ 1/
2*i*3^(1/2)*(1/6*(116+12*93^(1/2))^(1/3)2/3/(116+12*93^(1/2))^(1/3))
Запись этих корней в обычной нотации имеет вид:
y1 2
y2,3 3 1
s2 1 2s 1 4
,
6s
s2 1 4s 2 4
s2 1 4
4i 3
.
12s
12s
где s 1 3 116 2 12 93 .
Этот способ можно рекомендовать для получения аналитического
решения при заданных начальных условиях и произвольном коли!
честве лет.
22
Кроме того, он позволяет дать ответ на вопрос о структуре стада
– процентном соотношении численности коров, телок!двухлеток
и телок!однолеток при больших n, а именно их численность будет
соотноситься, как компоненты первого собственного вектора мат!
рицы А, т. е. x : y : z = 0,7710 : 0,3589 : 0,5261. Другими словами,
коров будет около 46,6%, телок!двухлеток – около 21,6%, телок!
однолеток – около 31,8%. Полученное выше решение [x, y, z] =
= [ 872 406 595] для n = 19 характеризуется как раз таким соотно!
шением.
Матричный способ. Он опирается на использование степенных
формул
Xn 1 A n 11X1; vn 1 CA n 11X1.
При работе в пакете MATLAB этот способ удобнее всего, посколь!
ку после ввода матрицы А решение получается с помощью одной един!
ственной команды:
A=[1 1 0;0 0 1;1 0 0]; v=[2 1 1]*A^19*[1;0;0].
В результате получаем ответ: v = 2745 – это численность стада на
двадцатый год.
Ниже для справки приведены некоторые степени матрицы А = [1 1 0;
0 0 1; 1 0 0], вычисленные в пакете MATLAB:
11 1 1 2
A 2 5 31 0 0 4 ,
31 1 0 4
6
7
12 1 1 2
A 3 5 31 1 0 4 ,
31 1 1 4
6
7
13 2 12
A 4 5 31 1 14 ,
32 1 14
6
7
16 3 3 2
A 6 5 33 2 1 4 ,
3
4
64 3 2 7
1595 406 277 2
A18 5 3277 189 129 4 ,
3406 277 189 4
6
7
14 3 2 2
A 5 5 32 1 1 4 ,
33 2 1 4
6
7
1872 595 406 2
A19 5 3406 277 189 4 .
3595 406 277 4
6
7
Любопытно, что все они обладают скрытой симметрией, которая
становится явной после перестановки второй и третьей строк.
Динамика изменения поголовья стада по годам в обычном и лога!
рифмическом масштабах иллюстрируется рис. 1.4. Звездочками по!
казана общая численность стада.
Вычисления проводились с помощью следующей MATLAB!про!
граммы:
X0=[1;0;0]; A=[1 1 0;0 0 1;1 0 0]; X1=X0;X=zeros(3,20);
for i=1:20, X(:,i)=A*X1; X1=X(:,i); end
x=X(1,:); y=X(2,:); z=X(3,:); v=2x+y+z; t=1:20;
23
Herd of Cows
Herd of Cows
3000
8
plot(t,x,t,y,t,z,t,v,'*')
plot(t,log(x),t,log(y),t,
log(z),t, log(v),'*')
6
2000
4
1000
2
0
0
0
4
8
12
16
20
0
4
8
Time
12
Time
16
20
Рис. 1.4
plot(t,x,t,y,t,z,t,v,’*’),grid,
plot(t,log(x),t,log(y),t,log(z),t,log(v),’*’)
1.4. Синтез разностных уравнений
для получения заданных функций
При компьютерном моделировании часто возникает необходи!
мость генерировать функции времени. Так, можно указать следую!
щие примеры:
1) получение возмущающих функций при решении неоднородных
разностных уравнений;
2) воспроизведение переменных коэффициентов при моделирова!
нии линейных разностных уравнений;
3) получение тестовых воздействий для исследования систем ав!
томатического регулирования.
Проще всего генерировать сигналы и функции, которые являются
решениями линейных разностных уравнений с постоянными коэф!
фициентами. К таким функциям относятся at , tn , sin 1 t, cos 1 t и
их комбинации.
Общий прием генерирования заданной функции y(t) предполагает
отыскание определяющего линейного разностного уравнения
yn 1k 1 an 21yn 1k 21 1 1 1 a1yk11 1 a0yk 2 0,
для которого она является решением.
24
(1.23)
При этом могут использоваться методы, аналогичные тем, что
применяются для получения определяющих дифференциальных
уравнений [11].
Рассмотрим три таких метода: метод последовательных сдвигов
(аналог метода последовательного дифференцирования), метод ха!
рактеристического полинома и метод списков.
Метод последовательных сдвигов
Пусть задана функция y = f(t) и требуется синтезировать разно!
стное уравнение (1.23), решением которого она является. Для отыс!
кания этого уравнения рассмотрим наряду с функцией y(t) “сдвину!
тые” функции y(t), y(t + h), y(t + 2h), ..., причем операцию сдвига
будем выполнять до тех пор, пока очередная функция не окажется
линейной комбинацией предыдущих:
y(t + nh) = a0 y(t) + a1 y(t + h) + ... + an–1 y(t + (n – 1)h).
Найденное соотношение и будет представлять собой искомое раз!
ностное уравнение.
Для проверки системы функций на линейную зависимость можно
использовать стандартные математические критерии: вычисление
определителя Грама, сингулярных чисел, определение ранга и т.п.
Пример 1. Найдем разностное уравнение, имеющее своим реше!
нием заданную функцию y(t) = et.
Р е ш е н и е . Выбирая шаг h и рассматривая функции y(t), y(t + h) = ehet,
убеждаемся в их линейной зависимости: y(t + h) = a y(t), где а = eh.
Следовательно, искомое разностное уравнение имеет вид
y(t + h) – eh y(t) = 0, y(0) = 1.
Пример 2. Пусть требуется получить функцию y(t) = sint.
Р е ш е н и е . Для отыскания определяющего разностного уравне!
ния рассмотрим “сдвинутые” функции y(t + h) и y(t + 2h), где кон!
станта h задает шаг дискретизации аргумента. Используя формулы
тригонометрии, можем записать
y(t) 1 sin t
1
y(t 2 h) 1 cos h 3 sin t 2 sin h 3 cos t 42cos h .
y(t 2 2h) 1 cos2h 3 sin t 2 sin 2h 3 cos t
1
(1.24)
Домножим эти равенства на коэффициенты, указанные справа, и
сложим:
y(t 1 2h) 2 2cos hy(t 1 h) 1 y(t) 3
2
3 (1 2 2cos h 1 cos2h)sin t 1 (sin2h 2 2sin h cos h)cos t.
25
Выражения в скобках тождественно равны нулю, поэтому окон!
чательно получаем
y(t 1 2h) 2 2y(t 1 h)cos h 1 y(t) 3 0.
Это – искомое разностное уравнение второго порядка. Начальные
условия для него находим из соотношений (1.24): y(0) = 0; y(h) = sin h.
При h = 0,1 получаем
y(t 1 0,2) 2 1,99008y(t 1 0,1) 1 y(t) 3 0,
y(0) 3 0,
y(h) 3 0,099833.
П р о в е р к а . Чтобы убедиться в правильности построенного урав!
нения, решим его при начальных условиях y(0) = 0, y(h) = sinh.
Выписываем характеристический полином и находим его корни:
z2 1 2z cos h 2 1 3 0,
z1,2 3 cos h 4 i sin h.
Вычисляем модуль и аргумент корней: 1 2 1,
общее решение
y(t) 1 c1 cos kh 2 c2 sin kh;
3 2 h и записываем
k 1 0, 1, 2, ... .
Подставляя начальные условия, находим, что с1 = 0, с2 = 1, т. е.
при t = kh полученное решение совпадает с требуемым.
Пример 3. Найдем методом последовательных сдвигов определя!
ющее разностное уравнение для функции y(t) = tet.
Р е ш е н и е . Выберем шаг h и перейдем к дискретному времени t = kh:
yk 1 khekh , yk11 1 (k 2 1)he(k11) h .
Раскрывая скобки, выражение для сдвинутой на такт функции
yk+1 можно переписать в форме
1
yk11 3 ykeh 15 1 4 26.
k8
7
Отсюда получаем разностное уравнение первого порядка
yk11 3 ayk 4 0,
1
a 4 16 1 5 27 eh ,
k9
8
y0 4 0.
Это нестационарное уравнение, так как коэффициент а зависит от
номера такта.
Чтобы получить стационарное уравнение, выпишем дополнитель!
но функцию yk12 :
2
yk12 3 (k 4 2)he(k12)h 3 15 1 4 26 e2h yk.
k8
7
26
Cовокупность трех функций kekh, (k + 1)ekh, (k + 2)ekh линейно
зависима, поэтому дальнейших сдвигов не требуется.
Сложим равенства
yk 1 khyk
1
yk11 1 yke h 2 ykeh
k
h 2
yk12 1 yke 2 yke2h
k
e2 h
32eh ,
1
домножив их на коэффициенты, указанные справа.
При этом правая часть обращается в нуль, и мы получаем искомое
разностное уравнение
yk12 1 2eh yk11 2 e2h yk 3 0,
y0 3 0,
y1 3 heh .
В частности, при h = 0,1 оно принимает вид
yk12 1 2,21yk11 2 1,22yk 3 0, y0 3 0,
y1 3 0,110517.
Метод характеристического полинома
При использовании этого метода требуется знать соответствие
между корнями характеристического уравнения и решениями соот!
ветствующего разностного уравнения (табл. 1).
Таблица 1
№
Вид корня
Характер решения
1
z1 = 1
c
2
z1 = a
cat
3
z1 = ... = zk = 1
c1 + c2 t + ... +ck tk – 1
4
z1 = ... = zk = a
(c1 + c2 t + ... +ck tk – 1 )at
5
z1 , 2 = a ± ib
rt(c1sinjt + c2 cosjt), r2 = a2 + b2; tgj = b/a
6
z1 , 2 = z3 , 4 = a ± ib
rt [(c1 + c2 t)sinjt + (c3 + c4 t)cosjt]
По заданному виду функций определяются необходимые корни
z1, ..., zn и составляется характеристическое уравнение
(z – z1)(z – z2) ... (z – zn) = 0.
От него легко перейти к разностному уравнению
y(t + nh) + an–1 y(t + (n – 1)h) + ... + a1 y(t + h) + a0 y(t) = 0,
27
где коэффициенты ai рассчитываются по известным корням zi с помо!
щью формул Виета
a0 1 (21)n z1 1zn, 1, an 11 1 2(z1 3 1 3 zn ).
Пример 4. Найдем определяющее разностное уравнение для фун!
кции y(t) = te– t методом характеристического полинома.
Р е ш е н и е . Перейдем к дискретному времени t = kh, h = 0, 1, 2, ... .
Нам требуется получить функцию вида yk 1 (c1 2 c2k)zkh , следователь!
но, характеристический полином должен иметь кратный ко!
рень z1 1 z2 1 e 1h . Восстанавливая полином по корням, получаем
(z 1 z1 )(z 1 z2 ) 2 z2 1 2e 1h z 3 e 12h 2 0.
Ему соответствует разностное уравнение второго порядка
yk22 1 2e1h yk 21 2 e 12h yk 3 0
c начальными условиями y0 1 0, y1 1 he 1h .
При h = 0,1 оно принимает вид
yk12 1 1,81yk11 2 0,819yk 3 0, y0 3 0,
y1 3 0,0904837.
Пример 5. Пусть необходимо воспроизвести функцию
1
y(t) 2 5t3t 3 2t sin t.
4
Р е ш е н и е . При t = kh первому слагаемому соответствуют корни
z1 1 z2 1 3h , а второму – корни z3,4 1 2 3 i4. Коэффициенты a и b нахо!
дят из соотношений 12 2 32 4 22h; tg5 4 3/ 1 4 1, откуда
z3,4 1 2 (1 2 i).
Характеристическое уравнение имеет вид
1z 3 3 2 1z
h
2
2
2
3 24z 5 62 7 0 .
После подстановки 1 2 2 , 32 2 4h и раскрытия скобок получаем
z4 1 (2 2 2 3h )z3 2 (4h 2 9h )z2 1 (12h 2 2 2 9h )z 2 36h 3 0.
Следовательно, искомое разностное уравнение имеет вид
yk1 4 1 (3h 2 2 2)yk13 2 (4h 2 9h )yk12 1 (12h 2 2 2 9h )yk 11 2 36h yk 3 0.
Начальные условия равны значениям функции y(t) при t = 0, h,
2h, 3h.
28
Метод списков
В некоторых случаях требуется воспроизвести не одну, а несколь!
ко функций времени y1 = f1(t), ..., ym = fm(t). Для этого можно ис!
пользовать описанные выше способы, что приводит к получению не!
скольких определяющих разностных уравнений разных порядков.
Однако при этом не гарантируется минимальность модели в смысле
числа используемых элементов задержки.
Для получения минимальной реализации будем искать ее описа!
ние в виде системы определяющих разностных уравнений первого
порядка
Xk11 1 AXk,
Yk 1 CXk,
X0 ,
(1.25)
где Х – вектор переменных состояния; Y – вектор выходных сигналов,
реализующих заданные функции времени; Х0 – вектор начальных ус!
ловий; А и С – постоянные матрицы, подлежащие определению.
Алгоритм получения системы (1.25) минимального порядка, обес!
печивающей выполнение равенств y1 = f1(t), ..., ym = fm(t), состоит в
следующем.
Ш а г 1. Составляем список S1 линейно независимых элемен!
тарных функций 1i (t), i 2 1, q, через которые линейно выражают!
ся заданные функции f1(t), ..., fm(t). В частности, можно принять
j1 = f1(t), ..., jm = fm(t).
Ш а г 2. Сдвигаем функции списка S1 на такт h, раскрываем выра!
жения j1 (t + h), ..., jq = (t + h) и составляем список S2 новых функ!
ций 1i (t), i 2 q 3 1, r, появившихся после такого сдвига и не выра!
жающихся линейно через функции списка S1.
Ш а г 3. Сдвигаем функции списка S2 на такт h, раскрываем вы!
ражения jq+1 (t + h), ..., jr = (t + h)и составляем список S3 новых
функций 1i (t), i 2 r 3 1, s, появившихся после этого сдвига и не вы!
ражающихся линейно через функции списков S1 и S2.
Шаги, аналогичные шагам 2 и 3, повторяются до тех пор, пока ока!
жется, что после очередного сдвига на k!м шаге новых функций, линей!
но независимых от предыдущих, не появилось. Если такого шага не
наступит, то это означает, что системы определяющих разностных урав!
нений вида (1.25) для заданных функций не существует.
Шаг k + 1. Размерность системы (1.25) берется равной общему
числу n функций 1i (t), вошедших в списки S1, S2, ..., Sk. Все эти
функции линейно независимы, чем гарантируется минимальное чис!
ло уравнений в системе. Cами эти функции принимаются за перемен!
ные состояния системы (1.25)
x1 = j1(t), x2 = j2(t), ..., xn = jn(t).
(1.26)
29
Записывая равенства (1.26) для момента времени t + h и выражая
их правые части через x1(t), ..., xn(t) (это возможно, так как они дол!
жны линейно зависеть от 11(t), ..., 1n (t) ), получим систему разно!
стных уравнений Xk 11 1 AXk. Матрицу С находим, выражая функ!
ции fi(t) через 11(t), ..., 1n (t) . Компоненты вектора начальных усло!
вий Х0 получаем из равенств (1.26) при t = 0.
Проиллюстрируем описанный алгоритм на примере.
Пример 6. Пусть требуется найти определяющие разностные урав!
нения в форме (1.25) для функций
x 1 t sin t,
y 1 t cos t,
t 1 0, h, 2h, ...,
представляющих параметрическое описание спирали Архимеда.
Решение.
Ш а г 1. Принимая 11 2 x, 12 2 y, получаем список S1:
S1 1 (t sin t, t cos t) 1 (21(t), 22 (t)).
Ш а г 2. Выполняем сдвиг функций 11(t), 12 (t) на такт
x(t 1 h) 2 (t 1 h)sin(t 1 h) 2 (t 1 h)(cos h sin t 1 sin h cos t),
y(t 1 h) 2 (t 1 h)cos(t 1 h) 2 (t 1 h)(cos t cos h 3 sin t sin h).
Новыми функциями здесь являются sint и cost, поэтому
S2 1 (sin t, cos t) 1 (23 (t), 24 (t)).
Ш а г 3. Выполняем сдвиг функций 13 (t), 14 (t) на такт
13 (t 2 h) 3 cos h sin t 2 sin h cos t, 14 (t 2 h) 3 cos h cos t 4 sin h sin t.
На данном шаге новых функций не появилось, поэтому дальней!
шие сдвиги не требуются.
Ш а г 4. Размерность системы (1.25) берем равной 4, равенства
(1.26) имеют вид
x1 1 t sin t, x2 1 t cos t, x3 1 sin t, x4 1 cos t.
После сдвига на такт с учетом формул, полученных на предыду!
щих шагах, приходим к следующей системе определяющих разно!
стных уравнений:
x1(t 1 h) 2 x1(t)cos h 1 x3 (t)h cos h 1 x2 (t)sin h 1 x4 (t)h sin h ,
x2 (t 1 h) 2 x2 (t)cos h 3 x1 (t)sin h 1 x4 (t)h cos h 1 x3h sin h ,
x3 (t 1 h) 2 x3 (t)cos h 1 x4 (t)sin h ,
x4 (t 1 h) 2 3x3 (t)sin h 1 x4 (t)cos h ,
30
y1(t) 1 x1(t),
y2 (t) 1 x2 (t).
Ей соответствует описание Xk11 1 AXk,
Yk 1 CXk с матрицами
1 cos h sin h h cos h h sin h 2
10 2
3 5 sin h cos h 5h sin h h cos h 4
30 4
11 0 0 0 2
A63
; C63
; X0 6 3 4 .
0
0
cos h
sin h 4
0
70 1 0 0 48
3
4
3 4
0
cos h 48
5 sin h
37 0
371 48
Моделируя его в любом математическом пакете и строя зависи!
мость y1 = f(y2), получим на дисплее график спирали Архимеда.
1.5. Задачи и упражнения
1. Показать эквивалентность уравнений (1.5) и (1.6).
2. Составить линейное однородное разностное уравнение с посто!
янными коэффициентами, решением которого является:
константа x(t) = a ; линейная функция x(t) = at + b ;
парабола x(t) = at2 ; косинусоида x(t) = соst.
Аргумент везде изменяется дискретно t = 0, 1, 2, ...
3. Как изменится решение однородного разностного уравнения,
если изменить все знаки начальных условий на противоположные?
Если вдвое изменить величину начальных условий?
4. Решить уравнение
x(t 1 2) 1 6x(t 1 1) 1 9x(t) 2 0 при
x(0)= 1 1, x(1) 2 2 .
Р е ш е н и е . Записываем характеристическое уравнение z2 + 6z +
+ 9 = 0. У него кратные корни z1,2 = – 3, следовательно, общее реше!
ние имеет вид x = (c1t + c2)(–3)t.
Находим постоянные коэффициенты: с2 = –1; c1 = 1/3.
О т в е т . x = –(–3)t + t(–3)t / 3 .
5. Найти общее решение разностного уравнения yk11 1 yk 2 2yk 21 3 0.
О т в е т . yk 3
1 2 2 1 C sin k4 5 C cos k42, 4 3 arctg
k
1
2
7.
6. Найти общее решение разностного уравнения
13yk 11 1 4yk 1 yk21 2 0.
О т в е т . yk 3
1
13
2 1 c sin k4 5 c cos k42, 4 3 6arctg 32 .
k
1
2
7. Найти решение разностного уравнения yk12 1 yk 2 0, y0 1 2,
y1 1 1.
31
1k
1k
3 sin .
2
2
8. Найти решение разностного уравнения
О т в е т . yk 2 2cos
yk11 1 4yk 2 yk21 2 6yk22 3 0, y0 3 6, y1 3 12,
y4 3 276.
О т в е т . yk 1 (21)k 3 2k 11 3 3k 11.
9. Найти решение системы разностных уравнений
x(t + 1) = x(t) – 3y(t)
y(t + 1) = 3x(t) + y(t).
У к а з а н и е . Собственные числа и собственные векторы матрицы
A в данном случае имеют вид:
112
112
z1,2 1 1 2 3i, H1 3 4 5 , H2 = H2 3 4 5 .
6
i
6i 7
7 8
10. Найти частное решение неоднородного уравнения
6yk 11 1 5yk 2 yk 21 3 2k.
О т в е т . yk 1 2k2k.
11. Найти общее решение неоднородного уравнения
5
yk11 1 yk 2 yk21 3 cos k.
2
О т в е т . yk 1 c12k 2 c221k 2
cos k
.
5
2
12. Найти общее решение неоднородного уравнения
2cos1 3
yk11 1 yk 2 5yk21 1 3yk 22 3 1.
k 1 2
О т в е т . yk 1 c1 2 c2k 2 c3 (33) 2 k .
8
13. Записать аналитические выражения для переменных x, y, z, v в
задаче Нарайаны.
32
2. ПРИМЕНЕНИЕ ЖОРДАНОВЫХ ЦЕПОЧЕК ВЕКТОРОВ
ПРИ РЕШЕНИИ ДИФФЕРЕНЦИАЛЬНЫХ
И РАЗНОСТНЫХ УРАВНЕНИЙ
2.1. Корневые векторы и циклические пространства
Понятие собственных векторов является одним из центральных в
линейной алгебре и находит многочисленные применения во многих
областях математики, физики и других естественных наук. Они ис!
пользуются при решении линейных дифференциальных и разностных
уравнений, при синтезе регуляторов в теории управления, при расче!
те форм упругих механических конструкций (балок, мостов, лета!
тельных аппаратов).
Процедуры вычисления собственных чисел и собственных векто!
ров входят в состав всех профессиональных математических паке!
тов, таких как MATHCAD, MATLAB, MAPLE и т. д.
В линейной алгебре собственные векторы матрицы А вводятся как
решения уравнения AHi 1 2iHi , где 1i – собственные числа матрицы.
Если все собственные числа различны, то матрица обладает полной
системой собственных векторов, из которых может быть составлен
базис пространства. В этом случае матрица А может быть приведена
преобразованием подобия к диагональному виду.
Ситуация осложняется, если среди собственных чисел матрицы
имеются кратные. В этом случае число собственных векторов может
быть меньше размерности пространства. Чтобы сформировать базис
пространства, к ним приходится добавлять дополнительные векто!
ры, которые называются обобщенными собственными или корневы!
ми векторами. Жорданова каноническая форма матрицы в этом слу!
чае имеет уже не диагональный, а клеточно!диагональный вид.
Процедуры вычисления размеров этих клеток и соответствующих
им корневых векторов достаточно сложны и плохо обусловлены. По
этой причине они отсутствуют во многих математических пакетах, в
том числе в таком мощном пакете, как MATLAB.
В то же время необходимость в их вычислении возникает при ре!
шении ряда задач. Например, известная форма записи общего реше!
1 1 AX :
ния системы дифференциальных уравнений X
33
X(t) 1 c1H1e11t 2 1 2 cnHn e1nt ,
где 1i, Hi – собственные числа и собственные векторы матрицы А,
справедлива только для матриц А с попарно различными собствен!
ными числами.
Возникает вопрос об обобщении этой формулы на случай кратных
собственных чисел. Очевидно, что оно потребует введения корневых
векторов и построения жордановых цепочек векторов, зависящих от
начальных условий системы. Аналогичный вопрос относится к систе!
мам разностных уравнений вида Xk 11 1 AXk и формулам (1.17), (1.18).
Ниже дается ответ на оба эти вопроса и выводятся формулы для
решения дифференциальных и разностных уравнений для матрицы
А с собственными числами произвольной кратности в терминах кор!
невых векторов. Сначала изложим необходимые математические све!
дения, в частности дадим определения минимальных многочленов
вектора, корневых векторов и жордановых цепочек векторов.
Минимальный многочлен вектора
Пусть дана матрица А. Возьмем произвольный вектор b и соста!
вим ряд векторов b, Ab, A2b, ... . Пусть m – наименьшее число, при
котором очередной вектор будет линейно зависим от предыдущих
A m b 1 23 m 11A m 11b 2 3m 12A m 12b 2 1 2 30 b.
(2.1)
Вводя многочлен 1(2) 3 2m 4 5 m 112m 11 4 1 4 5 0 , это равенство мож!
но переписать так:
1(A)b 2 0.
(2.2)
Всякий многочлен 1(2) , для которого имеет место равенство (2.2),
называется аннулирующим для вектора b (относительно оператора
А). Из всех аннулирующих многочленов для вектора b построенный
нами многочлен имеет наименьшую степень со старшим коэффици!
ентом 1. Он называется минимальным аннулирующим многочленом
вектора b или просто минимальным многочленом вектора b.
В теории управления матрица R 1 [b, Ab, ..., A n 11b] известна как
1 1 AX 2 bu. Число m совпадает с ее
матрица управляемости системы X
рангом rank R = m и называется индексом управляемости системы.
Если rank R < n, где n – размер матрицы А, то система будет неуправля!
емой – это известный критерий неуправляемости Калмана.
С помощью произвольного входного сигнала u(t) такую систему
можно вывести из состояния покоя, но она будет двигаться только в
m!мерном подпространстве n!мерного пространства состояний. Оно
называется подпространством управляемости.
34
Знание минимального многочлена вектора b позволяет сформиро!
вать такой входной сигнал u(t) конечной длительности, реакция на
который после его окончания будет равна нулю, т. е. он аннулируется
системой. Простой пример – последовательность из m + 1 прямоуголь!
ных импульсов, амплитуды которых равны коэффициентам аннулиру!
ющего многочлена. На этом принципе, в частности, основан один из
методов тестового диагностирования систем автоматического управле!
ния (метод Шрайбера или метод комплементарного сигнала1).
Циклические пространства
Вновь рассмотрим векторы b, Ab, ..., A m 11b, где m – порядок
минимального многочлена вектора b. Они линейно независимы и об!
разуют базис m!мерного пространства Rm, которое называется цик3
лическим.
Оператор А переводит первый вектор этого базиса во второй, вто!
рой – в третий и т. д. Последний базисный вектор переводится в ли!
нейную комбинацию базисных векторов согласно равенству (2.1).
Значит, и произвольный вектор из Rm будет переводиться операто!
ром А в Rm . Таким образом, циклическое пространство инвариантно
относительно действия оператора А.
Минимальный многочлен первого базисного вектора одновремен!
но будет минимальным многочленом пространства Rm . Отсюда
вытекает следующий критерий цикличности пространства.
Теорема. Пространство циклично относительно оператора А тог!
да и только тогда, когда число его измерений совпадает со степенью
его минимального многочлена.
Например, минимальный многочлен единичной матрицы имеет
вид 1(2) 3 2 4 1. Здесь m = 1, т. е. критерий цикличности не выполня!
ется даже для двумерной плоскости.
Циклическое пространство может расщепляться на циклические
подпространства. Простейшие циклические подпространства – соб!
ственные направления матрицы, их размерность равна единице. В
случае кратных собственных чисел появляются циклические под!
пространства большей размерности.
2.2. Жордановы цепочки векторов
Пусть у матрицы А размера n 1 n имеется собственное число l1 крат!
ности k и ему отвечает циклическое подпространство размерности k с
1
Мироновский Л.А. Диагностирование линейных систем методом компле!
ментарного сигнала //Приборы и системы. Управление, контроль, диагности!
ка. 2002. № 5. С. 52–57.
35
порождающим вектором b. Для этого вектора многочлен j(l) = (l – l1)k
будет минимальным многочленом.
Рассмотрим векторы
b1 = (A – l1E)k–1b, b2 = (A – l1E)k–2b, ..., bk = b.
Умножая их на A – l1E и учитывая, что (A – l1E)kb = 0, получаем
цепочку равенств:
(A – l1E)b1 = 0, (A – l1E)b2 = = b1 , ..., (A – l1E)bk = bk–1,
или
A b1 = l1b1, A b2 = l1b2 + b1 , ..., A bk = l1bk + bk–1.
(2.3)
Линейно независимые векторы b1, b2, ..., bk = b образуют цепочку
векторов, которая называется жордановой. В совокупности все жор!
дановы цепочки матрицы А образуют жорданов базис в Rn.
Корневые векторы
Векторы жордановой цепочки принадлежат множеству так назы!
ваемых корневых векторов.
О п р е д е л е н и е [7, с. 147]. Вектор Н называется корневым век3
тором матрицы А, принадлежащим ее собственному числу l, если
1 A 3 4E 2k H 5 0.
Наименьшее число k, при котором выполняется это равенство,
называется высотой корневого вектора.
Заметим, что при k = 1 получаем обычное определение собственно!
го вектора матрицы. Таким образом, понятие корневого вектора обоб!
щает понятие собственного вектора в случае кратных собственных
чисел.
Векторы жордановой цепочки – это совокупность корневых век!
торов различной высоты (от 1 до k включительно), получаемых на
основе одного порождающего вектора путем его последовательного
умножения на матрицу A 0 1 A 2 3E.
Разным порождающим векторам будут отвечать разные жордано!
вы цепочки, заканчивающиеся, однако, одним и тем же собственным
вектором b1 (с точностью до скалярного множителя).
Перечислим некоторые свойства корневых векторов.
1. Корневые векторы высоты k образуют линейное пространство
размерности k – инвариантное подпространство матрицы А.
2. Умножение корневого вектора высоты k > 1 на матри!
цу A 0 1 A 2 3E уменьшает его высоту на 1.
3. Следующие операции не меняют высоты корневого вектора:
– умножение на число;
36
– умножение на матрицу А;
– добавление линейной комбинации корневых векторов меньшей
высоты.
Отметим, что в некоторых книгах по линейной алгебре вместо на!
звания “корневой вектор высоты k” используется термин “обобщен!
ный собственный вектор ранга k”.
Опишем еще один подход к определению корневых векторов [4].
Пусть у матрицы А собственному числу l кратности два отвечает
единственный собственный вектор Н,
1G
так что АН = lН. Предположим, что в
результате малой вариации матрицы А
получается матрица А1 с близкими соб!
ственными числами l и l + e, которым
H
отвечают близкие собственные векторы
H+ 1 G
Н и Н + dG (рис. 2.1). Здесь e и d – вели!
чины одного порядка малости. Пред!
ставляется разумным в качестве обоб!
щенного собственного вектора матрицы
Рис. 2.1
А взять вектор G, вернее тот вектор, к
которому он стремится в пределе при 1 2 0.
Чтобы получить уравнение для этого вектора, рассмотрим равен!
ство
A1(H + dG) = (l + e)(Н + dG).
Раскрывая скобки, пренебрегая произведением малых величин и
учитывая равенство A1H = lН, получим
1
A1G 2 3G 4 H.
5
В пределе при 1 2 0 имеем A1 2 A, 1 2 const , так как величины e
3
и d одного порядка малости. Таким образом, вектор G является ненуле!
вым решением уравнения (A 1 2E)G 3 cH, которое при с = 1 имеет тот же
вид, что и уравнения (2.3) для векторов жордановой цепочки. Величи!
на константы с принципиальной роли не играет, поскольку обобщен!
ные собственные векторы, как и обычные собственные векторы, опреде!
ляются с точностью до произвольного коэффициента.
Пример 1. Найдем жорданову цепочку векторов матрицы
12 3 3 2
A 5 30 2 3 4 .
30 0 2 4
6
7
37
Она имеет единственный собственный вектор Н1 = [1 0 0]Т, отве!
чающий собственному числу l = 2 кратности k = 3. Следовательно,
жорданова цепочка будет состоять из трех векторов. Корневые век!
торы находим из уравнений А0Н2 = Н1, А0Н3 = Н2, где
1
T
10 3 32 3h22 1 3h23 2 1, 3h23 2 0 3 H2 2 [c2 3 0] ,
A0 5 A 6 7E 5 30 0 34 :
30 0 04 3h 1 3h 2 c , 3h 2 1 3 H 2 [c c2 1 ].
32
33
1
33
3
3
8
9
3
3 9
Здесь с2 и с3 – произвольные постоянные.
Если собственный вектор Н1 умножить на произвольный коэффи!
циент с1 ¹ 0, то проводя аналогичные выкладки, получим
1c1 2
H1 5 3 0 4 ,
304
7 8
1 c2 2
3c 4
H2 5 3 1 4 ,
334
73 0 84
1
2
3 c3 4
3c c 4
H3 5 3 2 6 1 4 .
33 94
3 c1 4
37 9 48
Варьируя коэффициенты с1, с2 , с3, можно получать различные
жордановы цепочки векторов. В частности, при с1 = 9, с2 = с3 = 0
имеем
19 2
H1 5 30 4 ,
3 4
70 8
10 2
H2 5 3 3 4 ,
3 4
70 8
102
H3 5 3 614 .
3 4
718
При решении практических задач жорданову цепочку бывает удоб!
нее строить не с начала, а с конца, выбирая в качестве Н3 конкрет!
ный корневой вектор высоты 3. Тогда неопределенность полностью
устраняется и векторы Н2 , Н1 находятся однозначно.
Например, взяв все компоненты вектора Н3 единичными, полу!
чим
112
H3 5 314 ,
314
6 7
16 2
H2 5 A 0H 3 5 3 3 4 ,
30 4
6 7
19 2
H1 5 A 0H2 5 30 4 .
30 4
6 7
Во всех случаях матрица корневых векторов H = [H1 H2 H3] ока!
зывается треугольной, что следует из треугольности матрицы А.
38
2.3. Дифференциальные уравнения
и жордановы цепочки векторов
Дифференциальные уравнения с кратными корнями
В теории линейных дифференциальных уравнений для записи об!
щего решения системы уравнений
1 1 AX,
X
X(0) 1 X0,
(2.4)
где А – постоянная матрица размера n 1 n ; X – вектор переменных,
наряду с формулой X(t) 1 e AtX0 используется формула [5, 13]
X(t) 1 c1H1e11t 2 1 2 cnHn e1nt ,
(2.5)
где li и Hi – собственные числа и собственные векторы матрицы А;
сi – произвольные постоянные, зависящие от начальных условий.
К сожалению, эта формула справедлива только для матриц А про!
стой структуры, когда среди собственных чисел нет кратных. Пока!
жем, как ее можно обобщить на случай кратных собственных чисел,
используя понятие жордановых цепочек векторов.
Пусть l – собственное число матрицы А кратности k. Будем искать
соответствующие ему компоненты решения в форме экспоненты e1t ,
домноженной на полином от t степени k – 1:
1
t2
t k11 2 2t
X(t) 3 5 H0 4 H1t 4 H2 4 1 4 Hk11
(2.6)
6e ,
2!
(k 7 1)! 9
8
где Hi – векторы, подлежащие определению. Факториальные коэф!
фициенты введены для удобства дальнейших выкладок и принципи!
альной роли не играют.
Подставим это выражение в систему (2.4). После дифференциро!
вания и сокращения на множитель e1t получим:
H1 3 H2t 3 1 3 Hk 11
1
tk 12
tk11 2
3 4 6 H0 3 H1t 3 1 3 Hk11
75
(k 8 2)!
(k 8 1)! 9
1
tk11 2
5 A 6 H0 3 H1t 3 1 3 Hk11
7.
(k 8 1)! 9
Приравняем коэффициенты при одинаковых степенях t в правой
и левой частях уравнения:
AH0 1 2H0 3 H1,
AH1 1 2H1 3 H2, ...,
AHk12 1 2Hk 12 3 Hk11,
AHk 11 1 2Hk11.
(2.7)
39
Мы получили уравнения, которые совпадают, с точностью до обо!
значений, с цепочкой уравнений (2.3). Следовательно, они описыва!
ют жорданову цепочку векторов, отвечающих кратному собственно!
му числу матрицы А. При этом Hk–1 – это собственный вектор матри!
цы А, а Н0 – корневой вектор высоты k.
Для определения вектора Н0 используются начальные условия
системы (2.4). Особенно просто это делается в случае, если у матри!
цы А, кроме собственного числа l кратности k = n, нет других соб!
ственных чисел (матрица с одной жордановой клеткой). Тогда фор!
мула (2.6) описывает общее решение системы, а корневой вектор Н0
совпадает с вектором начальных условий Н0 = Х0, что прямо следует
из равенства (2.6) при t = 0.
Пример 1. На рис. 2.2 приведена структурная схема системы, пред!
ставляющей собой соединение трех апериодических звеньев с одина!
ковыми передаточными функциями.
c
a
b
1
p –2
x3
3
1
x2
p –2
1
3
x1
p–2
Рис. 2.2
Требуется найти формулы, описывающие свободное движение этой
системы из начальных условий: х1 (0) = а, х2 (0) = b, х3 (0) = c.
Р е ш е н и е . Перейдем к описанию в пространстве состояний вида
(2.4):
12 3 3 2
1 5 30 2 34 X,
X
30 0 2 4
6
7
1a 2
X0 5 3 b 4 .
3c4
6 7
(2.8)
Матрица А этой системы такая же, как в предыдущем примере,
она имеет единственный собственный вектор и одно собственное чис!
ло l = 2 кратности три. Следовательно, решение системы может быть
записано в форме (2.6) при k = 3:
1
t2 2
X(t) 3 5 H0 4 H1t 4 H2 6 e2t .
28
7
40
(2.9)
Полагая t = 0, устанавливаем, что H0 = X0. Этот вектор порождает
жорданову цепочку векторов Н0, Н1, Н2, связанных формулами
(2.7):
1a 2
H0 5 3 b 4 , H1 5
3c 4
6 7
10 3 3 2
30 0 3 4 ,
30 0 0 4
6
7
2b 1 c 3
20 3 3 3
2c3
H0 6 3 4 c 5 , H2 6 40 0 3 5 , H1 6 9 405 .
4 0 5
40 0 0 5
405
7
8
7
8
7 8
Подставляя эти значения в формулу (2.9), получаем окончатель!
ный результат:
1a 5 3(b 5 c)t 5 4,5 ct2 2
3 1a 2
1b 5 c 2
1c 2 4
9 2 6 7 9 2t 6
7 e2t .
8
6
7
6
7
X(t) b 5 3t c 5 t 0
e b 5 3ct
6
7
6 0 7 2 60 7 99
88 6 c 7
c
6
7
Заметим, что описанную процедуру можно рассматривать как спо!
соб вычисления матричной экспоненты. Действительно, сравнивая
последнюю формулу с представлением решения в виде X(t) 1 e At X0,
легко заключаем, что в данном случае
e
At
11 3t 4,5t2 2
5 30 1
3t 4 e2t .
3
4
1 4
360 0
7
Завершая рассмотрение примера, приведем запись полученного
решения в скалярной форме:
x1 1 (a 2 3(b 2 c)t 2 4,5ct2 )e2t; x2 1 (b 2 3ct)e2t; x3 1 ce2t .
Выделение модальных компонент
Одно из эффективных применений формулы (2.5) связано с опре!
делением начальных условий, позволяющих подавлять в решении
все модальные компоненты, кроме одной. Это достигается, если в
качестве начальных условий брать один из собственных векторов
матрицы А. Например, полагая X0 = Hi, получим X(t) 1 Hie1i t , т. е. в
решении будет присутствовать только одна модальная компонента.
Аналогичный результат имеет место и в случае формулы (2.6),
если в качестве начальных условий брать корневые векторы различ!
41
ной высоты. Однако решение здесь оказывается неоднозначным из!
за возможного произвола в выборе корневых векторов.
Поясним эту ситуацию на примере матрицы А с единственной жор!
дановой клеткой.
Пример 2. Рассмотрим схему из трех последовательно соединен!
ных апериодических звеньев, показанную на рис. 2.3.
y0
z0
1
p +a
z
1
x0
y
p +a
1
x
p +a
Рис. 2.3
Ее описание в пространстве состояний имеет вид
03
2x0 3
2 1a 1
1 1 AX; A 6 4 0 1a 1 5 , X 6 4 y 5 .
X
0
4 05
40
0 1a 58
7
7 z0 8
Матрица А обладает единственным собственным числом l = – а
кратности три.
Выходной сигнал схемы при произвольных начальных условиях
описывается формулой
x 1 c1e 1at 2 c2te 1at 2 c3t2e 1at.
Требуется указать начальные условия x0, y0, z0, при которых в
этой формуле будет оставаться только одно из слагаемых (первое,
второе или третье соответственно).
Р е ш е н и е . Покажем, что искомыми начальными условиями яв!
ляются векторы жордановой цепочки
11 2
10 2
10 2
X1(0) 5 30 4 , X2 (0) 5 31 4 , X3 (0) 5 30 4 .
30 4
30 4
31 4
6 7
6 7
6 7
Действительно, эти векторы удовлетворяют равенствам
1 A 3 4E 2 X3 (0) 5 X2 (0); 1 A 3 4E 2 X2 (0) 5 X1(0); 1 A 3 4E 2 X1(0) 5 0,
10 1 0 2
30 0 1 4 , т. е. представляют собой жорданову
A
E
5
6
7
где l = – а;
30 0 0 4
8
9
T
цепочку с порождающим вектором X3 (0) = [0 0 1] .
42
Найдя выходной сигнал для каждого из трех начальных условий
1
x1 1 e 1at, x2 1 te 1at, x3 1 t2e 1at,
2
убеждаемся, что в каждом случае он содержит только одно слагаемое.
Таким образом, использование векторов из жордановой цепочки в
качестве начальных условий позволяет выделять одни компоненты
свободного движения и подавлять другие.
Если в качестве порождающего вектора взять X3 (0) = [1 1 1]T, то
картина несколько изменится:
11 2
X2 (0) 5 (A 6 7E)X3 (0) 5 31 4 ,
30 4
8 9
11 2
X1(0) 5 (A 6 7E)X2 (0) 5 30 4 .
30 4
8 9
При этом выходные сигналы для начальных условий Х1(0), Х2(0),
Х3(0) будут такими:
x1 1 e 1at; x2 1 e 1at 2 te 1at; x3 1 e 1at 2 te 1at 2 0,5 t2e 1at.
Видно, что в первом из них отсутствуют две модальные компонен!
ты, а во втором – одна.
В принципе в качестве порождающего вектора жордановой цепоч!
ки можно брать любой корневой вектор высоты три – им является
всякий вектор с ненулевой третьей компонентой. Однако только век!
тор X(0) = [0 0 1]T удовлетворяет условию задачи. Он порождает це!
почку ортогональных векторов, имеющих особенно простой вид и
полезные свойства.
Определение жордановой цепочки векторов
по начальным условиям
Для того чтобы применять формулу (2.6) при решении диффе!
ренциальных уравнений, нужно уметь определять корневые век!
торы Н0 , ..., Нk–1, которые зависят от начальных условий. Выше
было показано, что для систем с единственным собственным чис!
лом порождающий корневой вектор равен вектору начальных ус!
ловий Н0 = Х0, а остальные векторы Hi находятся по формулам (2.7).
В случае, когда у системы есть и простые, и кратные собственные
числа, процедура определения жордановой цепочки векторов, входя!
щих в формулу (2.6), усложняется. Теперь вектор Н0 и другие векторы
цепочки будут зависеть не только от начальных условий системы (2.4),
но и от всей совокупности собственных векторов матрицы А.
Опишем эту процедуру подробнее, приняв для определенности, что
у матрицы А есть собственное число l0 кратности k и n – k простых
43
собственных чисел l1 , ..., ln–k . Тогда формула для общего решения
системы (2.4) получается объединением формул (2.5) и (2.6):
X(t) 3 c1G1e11t 4 1 4 cn 2kGn 2ke1n1kt 4
5
tk 21 6 10t
4 7 H0 4 H1t 4 1 4 Hk 21
8e .
7
1 k 9 12 ! 8
(2.10)
Здесь Gi – собственные векторы, принадлежащие простым соб!
ственным числам; Hi – корневые векторы, принадлежащие собствен!
ному числу l0, причем Нk–1 – единственный собственный вектор, при!
надлежащий этому собственному числу.
Постоянные ci и векторы Hi зависят от начальных условий системы
и подлежат определению. Основой для этого служит соотношение
X(0) 1 c1G1 2 1 2 cn 1kGn 1k 2 H0,
(2.11)
получаемое из формулы (2.10) при t = 0.
Изложим процедуру отыскания решения (2.10) системы диффе!
ренциальных уравнений (2.4) в виде следующего алгоритма.
Построение жордановой цепочки векторов
для решения системы дифференциальных уравнений
Ш а г 1. Определяются собственные векторы G1, 1, Gn 1k матри!
цы А, принадлежащие простым собственным числам.
Ш а г 2. Определяются постоянные c1, 1, cn 1k . Для этого обе
части соотношения (2.11) умножают на матрицу A k0 , где A 0 1 A 2 30E
и решают полученную систему уравнений:
A k0X(0) 3 A k0 1 c1G1 4 1 4 cn 1kGn 1k 2.
Вводя обозначения G 1 [G1, 1, Gn 1k ], C 1 [c1, 1, cn 1k ]T, ее можно
переписать в компактной форме A k0GC 1 A k0X0.
Число n уравнений в этой системе превышает число n – k неизвест!
ных, однако количество линейно независимых уравнений равно n – k,
поэтому система имеет решение, и притом только одно.
Ш а г 3. Из соотношения (2.11) определяется неизвестный век!
тор Н0:
H0 1 X(0) 2 GC.
Ш а г 4. На основе порождающего вектора Н0 строится жорданова
цепочка векторов
H1 1 A 0H0; H2 1 A 0H1, 1, Hk11 1 A 0Hk12.
44
Процедура отыскания решения системы уравнений (2.4) завер!
шается подстановкой найденных значений ci, Gi и Hi в формулу (2.10),
описывающую решение системы дифференциальных уравнений (2.4)
при наличии корня произвольной кратности.
Пример 3. Рассмотрим линейную систему, схема которой приве!
дена на рис. 2.4. Требуется установить вид сигналов x1, x2, x3 при
свободном движении системы:
а) из единичных начальных условий: a = b = c = 1;
б) из произвольных начальных условий: x1(0) = а, x2(0) = b, x3(0) = c.
c
b
1
p +a
x3
1
p +1
a
x2
1
p +1
x1
Рис. 2.4
Решение.
Ш а г 1. Описание системы в пространстве состояний характери!
зуется матрицей
2 11 1 0 3
A 6 4 0 11 1 5 .
4 0 0 12 5
7
8
Она имеет собственный вектор G1, отвечающий простому собствен!
ному числу l1 = – 2, и собственный вектор Н1, отвечающий кратному
собственному числу l = – 1:
112
G1 5 3 614 ,
314
7 8
11 2
H1 5 30 4 .
30 4
7 8
Следовательно, решение находим в форме
X(t) 3 c1G1e 12t 4 1 H0 4 H1t 2 e 1t ,
где постоянную с1 и векторы Н0 , Н1 считаем пока неизвестными.
Ш а г 2. При t = 0 имеем равенство
X(0) 1 c1G1 2 H0.
10 1 0 2
Умножая его на матрицу A 20, где A 0 5 A 6 E 5 30 0 1 4 ,
30 0 714
8
9
45
получаем уравнение
c1A 20G1 1 A 20X(0).
Для варианта а) оно имеет вид:
10 0 1 2 1 1 2 10 0 1 2 112
c1 30 0 514 3 514 6 30 0 514 314 ,
30 0 1 4 3 1 4 30 0 1 4 314
7
87 8 7
87 8
откуда с1 = 1.
Ш а г 3. Определяем вектор
112 1 1 2 10 2
H0 5 X(0) 6 c1G1 5 314 6 3 614 5 32 4 .
3 4 3 4 3 4
718 7 1 8 70 8
10 2
Ш а г 4. Определяем вектор H1 5 A 0H0 5 32 4 . Видим, что он с точ!
30 4
6 7
ностью до постоянного коэффициента совпадает с собственным век!
тором, определенным на первом шаге.
Подставляя найденные значения в исходную формулу, получаем
искомое решение
112
1t 2
X(t) 5 3 614 e 12t 7 2 31 4 e 1t ,
314
30 4
8 9
8 9
или в скалярной записи
x1 1 e 12t 2 2te 1t,
x2 1 3e 12t 2 2te 1t,
x3 1 e 12t.
В случае б) аналогичные выкладки приводят к следующему ре!
зультату:
1 5 a 3 c 6 5b 4 c 6 2
516
X(t) 7 9 b 4 c 4 9 0 t 8 e 1t 4 c 9 31
e 12t .
91
77 9 0 9 0 88
Его скалярная запись имеет вид:
x1 1 (a 2 c 3 (b 3 c)t)e 1t 3 ce 12t ,
46
x2 1 (b 3 c)e 1t 2 ce 12t , x3 1 ce 12t .
2.4. Разностные уравнения и цепочки корневых векторов
Разностные уравнения с кратными корнями
Перейдем к решению системы разностных уравнений
X(t 1 1) 2 AX(t),
X(0) 2 X0 ,
или в другой системе обозначений:
Xk 11 1 AXk, X0 ,
(2.12)
где А – постоянная матрица размера n´n; X – вектор переменных; t –
дискретное время.
Наиболее часто используется показательная форма записи реше!
ния Xk 1 A kX0, однако она не всегда удобна, например, в тех случа!
ях, когда надо получить аналитические формулы для отдельных ком!
понент вектора Х.
Приведем другую форму записи решения, аналогичную формуле
(2.5) для системы дифференциальных уравнений:
Xk 1 c1H121k 3 1 3 cnHn2nk ,
(2.13)
где li , Hi – собственные числа и собственные векторы матрицы А; ci –
произвольные коэффициенты.
Эта формула справедлива только для случая простых собствен!
ных чисел и задача состоит в том, чтобы обобщить ее на случай крат!
ных собственных чисел.
Пусть матрица А имеет собственное число l0 кратности m. Будем
искать соответствующие ему компоненты решения в виде
Xk 1 (H0 2 H1k 2 1 2 Hm 11km 11)3k0.
(2.14)
После подстановки этого выражения в уравнение (2.12) и сокра!
щения на множитель lk, получим
1
2
3H 5 H1(k 5 1) 5 1 5 Hm 11(k 5 1)m 11 4 6 0 7 A H0 5 Hkk 5 1 5 Hm 11km 11 .
8 0
9
Раскроем скобки и приравняем коэффициенты при одинаковых
степенях k справа и слева:
AHm 11 1 20Hm 11,
m 11
AHm 12 1 20Hm 12 3 20 (m 4 1)Hm 11, 1, AH0 1 20 5 Hi .
(2.15)
i 20
Вводя обозначение A 0 1 A 2 30E, перепишем эти равенства в форме
47
A 0Hm 11 1 0,
A20Hm 12 1 0, 1, A m
0 H0 1 0.
Следовательно, Hm 11, 1, H0 – корневые векторы матрицы А
высоты 1, 2, ..., m соответственно. Обозначив H = [H0, ..., Hm–1],
запишем равенства (2.15) в матричном виде
AH 1 20HB,
где
11 0
31 1
3
B 5 31 2
31 1
1
36 1 Cm
0
0
1
1
2
Cm
1
0
02
1
0
04
1
0
0 4.
4
1 1 14
m 12
1 Cm
1 47
11
(2.16)
Треугольная матрица В имеет размер m´m, в ее строках стоят би!
номиальные коэффициенты. Вычитая из обеих частей равенства
(2.16) матрицу l0Н, получим эквивалентное матричное равенство
A 0H 1 20HB0,
где
10 0
31 0
3
B0 5 B 6 E 5 3 1 2
31 1
1
37 1 Cm
0
0
0
1
2
Cm
1
0
02
1
0
04
1
0
0 4.
4
1 1 14
m 12
1 Cm
0 48
11
(2.17)
Присоединив к нему равенство Н0 = Х0, которое следует из выра!
жения (2.14) при t = 0, получим систему линейных алгебраических
уравнений относительно неизвестных векторов H1, 1, Hm 11.
Опишем процедуру их определения, не требующую обращения
матриц.
Многократно умножая обе части равенства (2.17) слева на матри!
цу А0, получим совокупность соотношений
11
m 11
m 11
A 0H 1 20HB0, A20H 1 220HB20, 1, A m
0 H 1 20 H0B0 . (2.18)
11
m 11
Обозначим через h0 1 H0, h1 1 A 0H0 / 20, 1, hm 11 1 A m
0 H0 / 2 0
промасштабированные векторы жордановой цепочки, построенные
на основе порождающего вектора Н0 , а через l = [0 1 ... 1]T – первый
столбец матрицы В0.
48
Приравняем первые столбцы в каждом из матричных равенств
(2.18):
h1 1 Hl, h2 1 HB0l, 1, hm 11 1 HB0m 11l,
или в матричной форме
h = HB1,
11
где h 1 [h1, 1, hm 11], H 1 [H0, 1, Hm 11], B1 1 [l, B0l, 1, Bm
0 l].
Отбросив в матрице В1 первую строку (она нулевая), а в матри!
це Н – первый столбец (вектор Н0), получим систему уравнений
[h1, 1, hm 11 ] 1 [H1, 1, Hm 11]B2,
(2.19)
где В2 – нижнетреугольная матрица размера (m – 1) ´ (m – 1).
Из системы (2.19) легко последовательно определить все векторы
Hi, начиная с вектора Hm–1 (он с точностью до масштабного множи!
теля будет совпадать с собственным вектором матрицы А). Тем са!
мым векторы Hi окажутся выраженными через векторы hi с помощью
формулы вида
[H1, 1, Hm 11] 1 [h1, 1, hm 11]D.
(2.20)
Здесь D – нижнетреугольная матрица, ее элементы постоянны и
не зависят ни от матрицы А, ни от начальных условий системы.
Приведем матрицы D для кратности собственных чисел от m = 2 до
m = 5:
1
31
3 1
D4 5 3 6
3 2
31
73 3
11
D2 5 1, D3 5 3 1
36
7 2
02
1 4,
4
28
1 1
2
3 1
04
35
4
3 2
0 4 , D5 6 3 1
4
3 3
14
3 1
6 84
37 5 4
0
1
2
1
5
2
11
24
0
1
2
1
6
2
0
0
1
6
1
5
4
0 2
4
0 4
4
.
0 4
4
1 4
24 48
Видно, что каждая следующая матрица получается окаймлением
предыдущей. Элементы первого столбца матрицы D описываются
формулой di1 1 21/ i , а диагональные элементы – формулой dii 1 1/ i !
Это позволяет выписать в явном виде выражения для векторов Н1 и
Hm–1 при произвольном m:
49
(11)m 11
1
1
H1 2 h1 1 h2 3 h3 1 1 1
hm 11;
2
3
m 11
1
1
hm 11; hk 2 k A 0kH0.
Hm 11 2
(m 1 1)!
40
(2.21)
Пример 1. Найдем решение системы разностных уравнений
12 3 3 2
Xk11 5 30 2 34 Xk;
30 0 2 4
6
7
18 2
X0 5 38 4 .
38 4
6 7
В данном случае l1 = l2 = l3 = 2, m = 3, поэтому решение будет
иметь вид
Xk 1 (H0 2 H1k 2 H2k2 )2k.
Находим векторы Н0, h1, h2:
1152
1182
18 2
10 3 3 2
1
1
H0 5 X0 5 38 4 , A 0 5 30 0 34 , h1 5 A 0H0 3124 , h2 5 A 0h1 3 0 4 .
304
304
38 4
30 0 0 4
2
2
6 7
6 7
6 7
6
7
Векторы Н1 и Н2 вычисляем по формулам (2.21):
1152
1
H1 5 h1 6 h2 5 3124 ,
3 4
2
708
192
1
H2 5 h2 5 30 4 .
30 4
2
7 8
Следовательно, решение системы разностных уравнений имеет вид:
18 5 15k 5 9k2 2
3 18 2
4
1152
19 2
2
k
Xk 8 68 7 5 6127 k 5 60 7 k 9 2 6 8 5 12k 7 2k.
6
7
607
60 7
88 68 7
99
8
6
7
Определение цепочки корневых векторов
по начальным условиям
Корневые векторы Н1, ..., Hm–1, входящие в формулу (2.20), од!
нозначно определяются вектором Н0, порождающим жорданову це!
почку h1, ..., hm–1. Вектор Н0, в свою очередь, зависит от начальных
условий системы. Если у матрицы А системы разностных уравнений
(2.20) нет других собственных чисел, кроме l0, то эта зависимость
крайне проста: Н0 = Х0.
50
В более общем случае, когда у матрицы А кроме собственного чис!
ла l0 кратности m есть простые собственные числа l1, ..., ln–m, эта
зависимость становится более сложной. Формула для решения сис!
темы (2.12) принимает вид
1
2
Xk 3 c1G141k 5 1 5 cn 1mGn 1m4 nk 1m 5 H0 5 H1k 5 1 5 Hm 11km 11 4 k0, (2.22)
где Gi – собственные векторы матрицы А, отвечающие простым соб!
ственным числам; Hi – корневые векторы матрицы А, принадлежа!
щие собственному числу l0.
Постоянные сi и векторы Нi зависят от начальных условий и под!
лежат определению.
Соответствующая процедура, как и в случае дифференциальных
уравнений, включает выписывание равенства
X0 1 c1G1 2 1 2 cn 1mGn 1m 2 H0,
и исключение вектора Н0 путем умножения на матрицу
(2.23)
Am
0
3 1 A 4 5 0E 2 :
m
Am
0 X0 1 A 0 (c1G1 2 1 2 cn 1 mGn 1m ).
m
(2.24)
Из этой системы линейных уравнений определяются коэффици!
енты сi, после чего из соотношения (2.23) находится вектор Н0. Да!
лее на его основе строится жорданова цепочка векторов h1, ..., hm–1 и
из системы уравнений (2.19) находят векторы Нi.
Опишем эту процедуру в виде следующего алгоритма.
Определение корневых векторов
для решения системы разностных уравнений
Ш а г 1. Определяются собственные векторы G1, ..., Gn–m матри!
цы А, принадлежащие ее простым собственным числам.
Ш а г 2. Из системы (2.24) определяются постоянные c1, ..., cn–m.
Вводя обозначения G = [G1, ..., Gn–m], C = [c1, ..., cn–m], ее можно
записать в более компактной форме
m
Am
0 GC 1 A 0 X0.
Число n уравнений в этой системе на m превышает число неизвес!
тных из!за того, что часть уравнений является линейно зависимой (и
при решении может не учитываться).
Ш а г 3. Из соотношения (2.23) определяется вектор Н0:
H0 1 X(0) 2 GC.
Ш а г 4. На основе порождающего вектора Н0 строится промасш!
табированная жорданова цепочка векторов
51
h0H0,
!1
m 11
h1 1 A 0H0/20, ..., hm11 1 A m
0 H0 / 2 0 .
Ш а г 5. Из уравнения (2.19) или (2.20) определяются корневые
векторы H1, ..., Hm–1.
Процедура завершается подстановкой найденных значений ci, Gi
и Hi в формулу (2.22), описывающую решение системы разностных
уравнений (2.12) при наличии корня произвольной кратности.
Пример 2. (Динамика популяции рыб). Рассмотрим популяцию
рыб, живущих в водоеме, выделив три возрастные группы – однолеток
(мальков), рыб среднего и старшего возрастов. Заданы величины р1, р2
– вероятности дожития особями каждой возрастной группы до следую!
щего возраста, и числа а1, а2, а3, характеризующие среднюю плодови!
тость каждой возрастной группы.
Требуется выяснить, как изменяется со временем численность воз!
растных групп и каково будет их процентное соотношение через дос!
таточно большое время.
Р е ш е н и е . Обозначим начальную численность возрастных групп
через х1(0), х2(0), х3(0), а их численность спустя k лет – через х1(k),
х2(k), х3(k). Тогда динамику популяции можно описать системой трех
разностных уравнений:
x1(k 1 1) 2 a1x1(k) 1 a2x2 (k) 1 a3x3 (k),
x2 (k 1 1) 2 p1x1(k),
x3 (k 1 1) 2 p2x2 (k).
Переходя к матричной форме записи, получаем
X(k 5 1) 6 AX(k),
1 a1
A 6 3 p1
3
70
a2
0
p2
a3 2
0 4.
4
08
Примем для определенности а1 = 0, а2 = 9, а3 = 12, р1 = 1/3, р2 = 1/2.
Здесь компоненты вектора Х отражают возрастную структуру по!
пуляции, в которой выделено три возрастных группы. Числа 9 и 12
характеризуют среднюю плодовитость поколений (младшая группа
не дает приплода), величины 1/3 и 1/2 – доля рыб, выживающих в
младшей и средней группе и переходящих в следующую возрастную
группу.
Пусть в начальный момент имеется одна рыба старшего возраста,
т. е. X(0) 1 [0 0 1]T. Проводя прямые итерационные вычисления
на основе формулы:
52
9 122
1 0
Xk11 5 31/3 0
04 X ,
3 0 1/2 0 4 k
6
7
10 2
X0 5 30 4 ,
31 4
6 7
k 5 0, 1, 2, 1,
для первых восьми лет получаем следующий массив значений:
10 12 0 36 24 108 144 3722
X 5 30 0 4 0 12 8
36 48 4 .
31 0 0 2 0
6
4
18 47
6
Первая строка характеризует динамику изменения численнос!
ти мальков, вторая и третья – рыб среднего и старшего возрастов.
Из приведенного массива трудно усмотреть закон изменения числен!
ности популяции.
Чтобы установить аналитический вид этого закона, воспользуем!
ся приведенным выше алгоритмом.
Ш а г 1. Для определения собственных чисел матрицы А выпи!
шем ее характеристическое уравнение
1z
9 12
A 1 2E 3 1/3 1z 0 3 23 1 32 1 2 3 (2 1 2)(2 4 1)2 3 0.
0 1/2 1z
Матрица А системы имеет собственное число l1 = 2, которому от!
вечает собственный вектор G1 = [2 4 4 1]T, и собственное число l0 = – 1
кратности два. Следовательно, решение должно иметь вид:
Xk 11 3 c1G12k 4 1 H0 4 H1k 21 512 .
k
Ш а г 2. Постоянную с 1 находим из соотношения (2.24)
A 20X0 1 c1A 20G1 . Подставляя в него значения
9 122
1 1
A 0 5 31 3 1
04 ,
3 0 12 14
6
7
124 2
A 20X0 5 3 4 4 ,
314
6 7
124 2
A 20G1 5 9 3 4 4 ,
314
6 7
получаем с1 = 1/9.
Ш а г 3. Вектор Н0 определяем из соотношения (2.23):
20 3
224 3
2 124 3
1
1
H0 6 X0 1 c1G1 6 40 5 1 4 4 5 6 4 14 5 .
41 5 9 4 1 5 9 4 8 5
7 8
7 8
7
8
53
Ш а г 4. Определяем вектор h1:
2 1123
1
h1 6 1 A 0H0 6 4 4 5 .
3 4 12 5
7
8
Ш а г 5. Вектор Н1 находим из системы (2.20). Учитывая, что в
нашем случае D2 = 1, получаем Н1 = h1.
Следовательно, решение имеет вид:
3 6 56 7
624 7
6 56 7 4
k
1
4
2
Xk 4 2k 8 53 2 k 9 1 512 ,
8
9
9 1
3 51
89 2
9
или в скалярной форме записи
1
1
1
2
k
1
24 4 2k 5 6 1 8 6 6k 21 612 ,
9
k
1
x2 (k) 3 4 4 2k 6 2 1 8 6 6k 21 612 ,
9
k
1
x3 (k) 3 2k 5 1 8 6 6k 21 612 .
9
x1(k) 3
2
2
Fishes population
i
=
2
xi
1 100 %
x1 + x2 + x3
100
21
80
60
40
20
22
23
0
4
Время
8
Рис. 2.5
54
n
12
16
Очевидно, что с течением времени доминирующую роль будут иг!
рать экспоненциально возрастающие слагаемые, поэтому процент!
ный состав популяции будет определяться компонентами первого
собственного вектора.
Графики относительной численности рыб каждого возраста, оп!
ределяемые формулами 1i 2
xi
3 100%, i 2 1, 2, 3, приведены
x1 4 x2 4 x3
на рис. 2.5. Их установившиеся значения относятся, как 24:4:1.
Во многих практических приложениях желательно поддерживать
численность популяции на постоянном уровне (количество рыб в ак!
вариуме или пруду, численность деревьев в лесу). Полученное реше!
ние позволяет ответить на вопрос о том, какую часть популяции надо
для этого регулярно изымать, т. е. как установить научнообоснован!
ные квоты на вылов рыбы, нормы вырубки леса и т. п.
Ответ прост: для поддержания популяции на постоянном уровне
главное собственное число l матрицы А должно быть равно 1, если
же оно больше единицы, то доля изымаемой части должна состав!
лять 100 (l – 1)/l%. Например, при l = 2 следует регулярно изымать
половину популяции, а при l = 1,1 – около 9%.
2.5. Задачи и упражнения
1. Построить жорданову цепочку векторов для матриц
13 1 0
39 34 2
15
40 2 1
111 4 2
A1 6 4
, A 2 6 4 6 311 35 5 , A 3 6 4
5
4
3
0 0 2
3
4
5
7
8
4
7 37 13 6 8
47 31 31 31
02
05
,
15
5
1 58
взяв в качестве порождающего вектора b вектор с единичными ком!
понентами.
У к а з а н и е . Собственные числа матрицы А1: l1 = l2 = 7; матри!
цы А2: l1 = l2 = l3 = 0; матрицы А3: l1 = l2 = l3 = l4 = 2.
2. Построить корневые подпространства для матриц
1 31 1 1 2
A1 6 4 33 2 25 ,
4 31 1 1 5
7
8
11 1 02
A 2 6 4 34 32 1 5 ,
4 4 1 325
7
8
12
4 31
A3 6 4
0
4
47 1
3 0
0 1
3 2
2 31
32
25
.
35
5
0 58
О т в е т . А1: базис корневого подпространства для l = 0 – вектор
[0 1 – 1]T, для l = 1 – векторы [0 1 1]T, [0 1 0]T.
55
А2: базис корневого подпространства для l = –1 – векторы [1 0 0]T,
[0 1 0]T, [0 0 1]T.
А3: базис корневого подпространства для l = 2 – векторы [2 – 1 0 0]T,
[1 0 1 0]T, [2 0 0 1]T.
А4: базис корневого подпространства для l = –2 – вектор [0 1 0 –1].
3. Обобщить формулу (2.10) на случай двух кратных корней.
Как в этом случае получить систему линейных уравнений для оп!
ределения постоянных сi?
1 1 AX с мат!
4. Решить систему дифференциальных уравнений X
рицей
22 13
41 12
A64
0 0
4
470 0
4
2
2
1
16 3
14 5
135
5
1258
и найти начальные условия, при которых переменная х1 содержит
только одну из модальных компонент.
У к а з а н и е . Базис корневых подпространств для l = –1 – векто!
ры [1 1 0 0]T, [0 0 1 1]T. Базис корневых подпространств для l = 1 –
векторы [3 1 0 0]T, [0 –2 3 1]T.
5. Если порядок минимального полинома матрицы А меньше по!
рядка ее характеристического полинома, то один из корней l1, ..., ln–m
может совпадать с корнем l0. Как будет выглядеть формула (2.10) в
этом случае?
1 1 AX с мат!
6. Решить систему дифференциальных уравнений X
рицей
14 1 12
A 5 3 62 1 624 ,
31 1 44
7
8
у которой имеются две жордановы клетки, отвечающие собственно!
му числу l0 = 3. Принять Х0 = [a b c]T.
7. Обобщить формулу (2.12) на случай двух кратных корней. Как
в этом случае получить систему линейных уравнений для определе!
ния постоянных сi?
8. Решить систему разностных уравнений Xk+1 = AXk с матрицей
21 1 12
42 1 0
A64
1 0 1
4
740 11 2
56
03
25
15 .
5
1 85
У к а з а н и е . Все собственные числа матрицы А равны единице,
им отвечают две жордановы клетки второго порядка.
9. Если порядок минимального полинома матрицы А меньше по!
рядка ее характеристического полинома, то один из корней l1, ..., ln–m
может совпадать с корнем l0,. Как будет выглядеть формула (2.24) в
этом случае?
10. Решить систему разностных уравнений Xk+1 = AXk с матри!
цей
11 1 1
3 51 3 0
A63
51 0 51
3
73 0 51 51
02
14
14 .
4
1 84
У к а з а н и е . Все собственные числа матрицы А равны единице,
им отвечают две жордановы клетки третьего и первого порядка.
57
3. РАЗНОСТНЫЕ УРАВНЕНИЯ И ZAПРЕОБРАЗОВАНИЕ
3.1. Определение zAпреобразования
При исследовании непрерывных динамических систем с дискрет!
ным временем (они называются также импульсными системами) и
при решении разностных уравнений широко применяется математи!
ческий аппарат z!преобразования (другие названия: дискретное пре!
образование Лапласа, преобразование Лорана).
О п р е д е л е н и е . Пусть f(t) – функция целочисленного аргумента
t = 0, 1, 2 ... .
Ее z!преобразованием называется функция комплексного аргумен!
та z, определяемая формулой
F(z) 1 f (0) 2
1
f (1) f (2)
2 2 2 1 1 3 f (t)z2t .
z
z
t 30
(3.1)
Условная запись z!преобразования:
f (t) 1 F(z) .
Свойства z(преобразования
1. Линейность af1(t) 1 bf2 (t) 2 aF1(z) 1 bF2 (z) .
2. Сдвиг по времени
f (t 1 1) 2 z( F(z) 3 f0 ) f (t 1 2) 2 z2 ( F (z) 3 f0 3 z11f1) ,
f (t 1 n) 2 zn ( F(z) 3 f0 3 z11f1 3 ... 3z1(n 11)fn 11) .
3. Умножение на время
tf (t) 1 2z
dF (z)
.
dz
Приведем примеры z!преобразований простейших функций:
f(t) =
®F(z) = 1.
58
1 1
z
f (t) 1 1 2 F(z) 1 1 3 3 2 3 1 1
,
z z
z 41
2 3
z
f (t) 1 t 2 F(z) 1 1 3 3 2 3 1 1
,
z z
(z 4 1)2
2
e1 5 e1 6
z
f (t) 1 e 2 F(z) 1 1 3
,
37 8 311
z 9 z z 4 e1
z(z 4 cos )
z sin cos t 2 2
; sin t 2 2
.
z 4 2z cos 3 1
z 4 2z cos 3 1
1t
Таблица zAпреобразований
Оригинал f(t)
F(z)
Единичный импульс
1
1
(–1)t
t
t2
at
eat
sinw t
cosw t
shw t
chw t
atsinw t
atcosw t
t at
f(t+1)
f(t+2)
z
z 11
z
z 11
z
(z 1 1)2
z(z 1 1)
(z 2 1)3
z
z1a
z
z 1 ea
z sin 1
z2 2 2z cos 1 3 1
z(z 1 cos 2)
z2 1 2z cos 2 3 1
zsh1
z2 2 2zch1 3 1
z(z 1 ch2)
z2 1 2zch2 3 1
az sin 1
z2 2 2az cos 1 3 a2
z(z 1 a cos 2)
z2 1 2az cos 2 3 a2
az
(z 1 a)2
z (F(z)–f0)
z2(F(z)–f0–z–1f1)
59
При выводе последних формул удобно использовать теорему о диф!
ференцировании изображений (свойство 3). Она доказывается через
дифференцирование соотношения (3.1) по переменной z.
Заметим, что в тулбоксе Symbolic пакета MATLAB имеются ко!
манды ztrans и iztrans для выполнения прямого и обратного z–преоб!
разований. Например, в результате выполнения команд syms t,
ztrans(t^2) получим ans = z*(z+1)/(z–1)^3.
3.2. Решение разностных уравнений
с помощью zAпреобразования
Пусть дано разностное уравнение n!го порядка
y(t + n) + an–1y(t+n–1) +...+ a0y(t) = f(t)
с начальными условиями y(0) = y0; y(1) = y1; ...; y(n–1) = yn–1.
Алгоритм решения разностного уравнения
с помощью z(преобразования
Ш а г 1. Применить z!преобразование к исходному разностному
уравнению, заменяя y(t) на Y(z); y(t+1) на z(Y(z)–y0) и т. д.
Ш а г 2. Из полученного алгебраического уравнения выразить Y(z).
Ш а г 3. Выполнить разложение на простые дроби.
Ш а г 4. Пользуясь таблицей, выполнить обратное z!преобразование.
Рассмотрим примеры применения этого алгоритма.
Пример 1. Требуется решить разностное уравнение первого по!
рядка
(3.2)
y(t+1)–0,9y(t) = 0,1; y(0) = 5.
Ш а г 1. Применяем z!преобразование:
0,1z
z(Y(z)–5)–0,9Y(z) =
.
z 11
Ш а г 2. Выражаем Y(z):
5z 1 0,1z /(z 2 1)
.
z 2 0,9
Ш а г 3. Раскладываем на простые дроби
Y(z) 3
1 A
5(z 3 1) 4 0,1
B 2
5 z6
4
7,
(z 3 0,9)(z 3 1)
8 z 3 0,9 z 3 1 9
5z 3 4,9 5 Az 3 A 4 Bz 3 0,9B,
5 5 A 4 B; 3 4,9 5 3 A 3 0,9B; B 5 1; A 5 4,
z
z
4
Y(z) 5 4
.
z 3 0,9 z 3 1
Y(z) 5 z
60
Ш а г 4. Выполняем обратное преобразование
y(t) = 4(0,9)t+1.
(3.3)
Проверка для t = 0, 1, 2 по уравнениям (3.2) и (3.3) дает одинако!
вые результаты
y(0) = 5; y(1) = 4,6; y(2) = 4,24.
5
График полученного решения пока!
зан на рис. 3.1.
4
Пример 2. Пусть дано разностное
3
уравнение второго порядка
2
yn+2 +a1 yn+1+ a0 yn = un
с заданными начальными значениями
1
y0 и y1.
0
0
5
10 15 20 25 30
Ш а г 1. Применяем к нему z!пре!
Рис. 3.1
образование:
z2 3 Y(z) 5 y0 5 y1z11 4 6 a1z1 Y(z) 5 y0 2 6 a0Y(z) 7 U(z).
8
9
Ш а г 2. Выражаем Y(z):
Y(z) 2
y1z 1 y0z(z 1 a1 ) 1 U
z2 1 a1z 1 a0
.
Ш а г 3. Чтобы найти yn, надо разложить правую часть на про!
стые дроби.
Положим а1 = –1; а0 = 0,5; un = 0; y0 = 0; y1 = 1. Тогда получаем
Y(z) 1
z
.
z 2 z 3 0,5
2
Корни знаменателя – комплексные, следовательно, эта дробь –
простая.
Ш а г 4 . С помощью таблицы преобразований находим оригинал
yn = can sinwn. Постоянные с, а, w определяем из соотношений:
a2 1 0,5; 2a cos 2 1 1; c 1
1
,
a sin 2
откуда
a2
2
1
3 0,707; 4 2 ; c 2 2.
2
4
61
Окончательный вид решения:
yn 1 2 2 21n 2 sin(3 n /4).
Пример 3. (Задача о банковском вкладе). Мы открываем счет в бан!
ке, наш начальный взнос составляет В долларов. Каждый год добав!
ляем по b долларов, годовой процент равен k. Какая сумма окажется
на нашем счете через n лет? Какова будет прибыль?
Обозначая искомую величину через yn, можно записать следую!
щее разностное уравнение:
yn+1 = yn(1 + k) + b; y0 = B.
Заметим, что этому уравнению можно поставить в соответствие схе!
му моделирования, содержащую один элемент задержки ЭЗ (рис. 3.2, а)
а)
б)
B
b
yn+1
800
yn
ЭЗ
400
k+1
A
0
0
2
4
6
8
10
Рис. 3.2
Найдем решение двумя способами – с помощью характеристичес!
кого уравнения и с помощью z!преобразования.
П е р в ы й с п о с о б . Находим корень характеристического уравнения
z–(1+k) = 0; z1 = 1+k.
Решение находим в виде
y = yодн+ yчаст,
где yодн = C(1 + k)n; yчаст = –b/k.
Следовательно, yn = C(1 + k)n – b/k.
Неизвестный коэффициент С находим из начальных условий
y0 = B = C – b/k , C = B + b/k.
Таким образом, окончательно имеем
yn = (B + (b/k))(1 + k)n – (b/k) = B(1 + k)n + (b/k)((1 + k)n– 1).
Примерный вид графика решения показан на рис. 3.2, б.
Суммарная прибыль за n лет составит zn = yn – xn, где xn – общий
взнос. Она определяется формулой
zn = yn – (B + nb).
62
Пусть B = 100$, b = 50$, k = 5%, n = 10 лет. Подставим эти данные
в полученное решение:
1
2
yn 3 1100 4 1,05n 5 1000 3 1000 1,1 4 1,05n 5 1 .
При n = 10 получаем: y10 = 791,78, z10 = y10 – 600 = 191,78.
Таким образом, всего за 10 лет мы внесли 600$, получили в итоге
792$, прибыль составила 192$.
В т о р о й с п о с о б (z!преобразование).
Исходное разностное уравнение имеет вид yn+1 = (k + 1)yn + b; y0 = B.
Применяя к нему z!преобразование
yn 1 Y(z); yn 11 1 z(Y(z) 2 B),
получаем алгебраическое уравнение:
z(Y(z) 1 B) 2 (k 3 1)Y(z) 3 b
z
.
z 11
Выражаем из него Y(z):
Y(z) 5
z
z 4 1 5 z 2 B(z 4 1) 1 b 3 5 z Bz 4 B 1 b .
6 (z 4 k 4 1)(z 4 1) 7
z 4 k 41
(z 4 k 4 1)(z 4 1)
8
9
Bz 1 b
Раскладываем это выражение на две простые дроби
b
z
b z
Y(z) 3 16 B 4 27
5
.
k 9 z 5 k 51 k z 51
8
Возвращаемся к оригиналам
n
b
b
yn 5 38 B 6 49 1 k 6 12 7 .
k
k
Ответ получили тот же, что и первым способом.
Пример 4. Решим уравнение второго порядка с вещественными
корнями
x(t 1 2) 2 5x(t 1 1) 1 6x(t) 3 0; x0 3 1; x1 3 2.
Находим изображение
X(z) 3
z(z 1 5 2 2)
z(z 1 3)
z
3
3
;
z2 1 5z 2 6 (z 1 2)(z 1 3) z 1 2
x(t) 3 2t.
В данном случае имеет место сокращение нуля и полюса.
Рассмотрим еще несколько примеров.
63
Пример 5. Решим уравнение третьего порядка с кратными корнями
yn+3 – 5yn+2 + 8yn+1 – 4yn = 0; y0 = 0, y1 = 2, y2 = 1.
Находим изображение и выполняем разложение на простые дроби
Y(z) 3
z(2(z 1 5) 2 1)
2z2 1 9z
17z 7z
5z
3
3
2
1
.
2
2
z
1
1
z
1
2
z 1 5z 2 8z 1 4 (z 1 1)(z 1 2)
(z 1 2)2
3
Выполняем обратное z!преобразование
5
yn 1 27 3 7 4 2n 2 n 2n.
2
Пример 6. Найдем решение задачи об электрической цепи (см.
рис. 1.1) с помощью z!преобразования. Примем для определенности
R = 2r, u0 = 1. Эта цепь описывается разностным уравнением (1.7),
которое при R = 2r принимает вид
uk12 1 4uk 11 2 uk 3 0.
Применяем к нему z!преобразование:
z2U 1 zu0 1 zu1 1 3(zU 1 zu0 ) 2 U 3 0,
(z2 1 4z 2 1)U 3 z(u1 1 4u0 2 zu0 ).
Отсюда, учитывая что u0 = 1, находим
U(z) 3
z(z 1 4 2 u1 )
.
z2 1 4z 2 1
Обратное z!преобразование выполняем с помощью таблиц, для
чего последнюю формулу удобнее представить в виде суммы двух сла!
гаемых:
U(z) 2
z(u 1 2)
z(z 1 2)
3 2 1
.
z 1 4z 3 1 z 1 4z 3 1
2
Переходя к оригиналам, получаем формулу для расчета напряже!
ния в произвольной точке R!цепи:
uk = ch kw + C sh kw .
(3.4)
Константа w находится из условия ch w = 2, откуда w @ 1,32; С –
постоянная, зависящая от u1.
Проблема заключается в том, что значение u1 неизвестно и для
определения постоянной С нужно дополнительное условие. Обойти
эту проблему удается, рассмотрев последнее звено схемы рис. 1.1 (при
64
R = 2r оно представляет собой делитель напряжения с коэффициен!
том 1/3) и выписав краевое условие un = un–1 /3.
С учетом этого можно записать
3(ch nw + C sh nw) = ch (n – 1)w + C sh (n – 1)w,
откуда находим
C3
3ch n1 2 ch (n 2 1)1
.
sh (n 2 1)1 2 3sh n1
Воспользуемся тождествами гиперболической тригонометрии:
sh (x – y) = sh x ch y – ch x sh y; ch (x – y) = ch x ch y – sh x sh y.
Раскрывая с помощью этих формул ch (n – 1)w и sh (n – 1)w в числи!
теле и знаменателе и учитывая, что ch w = 2, sh1 2 3, получаем
C34
3shn1 2 chn1
.
shn1 2 3chn1
Подставляя это выражение в формулу (3.4) и выполняя преобра!
зования, получаем окончательный результат:
uk 4
sh(n 1 k)2 3 3ch(n 1 k)2
.
shn2 3 3chn2
В частности, для k = 0 (вход схемы) имеем u0 = 1, а для k = n (выход
схемы) получим
un 1
1
3
chn2 3
shn2
3
.
Если схема содержит только одно звено (простой делитель напря!
жения), то
1
u1 1
ch2 3
3
sh2
3
1
1 .
3
3.3. Цифровые интеграторы
и дискретные передаточные функции
Для описания цифровых линейных систем используются разно!
стные уравнения различных порядков и дискретные передаточные
функции Q(z).
65
О п р е д е л е н и е . Дискретной передаточной функцией линейной
системы с входом x(t) и выходом y(t) называется отношение zAпреоб!
разования выходного сигнала к zAпреобразованию входного сигнала
при нулевых начальных условиях Q(z) = Y(z)/X(z).
Аналоговый интегратор описывается дифференциальным уравне!
нием первого порядка y1 1 x и передаточной функцией Q(p) = 1/p. Про!
стейший цифровой интегратор описывается разностным уравнением
первого порядка, которое реализует численное интегрирование по
методу Эйлера (методу прямоугольников). Несколько большую точ!
ность обеспечивает цифровой интегратор, использующий формулу
трапеций.
Найдем дискретную передаточную функцию цифрового интегра!
тора для двух вариантов его реализации.
Метод Эйлера. Согласно этому методу интегрирование прибли!
женно выполняется с помощью схемы, приведенной на рис. 3.3.
x
Т
ЭЗ
y
Рис. 3.3
Первый блок схемы осуществляет усиление входного сигнала в Т
раз, второй – задержку на время Т, где Т – шаг дискретного времени t
= nT, n = 0, 1, 2, ... . Схеме соответствует разностное уравнение:
y((n+1)T)–y(nT) = Тx(nT) , или yn+1 – yn = Т xn.
Применяя z!преобразование, находим дискретную передаточную
функцию
Q(z) = Y(z) / X(z) = T / (z–1).
Ее можно было получить из непрерывной передаточной функции
интегратора Q(p) = 1/p с помощью замены переменных
p = (z–1) / T.
(3.5)
Вывод тех же формул через дифференциальное уравнение ин!
тегратора основан на замене производной отношением приращений:
y1 1 x;
(yn 11 2 yn )/T 1 xn .
Метод трапеций. Согласно этому методу дискретная передаточ!
ная функция интегратора выбирается в виде Q(z) 2 T 3 z 1 1 , что соот!
2 z 41
2 z 11
ветствует замене p 2 3
.
T z 41
66
Это стандартная формула замены при переходе от Q(p) к Q(z) по мето!
ду трапеций. При этом обе передаточные функции Q(p) и Q(z) будут дроб!
но!рациональными передаточными функциями одного порядка.
Операторный вывод этой формулы основан на приближенном со!
отношении:
z 2 e pT 2
e
e
pT
2
1
pT
2
pT
2 2 2 1 pT ,
3
pT 2 4 pT
14
2
11
где e–pT – известная из теории преобразования Лапласа передаточная
функция блока запаздывания на время Т.
Выразим отсюда р через z: 2+pT = 2z – pTz; Tp(1 + z) = 2(z – 1),
откуда
p2
2 z 11
.
T z 31
(3.6)
Другой вывод основан на подстановке в дифференциальное урав!
y 1y
x 2 xn
, в резуль!
нение интегратора y1 1 x значений y1 3 n 11 n ; x 4 n 11
T
2
тате чего получаем разностное уравнение
yn 11 1 yn 2
T
(xn 11 2 xn ).
2
На рис. 3.4, а в качестве иллюстрации показан результат интегриро!
вания синусоидального сигнала с помощью трех цифровых интеграто!
ров, использующих методы Эйлера (кусочно!постоянный сигнал), тра!
пеций (кусочно!линейный сигнал) и Рунге!Кутты (непрерывная линия).
Графики получены путем моделирования в пакете SIMULINK с
помощью схемы, приведенной на рис. 3.4, б (она соответствует зна!
чению Т = 1).
Формулы (3.5) и (3.6) позволяют осуществлять приближенный
переход от непрерывной передаточной функции любого объекта Q(p)
к его дискретной передаточной функции Q(z). Шаг квантования вре!
мени Т должен выбираться с учетом динамики объекта и быть по край!
ней мере на порядок меньше его постоянных времени или периода
собственных колебаний.
Пример 7. Рассмотрим колебательное звено с передаточной функ!
цией
Q( p) 1
1
.
p 2 k2
2
67
а) 7
6
5
4
3
2
1
0
–1
0
10
20
30
40
б)
1
?
Integrator
Scope4
1
z
Sine Wave
Unit Delay
0.5
Gain
1
z
Unit Delay
Рис. 3.4
Его импульсная весовая функция имеет вид синусоиды q(t) = sin kt
с периодом 2p/k. Найдем дискретную передаточную функцию этого
колебательного звена.
Шаг Т выберем таким образом, чтобы на период собственных ко!
лебаний приходилось 20 точек дискретного времени: T = 0,1p/k. Ис!
пользуя, согласно формуле (3.5), замену p = (z–1)/T, получаем
68
Q1(z) 1
T2
T2
1 2
.
2
2 2
(z 2 1) 3 k T
z 2 2z 3 k2T 2 3 1
Применение формулы трапеций (3.6) дает более точный резуль!
тат:
Q2 (z) 2
T 2 (z 1 1)2
T 2z2 1 2T 2z 1 T 2
2 2 2
.
2
2 2
2
4(z 3 1) 1 k T (z 1 1)
(k T 1 4)z2 3 (8 3 2k2T 2 )z 1 k2T 2 1 4
Учитывая, что kТ = 0,314; к2 Т2 = 0,0987, получаем
Q1(z) 2
1
0,0987
0,07665 z2 1 2z 1 1
;
Q
(
z
)
2
.
2
k2 z2 3 2z 1 1,0987
k2
z2 3 1,928z 1 1
В первом случае свободный член в знаменателе немного больше
единицы, следовательно, полюсы лежат вне единичной окружности
и решение (весовая функция дискретной системы) будет медленно
расходиться.
Во втором случае полюсы передаточной функции лежат точно на
единичной окружности в плоскости z, поэтому решение будет носить
чисто колебательный характер.
Рассмотрим случай матричного описания дискретной системы
Xk11 1 AXk 2 BUk,
Yk 1 CXk 2 DUk,
(3.7)
где Xk 1 R n – вектор состояния; Uk, Yk 1 R – входной и выходной сиг!
налы, A, B, C, D – матрицы соответствующих размеров.
Для того чтобы найти дискретную передаточную функцию такой
системы, применяем z!преобразование:
zX(z) 1 AX(z) 2 BU(z);
Y (z) 1 CX(z) 2 DU(z).
Выразим из первого уравнения X(z) и, подставляя его во второе:
Y (z) 1 [C(zE 2 A) 11 B 3 D]U(z).
Следовательно, дискретная передаточная функция системы име!
ет вид
Q(z) 1 C(zE 2 A) 11 B 3 D.
(3.8)
Для систем с одним входом и одним выходом она представляет
собой отношение двух полиномов
Q(z) 2
bn 11zn 11 1 1 1b1z 1 b0
,
zn 1 an 11zn 11 1 1 1a1z 1 a0
(3.9)
69
причем знаменатель совпадает с характеристическим полиномом
матрицы А.
Формула (3.8) справедлива как для скалярных систем, так и для
систем с несколькими входами и выходами. В общем случае С и В –
прямоугольные матрицы, Q(z) – матричная передаточная функция,
ее элементы – рациональные дроби типа (3.9).
Пример 8. Дискретная система задана уравнениями
yk 11 1 xk,
xk 11 1 25xk 2 6yk 3 7uk.
Найти передаточную функцию от входа u до выхода у и импульс!
ную весовую функцию.
В данном случае имеем
1y2
X 3 4 5,
7x 8
10 12
A34
,
7 66 6585
10 2
B 3 4 5,
77 8
C 3 [0
1].
Выписываем матрицу zE – А и преобразуем ее
5 z 31 6
zE 3 A 7 8
,
6 z 4 59
1 zE 3 A 211 7
1
5 z 4 5 16
.
(z 4 2)(z 4 3) 8
36 z 9
Дискретную передаточную функцию получаем по формуле (3.8):
Q(z) 1
7z
.
(z 2 2)(z 2 3)
Импульсную весовую функцию находим с помощью обратного z!
преобразования
qk 3 7[1 42 2 4 1 43 2 ].
k
k
Все эти выкладки можно сделать в тулбоксе Symbolic пакета
MATLAB. Соответствующая программа имеет вид:
syms z;
A = [0 1;6 5];B = [0; 7];C = [0 1];
Q = C*inv(z*eye(2)А))*В;
q = iztrans(Q);
Q,q
Q = 7*z/z^2+5*z+6),q = 7*(2)^n7*(3)^n
70
% описание в пространстве состояний
% передаточная функция
% весовая функция
% результат
3.4. Задачи и упражнения
1. Найти функцию f(k), если ее z!преобразование имеет вид
F(z) 2
2z2 1 (e3 1 z)z
.
z 1 (2 3 e3 )z 1 2e3
2
Р е ш е н и е . Выполняем разложение на простые дроби
F(z) 1
z
z
2
.
z 3 2 z 3 e3
Переходя к оригиналам, получаем ответ:
f (k) 1 2k 2 e3k.
2. Решить разностное уравнение
yk12 1 2eh yk 11 2 e2h yk 3 0, y0 1 0, y1 1 heh ,
используя zAпреобразование.
3. Найти передаточную функцию и реакцию дискретной системы
yk+2 + yk+1 – 6yk = 3uk+1 – 3uk
с нулевыми начальными условиями на линейно возрастающее вход!
ное воздействие uk = k при k ³ 0 (uk = 0 при k < 0).
Р е ш е н и е . Применяем z!преобразование:
Y (z) 2
3(z 1 1)
3z 1 3
U(z) 2
U(z) 2 Q(z)U(z),
(z 3 3)(z 1 2)
z 3z16
2
где Q(z) – дискретная передаточная функция.
Так как U(z) 1
z
3z
, то Y (z) 1
.
(z 2 1)2
(z 2 3)(z 3 2)(z 3 1)
Далее следует выполнить разложение на простые дроби и перейти
к оригиналам.
1
3
3
(23)k 11 3 2k 2 .
20
5
4
4. Цепная схема представляет собой соединение n идентичных
Т!образных резистивных звеньев (рис. 3.5). Найти ее выходное на!
пряжение un, если u0 = 1В и все резисторы одинаковы.
У к а з а н и е . Эта задача отличается от примера 6 только краевым
условием un = un–1 /2, откуда C = –thnw.
О т в е т . yk 1 2
71
u0
u1
un–1
un
...
Рис. 3.5
О т в е т . Формула для выходного напряжения имеет вид
un 1
1
,
chn2
ch2 1 2.
5. Найти дискретную передаточную функцию системы, заданной
матрицами описания в пространстве состояний:
10 12
A34
,
61 175
112
B 3 4 5,
617
C 3 [0
1],
d 3 1.
О т в е т . Дискретная передаточная функция имеет вид
Q(z) 1
72
z2
.
z2 2 z 2 1
4. СТРУКТУРНЫЕ МОДЕЛИ РАЗНОСТНЫХ УРАВНЕНИЙ
4.1. Модели уравнений с начальными условиями
При компьютерном моделировании разностных уравнений можно
использовать численные и структурные методы. Программа вычис!
лений в первом случае представляет собой простой цикл, обеспечива!
ющий расчет очередного значения решения на основе уже рассчитан!
ных предыдущих. Такой способ удобен для решения разностных урав!
нений с заданными начальными условиями, но непригоден для ре!
шения краевых задач, таких, как рассмотренная ранее задача об элек!
трической цепи.
Альтернатива состоит в использовании структурного моделиро!
вания, когда разностному уравнению сопоставляется некоторая
структурная схема на сумматорах, элементах задержки и других бло!
ках. Возможности для такого моделирования имеются в пакетах
VISSIM, SIMULINK и ряде других. Остановимся подробнее на двух
подходах к структурному моделированию разностных уравнений.
Согласно первому хорошо известному подходу схема моделирова!
ния составляется так же, как это делается при моделировании диф!
ференциальных уравнений – методом понижения производной (ме!
тодом Кельвина) [9]. Если, например, дано разностное уравнение вто!
рого порядка
x(t 1 2) 1 a1x(t 1 1) 1 a0x(t) 2 bu(t),
(4.1)
то его сначала разрешают относительно “старшего” члена, в нашем
случае это x(t + 2):
x(t 1 2) 2 3a1x(t 1 1) 3 a0x(t) 1 bu(t).
(4.2)
Затем, предполагая, что член x(t + 2) известен, пропускают со!
ответствующий сигнал через два элемента задержки ЭЗ1 и ЭЗ2
(рис. 4.1), каждый из которых задерживает свой входной сигнал на
один такт. Выходные сигналы элементов задержки складывают в со!
ответствии с уравнением (4.2) и получают сигнал x(t + 2), подавая
который на вход первого элемента задержки, замыкают схему моде!
лирования.
73
b
u(t)
–a 0
1
x(t+2)
x(t)
x(t+1)
ЭЗ1
ЭЗ 2
–a 1
Рис. 4.1
Для разностного уравнения n!го порядка потребуется схема, со!
держащая n элементов задержки. Схема моделирования, построен!
ная таким образом, подобна схеме моделирования дифференциаль!
ных уравнений и отличается только использованием элементов за!
держки вместо интеграторов. Это отличие облегчает программную
реализацию схем, но приводит к серьезным практическим трудно!
стям при их аппаратной реализации ввиду отсутствия точных и не!
дорогих элементов со временем задержки порядка секунды.
Второй подход удобен тем, что не требует использования элемен!
тов задержки, а может быть реализован с помощью усилителей и сум!
маторов. Поясним его сущность на примере моделирования однород!
ного уравнения n!го порядка
x(t 1 n) 1 an 11x(t 1 n 2 1) 1 ... 1 a1x(t 1 1) 1 a0 x(t) 3 0 .
(4.3)
Вводя индексацию переменных, перепишем его в форме
xk 1n 3 4 1 an 21xk 1n 21 5 ... 5 a0xk 2 .
(4.4)
Если нам известны начальные значения функции x(t)
x(0) 1 x0,
x(1) 1 x1, ..., x(n 2 1) 1 xn 11 ,
то, подставляя эти значения в формулу (4.4) при k = 0, можно найти
x(n). Далее, полагая в формуле (4.4) k = 1 и зная x(1), ..., x(n), можно
найти x(n+1); полагая k = 2, можно найти x(n + 2) и т. д.
Таким образом, используя однотипные операции суммирования,
можно шаг за шагом вычислить любое требуемое число точек решения.
Этот метод решения разностных уравнений получил название «метода
шагов». Соответствующая схема моделирования показана на рис. 4.2.
x0
x1
x1 a0
x2 a
1
a0
a1
...
1
–
xn
...
xn–1
a n–1
a n–1
1
–
x2
x3
a0
a1
xn+ 1 ...
a n–1
Рис. 4.2
74
1
–
xk
xk+1
xn+2
...
a0
a1
...
a n–1
1
–
xn+k
Число сумматоров в ней определяется требуемым числом точек
решения (числом шагов), а число входов у каждого сумматора равно
порядку разностного уравнения. Решение снимается путем последо!
вательного подключения измерительного прибора к выходам сумма!
торов (вручную или с помощью специального коммутатора).
Для иллюстрации моделирования разностных уравнений на анало!
говых сумматорах обратимся к примерам, приведенным в разд. 1.
Геометрическая прогрессия 1, a, a2, a3, ... (см. разд. 1, пример 2)
описывается разностным уравнением первого порядка xk+1 = axk , x0 = 1.
Оно моделируется с помощью цепочки последовательно соединен!
ных усилителей (рис. 4.3).
x0
x1
a
x2
a
x3
a
xk
a
...
Рис. 4.3
Арифметическая прогрессия, уравнение электрической цепи и ряд
Фибоначчи (см. разд. 1, примеры 3, 4, 5 ) описываются однородным
разностным уравнением второго порядка вида
xk 12 1 a1xk 11 1 a0xk 2 0 .
(4.5)
Схема моделирования такого уравнения при заданных начальных
условиях x0, x1 приведена на рис. 4.4.
...
x1
a1
–1
x0
a0
x2
a0
a1
–1
x3
a1
a0
x4
–1
a0
...
a1
–1
xk
Рис. 4.4
4.2. Пример моделирования в пакете SIMULINK
Проиллюстрируем методику моделирования методом шагов на
примере разностного уравнения второго порядка
(4.6)
1,9x(t + 2) – 0,81x(t + 1)+ 0,9x(t) = 0, x0 = 1, x1 = –1.
Требуется получить его решение на интервале 0 £ t £ 12.
Переписываем исходное уравнение в виде
x(t + 2) = 1,426x(t + 1) – 0,474x(t).
75
Используя библиотеку блоков SIMULINK, создадим схему, пока!
занную на рис. 4.5, и оформим ее в виде отдельного блока с помощью
процедуры маскирования. Его входные сигналы x(t) и x(t + 1), вы!
ходные сигналы x(t + 1) и x(t + 2).
Gain1
In1
1
–0.474
2
1.426
In2
Out1
1
Sum
2
Out2
Gain2
Рис. 4.5
Чтобы получить решение в 12 точках, построим схему из 12 таких
блоков, соединенных последовательно (рис. 4.6). Блок “To
Workspace” обеспечивает передачу результатов в рабочую область
MATLAB для их последующей обработки.
На выходах схемы получаем искомые значения: x = [ 1; –1; –0,9;
0,09; 0,891; 0,7209; –0,1531; –0,7866; –0,5702; 0,1948; 0,6885;
0,4443; –0,2197],
X0
1
–1
X1 Блок1
Блок2
X12
X2
Mux
To Workspace
x
Рис. 4.6
График полученного решения приведен на рис. 4.7. Для его пост!
роения использовались команды MATLAB:
t = 0:12; tt = 0:.25:12; xx = interp1(t,x,tt,’spline’);
plot(tt,xx),grid; plot(t,x,t,x,’o’,tt,xx),grid
76
1
1
0,5
0
– 0,5
–1
–1,5
0
4
8
12
Рис. 4.7
Аналогичным образом выполняется структурное моделирование
нелинейных разностных уравнений и краевых задач.
4.3. Модели уравнений с краевыми условиями
Аналитическое решение краевых разностных уравнений практи!
чески ничем не отличается от решения уравнений с заданной началь!
ной последовательностью. Однако с точки зрения компьютерного
моделирования или численного решения эти задачи существенно раз!
личаются. В частности, для краевых задач не удается непосредствен!
но применить описанный выше метод шагов (подумайте, почему?) и,
следовательно, связанные с ним итерационный алгоритм решения и
схему моделирования по методу Кельвина (см. рис. 4.1).
Трудности, возникающие при решении краевых задач, преодолева!
ют, сводя решение уравнения с краевыми условиями к многократному
решению того же уравнения с разными начальными последовательнос!
тями, подобно тому как это делается при решении краевых дифферен!
циальных уравнений. Использование схемы моделирования на анало!
говых сумматорах (см. рис. 4.2) существенно упрощает процесс нахож!
дения решения, поскольку мы получаем все точки решения одновре!
менно. Например, решение краевой задачи методом пристрелки сводит!
ся к тому, чтобы, изменяя значения начальной последовательности,
добиться выполнения краевых условий в заданных точках.
77
Проиллюстрируем это на примере.
Пусть требуется найти решение уравнения (4.5) с краевыми усло!
виями
x(0) = x0, x(k) = xk.
Схема моделирования этого уравнения показана на рис. 4.4. Для
получения искомого решения достаточно на вход первого сумматора
подать известную величину x0 и плавно изменять величину x1 до тех
пор, пока напряжение на выходе схемы не станет равным заданной ве!
личине xk . На этом моделирование задачи завершается.
Менее эффективен описанный метод при решении задач высокого
порядка, когда неизвестны не один, а несколько членов начальной пос!
ледовательности. В таких случаях приходится одновременно варьиро!
вать несколько входных сигналов схемы моделирования, добиваясь
того, чтобы выходные сигналы нескольких сумматоров (по количеству
краевых условий) приняли заданные значения.
Причина затруднений состоит в том, что схема моделирования не
вполне адекватна моделируемой задаче, поскольку мы строим схему для
решения уравнения с заданной начальной последовательностью, а ре!
шаем уравнение с краевыми условиями. Проблема построения алгорит!
ма или устройства, наиболее полно отражающего структуру задачи –
одна из основных проблем вычислительной техники, и решить ее удает!
ся далеко не всегда.
В рассматриваемом случае решение дает метод, который можно на!
звать методом краевой структурной схемы: согласно этому методу
структурная схема моделирования должна отражать не только само
уравнение, но и вид краевых условий.
Структурная схема (см. рис. 4.2) была построена по уравнению (4.3),
разрешенному относительно старшего слагаемого. Это позволило легко
учитывать заданную начальную последовательность. Разрешая то же
уравнение относительно других слагаемых, будем получать другие ва!
рианты структурной схемы, позволяющие учитывать различные виды
краевых условий.
Например, уравнение (4.5) можно разрешить относительно среднего
члена
xk 11 1 2xk 3 4xk 12, 2 1 5
a0
1
, 415
a1
a1
и использовать это выражение для составления схемы моделирования.
Мы по!прежнему получим цепочку последовательно включенных сум!
маторов, однако теперь каждая точка решения формируется не на осно!
ве двух предыдущих, как это было на рис. 4.4, а на основе двух соседних
78
точек – предыдущей и последующей (рис. 4.8). Благодаря этому в схеме
моделирования появились контуры обратных связей, обеспечивающие дви!
жение информации как слева направо (учет начальной точки x0), так и
справа налево (учет конечной точки xk). При этом отпадает необходимость
в каком бы то ни было подборе начальной последовательности – решение
формируется сразу с учетом всех краевых условий.
x0
2
1
1
x1
1
2
1
x2
2
1
1
x3
1
x4
1
1 1
2
1
xk–2
2
1
1
xk –1
xk
Рис. 4.8
Увеличение порядка моделируемого уравнения не вызывает прин!
ципиальных трудностей. Пусть требуется решить уравнение (4.3) с
краевыми условиями, часть из которых задана в начале интервала
решения (n1), а часть – в конце (n2), причем n1 + n2 = n:
x(1) 1 u1, ..., x(n1) 1 un1 , x(k 2 n1) 1 v1, ..., x(k) 1 vn2 .
Разрешим уравнение (4.3) относительно члена x(t + n1) и по полу!
ченному выражению составим схему моделирования. Она будет отли!
чаться от схемы, показанной на рис. 4.3, тем, что у каждого сумматора
будет n1 входов, соединенных с предшествующими сумматорами и n2
входов, соединенных с последующими сумматорами. У крайнего левого
сумматора на первую группу входов будут поданы напряжения
u1,..., un1 , а у крайнего правого сумматора на вторую группу входов бу!
дут поданы напряжения v1, ..., vn2 . Благодаря наличию в схеме пря!
мых и обратных связей на выходах всех сумматоров мгновенно устано!
вится решение, удовлетворяющее краевым условиям.
Аналогичным образом может быть выполнено структурное мо!
делирование нестационарных и нелинейных разностных уравне!
ний. В случае нестационарных уравнений коэффициенты a, b на
входах разных сумматоров в схеме рис. 4.8 будут иметь разные
значения, а в случае нелинейных уравнений в схеме появятся, кроме
усилителей и сумматоров, нелинейные блоки (логарифмические, три!
гонометрические, произведения, деления и др.).
4.4. Матричное решение разностных уравнений
Структурные модели, описанные выше, позволяют получать чис!
ленное решение разностных уравнений с начальными и краевыми
79
условиями. Рассмотрим еще один способ решения подобных задач.
Он опирается на матричное представление линейных разностных
уравнений. Изложим его суть на примере линейного неоднородного
уравнения второго порядка
a2xk 12 1 a1xk11 1 a0xk 2 fk
(4.7)
с краевыми условиями х0 = а, хn+1 = b.
Для того чтобы найти решение этого уравнения на интервале
0 1 k 1 n 2 1, выпишем систему линейных алгебраических уравнений
a2x2 1 a1x1 1 a0x0 2 f0,
a2x3 1 a1x2 1 a0x1 2 f1,
1111111111
a2xn 1 a1xn 11 1 a0xn 12 2 fn 12,
a2xn 21 1 a1xn 1 a0xn 11 2 fn 11.
(4.8)
Перейдем к матричной форме записи:
AX 5 F,
1 a1 a2 0
3a a a
1
2
3 0
0 a0 a1
3
A5
31 1 1
30 0 0
30 0 0
63
1 0 0 02
1 0 0 04
4
1 0 0 04
,
1 1 1 14
1 a0 a1 a2 4
1 0 a0 a1 447
1 x1 2
3 x 4
3 2 4
X 5 3 2 4,
3xn 11 4
63 xn 47
2 f0 1 a0a 3
4
5
f1
4
5
F64 1
5.
4 fn 12 5
47 fn 11 1 a2b 58
Матрица этой системы трехдиагональная, причем на диагоналях
стоят коэффициенты разностного уравнения. Решение дается фор!
мулой X 1 A 11F, компьютерная реализация которой не вызывает
трудностей.
Изложенный способ можно применять для решения линейных
разностных уравнений произвольного порядка, в том числе неста!
ционарных. При этом матрица А останется ленточной (ширина лен!
ты определяется порядком разностного уравнения). Тип краевых ус!
ловий также может быть любым, он влияет лишь на вид вектора F и
80
матрицы А (так, в случае начальных условий она будет треуголь!
ной).
Можно сказать, что структурные схемы, приведенные на рис. 4.1,
4.4, 4.8, отвечали различным численным алгоритмам решения систе!
мы (4.8). Например, схема рис. 4.4 реализует рекуррентный способ ре!
шения этой системы при известных значениях х0, х1, а схема рис. 4.8 –
процедуру решения краевой задачи методом простой итерации. На мат!
ричном языке итерационный процесс описывается формулой
Xn 11 3 1 E 4 A / a1 2 Xn 5 F / a1,
которая опирается на соотношение a1X 1 a1X 2 AX 3 F. Для сходи!
мости процесса нужно, чтобы спектральный радиус (или какая!ни!
будь из норм) матрицы E 1 A / a1 была меньше единицы. В частности,
при использовании первой (строчной) нормы получаем следующее
достаточное условие сходимости:
a0 1 a2 2 a1 .
Пример. Решим матричным способом задачу Фибоначчи о кроли!
ках. Она сводится к разностному уравнению (1.8)
xk 12 1 xk11 1 xk 2 0, x0 2 x1 2 1.
Его матричное представление имеет вид АХ = F, где
X 1 [x2, x3, ..., x12 ]T, А – треугольная матрица 11!го порядка с еди!
ницами на главной диагонали и с элементами –1 на двух нижележа!
щих диагоналях. Два первые элемента вектора F равны 2 и 1, а ос!
тальные равны нулю.
Сформируем указанные матрицы в пакете MATLAB:
F=[2 1 0 0 0 0 0 0 0 0 0]’; c=[1 1 1 0 0 0 0 0 0 0 0]’;
r=[1 0 0 0 0 0 0 0 0 0 0]; A=toeplitz(c,r);
11 0 0 0 0 0 0 0 0 0
3 51 1 0 0 0 0 0 0 0 0
3 51 51 1 0 0 0 0 0 0 0
3
3 0 51 51 1 0 0 0 0 0 0
3 0 0 51 51 1 0 0 0 0 0
A 6 3 0 0 0 51 51 1 0 0 0 0
3 0 0 0 0 51 51 1 0 0 0
3
3 0 0 0 0 0 51 51 1 0 0
3 0 0 0 0 0 0 51 51 1 0
3 0 0 0 0 0 0 0 51 5 1 1
3
73 0 0 0 0 0 0 0 0 51 51
02
04
04
4
04
04
04 .
0 44
04
04
04
4
1 84
81
Решение получаем с помощью команды X=inv(A)*F:
X’=2 3 5 8 13 21 34 55 89 144 233.
Последняя компонента вектора Х указывает количество кроли!
ков в конце года х12 = 233.
4.5. Задачи и упражнения
1. Составьте блок!схему программы решения разностного уравне!
ния методом шагов. Сравните ее со схемой на рис. 4.2.
2. Каким образом можно “смоделировать” решение разностного
уравнения, используя всего один решающий усилитель? Что для этого
нужно?
3. Составьте схему моделирования уравнения Фибоначчи.
4. Сформулируйте краевые условия в задаче об электрической цепи
(подразд. 1.1, пример 4 ).
5. В чем состоит метод краевой структурной схемы? Составьте с
помощью этого метода схему моделирования электрической цепи (под!
разд. 1.1, пример 4, подразд. 1.2, пример 2, подразд. 3.2, пример 6).
6. Составьте три варианта структурных схем для моделирования
уравнения (4.5).
7. Найдите аналитическое решение разностного уравнения (4.6) и
сравните его с графиком (см. рис. 4.7), полученным путем моделиро!
вания в SIMULINK.
8. Решите матричным способом задачу об электрической цепи
(рис. 1.1), рассмотрев режимы холостого хода и короткого замы!
кания.
9. Выполните задание п. 8 для электрической цепи, приведенной
на рис. 3.5.
82
5. НЕЛИНЕЙНЫЕ РАЗНОСТНЫЕ УРАВНЕНИЯ
5.1. Источники нелинейных разностных уравнений
Нелинейные разностные уравнения часто встречаются при ис!
следовании и моделировании различных прикладных задач в тех!
нике, экономике, социологии, демографии и других дисциплинах.
В математике нелинейные разностные уравнения появляются, в
частности, при анализе сходимости различных итерационных про!
цессов.
Отметим две формы записи нелинейных разностных уравнений: в
виде одного уравнения порядка k
(5.1)
F(xn+k, ..., xn) = 0
и в виде системы уравнений первого порядка (описание в простран!
стве состояний)
(5.2)
Xn+1 = Ф(Xn),
k
где X 1 R – вектор состояния; Ф – нелинейная векторная функция.
Общей теории исследования нелинейных разностных уравнений
нет. При их изучении стараются ответить на качественные вопросы,
такие как отыскание стационарных точек, анализ устойчивости, пе!
риодичности, наличия предельных циклов и аттракторов.
В последние десятилетия интересные результаты были достигну!
ты в ходе исследования очень простых математических моделей, опи!
сываемых разностными уравнениями первого порядка:
xn 11 1 f (xn , r ), x1 1 a, n 1 1, 2, 3, ...,
(5.3)
где r – некоторый параметр.
Уравнение (5.3) можно рассматривать как явную формулу, кото!
рая позволяет по числу xn находить следующее число xn+1 и таким
образом определяет всю последовательность {xn}. Можно сказать, что
соотношение (5.3) соответствует некоторому итерационному процес!
су. Поэтому последовательность чисел {xn} часто называют итераци!
ями начальной точки x1 1 a.
Простейшие уравнения вида (5.3) были рассмотрены в разд. 1.
При f (xn , r ) 1 rxn соотношение (5.3) определяет геометрическую про!
83
грессию со знаменателем r. При f (xn , r ) 1 r 2 xn – арифметическую
прогрессию.
Однако если f(x, r) простейшая нелинейная функция аргумента х,
например, квадратичная парабола, то свойства последовательности
{xn} могут оказаться совершенно необычными.
Приведем, следуя книге [6], несколько типичных ситуаций, в ко!
торых уравнения вида (5.3) выступают как математические модели
либо возникают при численном решении уравнений.
Моделирование процессов с дискретным временем
Допустим, исследуется динамика численности популяции неко!
торых животных и существует возможность оценивать эту числен!
ность только раз в год. Тогда естественно временную переменную счи!
тать дискретной, которая может принимать только целые действи!
тельные значения п = 1, 2, ... . Численность популяции будет выра!
жаться функцией дискретного аргумента x(n) или, как ее часто обо!
значают, xn. Последовательность x1, x2, ... будем обозначать {xn}.
Можно считать, что xn – численность популяции в год с номером п.
Ясно, что численность популяции в данный год xn+1 зависит от того,
сколько животных было год назад, т. е. от величины xn. В простейшем
случае, когда величина xn+1 зависит только от численности в предыду!
щий год, мы приходим к математической модели вида (5.3). В этой мо!
дели f – непрерывная однозначная функция своих аргументов, r – пара!
метр, который зависит от того, какую конкретную задачу мы рассмат!
риваем, а – начальное значение, первый член в последовательности {xn}.
Часто используется функция f вида rx(N–х):
xn+1 = rxn(N–хn), 0 1 xn 1 N.
(5.4)
Формула (5.4) показывает, что если rN > 1, то численность вида
растет, пока она мала xn 11 N. В силу ограниченности ресурсов числен!
ность животных начинает убывать, когда животных становится слиш!
ком много. Это учитывается с помощью ограничивающего квадра!
тичного члена. Удобно сделать замену переменных xn 1 yn N, r 1 q / N.
При этом формула (5.4) приобретает вид
yn+1 = qyn(N–yn), 0 1 yn 1 N.
(5.5)
Отображения (5.4) и (5.5) часто называют логистическим отобра!
жением.
Нас интересует вопрос о том, что произойдет с различными вида!
ми по прошествии достаточно долгого времени, т. е. каковы аттрак!
торы изучаемого отображения. Для ответа на него в этой простейшей
84
модели достаточно выяснить, какой будет последовательность
{xn}, n 1 2 при различных значениях r.
Дискретизация дифференциальных уравнений
Простейшая модель изменения численности популяции, предло!
женная Мальтусом более 200 лет назад, описывается линейным диф!
ференциальным уравнением первого порядка x1 1 ax. В соответствии
с ней численность популяции должна неограниченно расти. Для того
чтобы учесть ограниченность ресурсов, доступных популяции, и внут!
ривидовой отбор, естественно ввести нелинейный ограничивающий
член. Это приводит к уравнению x1 1 ax(1 2 x).
Единицы измерения выбраны так, чтобы предельная численность
популяции равнялась 1. Если 0 < х(0) < 1, то х(t) ® 1 при t ® ¥. В
этом нетрудно убедиться, проводя качественное исследование этой
модели. Неподвижная точка х = 0 является неустойчивой, точка
х = 1 устойчивой.
Применим к этому уравнению метод Эйлера, заменив производ!
ную приращением функции
x((k 1 1)h) 2 x(kh) 3 ahx(kh)(1 2 x(kh)),
где h – шаг дискретного времени.
Обозначая x(kh) 1 xk, получаем
xk11 1 2xk 3 4xk2, 2 1 ah 5 1, 4 1 ah.
Мы вновь получили одномерное разностное уравнение вида (5.4).
Анализ показывает, что его свойства похожи на свойства дифферен!
циального уравнения только при небольших h. Когда шаг по време!
ни h превышает некоторое критическое значение, поведение этих двух
объектов начинает качественно отличаться.
Решение нелинейных алгебраических уравнений
Стандартным методом численного решения нелинейных алгебра!
ических уравнений вида F(x) = 0 является метод простой итерации.
В соответствии с ним для нахождения корней исходное алгебраи!
ческое уравнение переписывается в виде x = f(x), и затем ему сопос!
тавляется разностное уравнение
xn+1 = f(xn), n = 0, 1, 2, ... .
При определенных условиях решение этого уравнения образует пос!
ледовательности {xn}, сходящиеся к корню исходного уравнения. Со!
ответствующий итерационный процесс иллюстрируется рис. 5.1, а.
85
Другой стандартный метод численного решения алгебраического
уравнения вида x = f(x) – метод Ньютона (метод касательных). В нем
последовательность {xn} строится на основе формулы
xn 11 2 xn 3 f (xn )/ f 1(xn ) ,
т. е. для получения следующей точки проводится касательная к кри!
вой в предыдущей точке.
Графическое изображение итерационного процесса по методу Нью!
тона показано на рис. 5.1, б.
f(x)
f(x)
x3
x3
x2
x2
x3 x1
x
x4 x3
x2
x1
x
Рис 5.1
В обоих случаях возникают разностные уравнения первого поряд!
ка вида (5.3).
Обработка экспериментальных данных
Пусть некоторый сложный процесс, развивающийся во времени,
характеризуется функцией х(t). Для построения его математической
модели поступим следующим образом. Выделим локальные макси!
мумы функции х(t). Первый из них обозначим М1, второй – М2, k!й –
Мk. Первый максимум достигается в момент t1, второй – в момент t2, и
т. д. На плоскости {Мk, Мk+1} будем откладывать точки с координата!
ми (Мk, Мk+1), т. е., первая точка будет (М1, М2), вторая (М2, М3).
Оказывается, для некоторых колебательных химических реак!
ций, математических моделей гидродинамики и ряда других систем
точки {Мk, Мk+1} с высокой точностью ложатся на однозначные не!
прерывные кривые Мn+1 = f(Мn).
Наличие такой функции позволяет строить простые феноменоло!
гические модели изучаемых явлений. Они дают возможность по пре!
дыдущим значениям локальных максимумов предсказывать следу!
ющие, т. е. прогнозировать дальнейший ход процесса, исходя из его
86
предыстории. В обсуждаемом случае мы можем предсказать величи!
ну Мn+1, если известен максимум Мn и функция f.
Таким образом, нелинейные разностные уравнения представляет
собой удобный математический аппарат для описания многих про!
цессов и явлений.
5.2. Решение нелинейных разностных уравнений
Для нелинейных разностных уравнений, так же как для диффе!
ренциальных и алгебраических уравнений, не существует общих ме!
тодов решения, они разработаны только для специальных видов урав!
нений. Так, в некоторых случаях удается заменой переменных свес!
ти нелинейное разностное уравнение к одному или нескольким ли!
нейным.
Остановимся на одном из таких случаев, когда разностное уравне!
ние имеет вид
yk11 2
ayk 1 b
,
cyk 1 d
(5.6)
где a, b, c, d – постоянные коэффициенты.
Это разностное уравнение вида (5.3) с дробно!линейной функци!
ей f. Его нелинейность связана с наличием операции деления в пра!
вой части.
Уравнение (5.6) можно переписать в эквивалентной форме
dyk11 1 cyk11yk 2 ayk 3 b.
(5.7)
Теперь нелинейность связана с наличием члена, содержащего про!
изведение переменных.
Изложим два способа решения таких уравнений.
С п о с о б 1 (замена переменной).
Выполняя в уравнении (5.6) замену переменных по формуле
yk 1 2 3
1
,
zk
получим
12
1
zk11
1
)2b
1zk11 2 1 a(1 zk 2 1) 2 bzk
zk
3
.
3
,
1
zk11
c(1zk 2 1) 2 dzk
c(1 2 ) 2 d
zk
a(1 2
87
Избавимся от знаменателей, выполнив переход к уравнению типа
(5.7) и потребуем, чтобы коэффициент при нелинейном члене zk zk+1
равнялся нулю:
12c 2 1d 3 a1 2 b,
12c 2 1(d 4 a) 4 b 3 0.
Положим с = 1 (это всегда можно сделать, поделив числитель и
знаменатель правой части формулы (5.6) на с): 12 2 1(a 2 d) 2 b 3 0.
Если дискриминант этого квадратного уравнения положителен
D 1 (a 2 d)2 3 4b 4 0,
(a 2 d)2 4 24b,
то существует вещественное число a, для которого указанная замена
переменных приводит к линейному уравнению для z.
Найдя a и выполнив обратную замену переменных, получим реше!
ние исходного нелинейного разностного уравнения (5.6).
Пример 1. Требуется решить нелинейное разностное уравнение
yk11 2
2yk 1 1
.
yk
Здесь a = 2, b = –1, c = 1, d = 0. Выписываем квадратное уравнение
для определения a:
1 2 2 21 3 1 4 0.
Оно имеет единственный корень a = 1. Делаем замену переменных
yk 1 1 2
z
1
1
: 11
223 k .
zk
zk11
1 1 zk
Приводим к общему знаменателю
(zk11 1 1)(zk 1 1) 2 zk11(zk 1 2).
Раскрывая скобки, получаем линейное разностное уравнение
zk11 1 zk 2 1.
Его общее решение имеет вид zk 1 C 2 k, откуда находим решение
исходного нелинейного разностного уравнения
yk 1 1 2
1
,
C2k
где С – произвольная постоянная.
С п о с о б 2 (переход к системе линейных уравнений).
88
Представим переменную yk в виде отношения двух новых пере!
менных
yk 1
uk
.
vk
Тогда уравнение (5.6) принимает вид
uk11 auk 1 bvk
2
.
vk11 cuk 1 dvk
Приравнивая отдельно числители и знаменатели, получаем сис!
тему двух линейных разностных уравнений
uk 11 1 auk 2 bvk,
vk 11 1 cuk 2 dvk ,
или в матричной форме
Xk11 3 AXk,
1u 2
Xk 3 4 k 5 ,
6vk 7
1a b 2
A34
.
6 c d 75
Ее решение может быть записано с помощью одной из формул
X k 1 A kX 0 ,
Xk 1 c1H1z1k 2 c2H2z2k,
где X0 –вектор начальных условий; zi и Hi – собственные числа и
собственные векторы матрицы А.
Окончательное решение получаем в виде отношения компонент
вектора Хk.
Пример 2. Найдем указанным способом решение предыдущего
примера. Выписываем эквивалентную систему линейных разностных
уравнений
3uk 11 1 2uk 2 vk
uk 11 2uk 2 vk
1
4 5
vk11
uk
6vk11 1 uk.
Матричная запись этой системы:
2c 3
22 113
Xk11 4 5
Xk , X0 4 5 1 6 .
71 0 86
7 c2 8
Вычислим последовательные степени матрицы А:
33 12 4
3 4 13 4
3 5 14 4
3k 2 1 1k 4
A2 5 6
, A3 5 6
, A4 5 6
, ...., A k 5 6
.
7
7
7
2
1
1
3
1
2
4
1
3
1 1 k 79
8
9
8
9
8
9
8 k
89
Следовательно, решение линейной системы имеет вид
4uk 1 c1(k 2 1) 3 c2k,
5k 2 1 3k 6
Xk 1 8
X0; 7
9
1 3 k
k
vk 1 c1k 2 c2 (1 3 k).
Решение исходного уравнения получаем как отношение uk и vk:
yk 3
где c 1
c1(k 1 1) 2 c2k
c
311
,
c1k 2 c2 (k 2 1)
ck 1 1
c1
2 1 – произвольная постоянная, зависящая от начального
c2
условия.
Видим, что решения, полученные двумя способами, совпадают.
а)
u0
б)
u1
u2
un
...
u1
u0
1
1
R n–1
Рис. 5.2
Пример 3. Рассмотрим задачу об определении входного сопротив!
ления цепной схемы, содержащей n одинаковых Г!образных звеньев
(рис. 5.2, а), если все сопротивления в схеме имеют величину 1 Ом.
Непосредственное вычисление входного сопротивления для слу!
чаев n = 1, 2, 3, ... дает:
R1 = 2; R2 1 1 2
53
2 5
13
1 ; R3 1 1 2
1 ; ...
3 3
12 5 3 8
Мы видим, что в полученных дробях фигурируют отношения чи!
сел Фибоначчи 1, 1, 2, 3, 5, 8, 13, 21, ... . Можно предположить, что
справедлива общая формула Rn 2
12n
, где 1n – n!е число Фибонач!
12n 11
чи. Тогда в пределе при n ® ¥ получаем, что входное сопротивление
схемы стремится к золотому сечению
R1 2 3,
90
32
11 5
4 1,618.
2
Приведем аналитический вывод формулы для входного сопротив!
ления. Пусть найдено входное сопротивление Rn–1 схемы, содержа!
щей n–1 звено. Добавляя одно звено на входе, получим схему рис.
5.2, б. Находим ее входное сопротивление Rn:
Rn 1 1 2
Rn 11
.
Rn 11 2 1
Мы получили разностное уравнение с дробно!линейной функцией
в правой части.
Для его решения воспользуемся способом 2:
Rn 2
2Rn 11 1 1 un
2 .
Rn 11 1 1 vn
Отсюда получаем систему двух линейных разностных уравнений
3un 1 2un 11 2 vn 11,
4
5vn 1 un 11 2 vn 11.
Находим характеристическое уравнение этой системы, его корни
и собственные векторы матрицы А:
12 12
A36
;
81 179
z1,2 5
A 4 zE 3 z2 4 3z 5 1,
32 5
3 14
; H1 5 6 7 ;
2
91 314
H2 5 6 7 .
9 81
Отметим, что корни связаны с золотым сечением соотношениями
z1 1 1 2 3 4 2,618, z2 1 2 5 3 4 0,382, z1 2 z2 1 1, z1 5 z2 1 23 5 1.
Общее решение системы имеет вид:
1un 2
132 n
112 n
6v 7 4 c1 691 7
z1 5 c2 69 837
z2 .
9 n
Коэффициенты с1 и с2 находим из начальных условий и0 = 2, v0 = 1:
4c11 2 c2 3 2,
21 2 1
251
c1 3 2
, c2 3 2
.
6
1 21
1 21
7c1 5 c21 3 1,
91
Отсюда
un 4
(21 2 1)1 n 2 3 1 n
z1 2 2
z2 ,
12 2 1
1 21
vn 4
21 2 1 n (2 3 1)1 n
z1 3 2
z2 .
12 2 1
1 21
Формулу для Rn получаем как отношение
Rn 4
un
(21 2 1)1 z1n 2 (2 3 1)z2n
4
.
vn (21 2 1) z1n 3 (2 3 1)1 z2n
Упростим эту формулу, используя уравнение золотого сече!
ния 12 2 1 2 1 3 0, вытекающие из него тождества
1 1 2 3 22,
242 3
1
,
22
22 1 1 3 23,
а также приведенные выше соотношения z1 1 1 2 3, z2 1 2 4 3 :
Rn 4
1
1 2 12(n 11) 3 121 2 122(n 11)
4 2n .
2(n 11)
22(n 11)
12n 21
1
31
Последнее равенство следует из формулы Бине для чисел Фибо!
наччи, приведенной в разд. 1 (пример 10).
Случаи, когда нелинейные разностные уравнения удается решить
аналитически, сравнительно редки. На практике обычно приходит!
ся использовать численные методы либо структурные модели, по!
добные рассмотренным в разд. 4. При этом остается открытым воп!
рос о сходимости таких методов и правильности получаемых резуль!
татов. Поскольку компьютерное решение по многим причинам мо!
жет оказаться неточным или неверным, не следует жалеть времени и
усилий на его проверку (обычно эта задача не менее трудоемка, чем
само решение).
Пример 4 (Динамика эпидемии [12] ). В городе с 20 000 жителей
появляются 50 инфекционных больных, что вызывает эпидемию.
Предположим, что прирост больных за день пропорционален числу
контактов больных и здоровых, т.е. произведению числа здоровых
(еще не переболевших и не приобретших иммунитет) на число боль!
ных. Коэффициент пропорциональности (он характеризует скорость
распространения эпидемии) примем равным 10–4.
Спрашивается: как развивается эпидемия – как изо дня в день
меняется число больных? На какой день будет максимальное число
заболевших?.
Р е ш е н и е . Обозначим через х число больных, через у – число здо!
ровых, через k – коэффициент пропорциональности.
92
В соответствии с условиями задачи имеем два уравнения
x(n 1 1) 2 kx(n)y(n), k 2 0,0001,
y(n 1 1) 2 y(n) 3 kx(n)y(n).
Первое из них характеризует число заболевших в очередной (n + 1)!й
день, второе – число здоровых (еще не болевших) к (n + 1)!му дню.
Система разностных уравнений нелинейная из!за наличия произ!
ведения x(n)y(n), это не позволяет решить их аналитически.
Воспользуемся методом численного моделирования. Мы имеем
задачу с начальными условиями х(0) = 50, у(0) = 20 000, это позво!
ляет шаг за шагом вычислять изменение числа больных по дням,
полагая n = 1, 2, 3 ... .
Например, при n = 1 получаем
x1 1 0,0001 2 20 000 2 50 1 100; y1 1 20 000 3 100 1 19 900.
Программа для моделирования этой задачи в пакете MATLAB име!
ет вид:
x(1)=50; y(1)=20000;
% Количество больных и здоровых людей
k=10^(4); n = 12;
for i=1:n,
x(i+1)=k*x(i)*y(i);
y(i+1)=y(i)k*x(i)*y(i);
end;
bar(0:12,x);grid
Полученные результаты отображены столбчатой диаграммой на
рис. 5.3. Из нее видно, что критическая точка – восьмой день (3972
больных). На двенадцатый день (спад эпидемии) в городе остается
105 больных.
4000
3000
2000
1000
0
0
2
4
6
8
10
12
Рис. 5.3
93
Для оценки качества решения выполним моделирование той же
задачи в непрерывном времени. Она описывается системой обыкно!
венных дифференциальных уравнений второго порядка
x1 1 2kxy, y1 1 kxy 2 y,
где k 1 0,0001, х(0) = 50, у(0) = 20000.
Программа моделирования этой системы уравнений в пакете
MATLAB имеет вид:
T=[0 12]; X0=[20000 50]';
[t,X]=ode23('epid',T,X0);
plot(t,X(:,2)), grid
% Период наблюдения и начальные значения
% Моделирование
В команде ode23 используется вспомогательная функция xdot
function xdot=epid(t,x)
k=0.0001;
x1=k*x(1)*x(2);
y1=k*x(1)*x(2)x(2);
xdot=[x1; y1];
Результат моделирования эпидемии в непрерывном времени по!
казан плавной кривой на том же графике. Здесь критическая точка –
шестой день (3119 больных), спад эпидемии – на двенадцатый день
(209 больных). Сравнение графиков показывает, что погрешность
результатов компьютерного моделирования составляет около 20%.
Она может быть уменьшена путем корректировки дискретной мате!
матической модели.
Пример 5. Члены последовательности {xn} связаны рекуррентной
зависимостью
xn = 1 – n xn–1.
(5.8)
Она представляет собой нестационарное линейное разностное урав!
нение первого порядка.
Требуется составить программу для вычисления членов этой пос!
ледовательности при 1 1 n 1 50 , если известно, что x1 = 1/e.
Очевидный путь состоит в построении программы с простым цик!
лом, вычисляющей члены x2, x3 и т. д., непосредственно по формуле
(5.8). Приведем ее вариант, написанный на языке пакета MATLAB:
x=zeros(1,50); x(1) = exp(1);
for i=2:50,
x(i)=1i*x(i1);
end;
Вычисление первых четырех членов дает:
x1 = 1/e = 0,3678794, x2 = 1–2 x1 = 0,2642411, x3 = 1–3 x2 =
= 0,2072766, x4 = 1–4 x3 = 0,1708934.
94
Видно, что мы имеем дело с монотонно убывающей последователь!
ностью, более того, можно показать, что она стремится к нулю. Од!
нако, производя вычисления с двойной точностью в пакете MATLAB,
при n = 20 получаем отрицательное число, хотя в самом деле все чле!
ны исходной последовательности строго положительны.
Попытка проверить качество вычислений путем подстановки вы!
численных членов последовательности в формулу (5.8) только вво!
дит в заблуждение: равенство будет выполнено с высокой степенью
точности при всех n, в том числе и при n = 20. Таким образом “очевид!
ный” алгоритм оказывается практически неработоспособным. Это
лишний раз говорит о необходимости критически относиться к ре!
зультатам компьютерных вычислений. Более подробный разбор это!
го примера содержится в учебном пособии [10].
5.3. Анализ стационарных точек и точек бифуркации
В тех случаях, когда нелинейные разностные уравнения не имеют
аналитического решения, наряду с численными методами применя!
ют методы качественного анализа. Один из важных вопросов, с кото!
рого обычно начинают исследование разностных уравнений (5.2) или
(5.3) – это анализ стационарных точек. Так называются точки, удов!
летворяющие уравнению X = Ф(X) в первом случае и x 1 f (x, r ) – во
втором. Они характеризуют положения равновесия системы, когда
ее следующее состояние совпадает с предыдущим. Отметим, что ста!
ционарные точки разностного уравнения совпадают с так называе!
мыми неподвижными точками отображений Ф(X) и f (x) .
Различают два типа стационарных точек: притягивающие (устой!
чивые) и отталкивающие (неустойчивые). Поясним их смысл на при!
мере поиска корней нелинейного алгебраического уравнения мето!
дом простых итераций. Как уже отмечалось, в соответствии с этим ме!
тодом для численного нахождения корней исходное алгебраическое
уравнение F(x) = 0 переписывается в виде x = f(x), и затем ему сопос!
тавляется разностное уравнение xn+1 = f(xn), n = 0, 1, 2, ... . Тогда
корни функции F – это неподвижные точки функции f.
Пример 6. Рассмотрим кубическое уравнение
2(х – 1)2х – х2 = 0.
Оно имеет три вещественных корня 0, 1/2, 2.
x2
Для отыскания корней перепишем уравнение в виде x 1
и
2(x 2 1)2
организуем простые итерации по формуле
xn 11 1
xn2
2(xn 2 1)2
.
95
Это разностное уравнение первого порядка. Для определения его
стационарных точек (неподвижных точек функции f) построим гра!
фик в плоскости (y, х), где y 1
x2
(рис. 5.4).
2(x 2 1)2
y
2,5
2
1,5
1
x
0,5
–1
0
1
2
3
4
5
Рис.5.4
Он пересекается с биссектрисой y = х в трех точках (0,0), (1/2, 1/2),
(2,2). Анализ итерационного процесса показывает, что первая из них –
притягивающая, две другие – отталкивающие.
Если, например, взять начальное значение х0 = –1, то вычисляя сле!
дующие значения хi по формуле (5.3), получим х1 = 1/8, х2 = 1/98, ...,
т. е. процесс быстро сходится к стационарной точке (0, 0).
Если же взять начальное значение вблизи точки 2, например, х0 =
=1,9, то последующие значения будут удаляться от нее: х1 = 1,82,
х2 = 2,46, ..., т. е. эта стационарная точка – неустойчивая (отталки!
вающая). На рис. 5.4 соответствующий итерационный процесс по!
казан раскручивающейся прямоугольной спиралью.
Аналогично устанавливается неустойчивость стационарной точ!
ки (0,5; 0,5). Заметим, что область сходимости к точке 0 – полу!
ось 12 3 x 3 1/2. Следовательно, описанный вариант метода простой
итерации позволит найти только один из трех корней.
Приведем общее правило (критерий), позволяющее определять
устойчивость стационарных точек разностного уравнения (5.3). Оно
опирается на вычисление значения производной функции f(x, r) в
стационарной точке и состоит в следующем.
Критерий устойчивости стационарной точки. Пусть хi – стаци!
онарная точка уравнения xn+1 = f(xn). Если f 3 1 xi 2 4 1, то эта точка –
устойчивая (притягивающая), если же f 3 1 xi 2 4 1, то стационарная
точка xi – неустойчивая (отталкивающая). Случай, когда модуль
производной функции f(x) в стационарной точке равен нулю, требует
96
дополнительного исследования (анализа производных высших по!
рядков).
Применим этот критерий к рассмотренному примеру 6. В этом при!
мере функция f(x) и ее производная имеют вид:
f (x) 3
1x
x2
, f 2(x) 3
.
2(x 1 1)
(x 1 1)3
Вычисляя значение производной в стационарных точках х = 0;
0,5; 2, получаем:
f 1(0) 2 0, f 1(0,5) 2 34, f 1(2) 2 32.
В соответствии с критерием заключаем, что стационарная точка
х1 = 0 – устойчива, а точки х2 = 0,5; х3 = 2 – неустойчивы.
Наряду с исследованием стационарных точек нелинейных систем
важное значение имеет анализ точек бифуркации. Так называют зна!
чения параметров системы, при которых происходит качественное
изменение ее поведения, например, потеря устойчивости, возникно!
вение колебательного режима, появление предельного цикла и т. д.
Знание стационарных точек и точек бифуркации несет большую
информацию о системе и возможных типах ее поведения, что позво!
ляет в ряде случаев избежать трудоемкой процедуры нахождения
численного решения нелинейных уравнений.
Обратимся к разностному уравнению (5.3). В данном случае точ!
ками бифуркации являются критические значения параметра r, при
которых будут появляться (или исчезать) стационарные точки сис!
темы, либо меняться тип их устойчивости. Проиллюстрируем мето!
дику их отыскания на примере анализа логистического отображения
(5.4).
Пример 7 [Динамика популяции (логистическое отображение)].
Рассмотрим один вид животных или насекомых, живущий на ог!
раниченном пространстве. Обозначим через х текущую нормирован!
ную численность популяции, тогда в соответствии с классической
математической моделью ее динамика будет описываться разностным
уравнением (5.5):
xn+1 = rxn (1–xn).
Оно означает, что численность популяции в n+1!м году пропорци!
ональна численности популяции в n!м году и свободной части жиз!
ненного пространства (1–xn), значение xn изменяется от 0 до 1. По!
ложительный коэффициент r зависит от плодовитости популяции,
условий питания и реальной площади.
97
Найдем стационарные точки этой модели, рассмотрев алгебраи!
ческое уравнение x = rx(1–x). У него есть очевидный корень x1 = 0.
Если численность популяции в начальный момент равна 0, то и она
1
будет равняться 0 и далее. Второй корень равен x2 1 1 2 . Таким об!
r
разом, x1 и x2 – стационарные точки данной задачи. Проанализируем
их устойчивость и найдем точки бифуркации (критические значения
параметра r). Это позволит ответить на вопрос, как будет себя вести
популяция при разных значениях r.
Рассмотрим отдельно случаи 0 < r < 1, 1 < r < 3, 3 < r < 4, полагая
везде, что 0 1 x 1 1.
С л у ч а й 1 (0 < r < 1). В этом случае имеется одна стационарная
точка х1 = 0 (корень х2 отрицателен). Чтобы проанализировать ее
устойчивость, нарисуем кривую y 1 f (x) , где f(x) = rx(1–х) и прямую
y = x (рис. 5.5). Отложим х1 по оси абсцисс, проведем вертикаль до
пересечения с кривой y 1 f (x) (точка А), затем из нее горизонталь до
пересечения с осью х. Полученную точку пересечения обозначим че!
рез х2. Легко проверить, что x2 1 f (x1 ). Взяв точку х2 за начальную и
повторив все те же операции, получим х3, затем х4 и т. д. Эта проце!
дура называется построением лестницы Ламерея. Она позволяет гра!
фически строить члены последовательности {xn}.
y
0,4
0,3
r=1
B
A
0,2
C
0,1
0
x
0
0,2
0,6
1
Рис. 5.5
В данном случае эта последовательность стремится к нулю (спуск
по лестнице Ламерея приводит в начало координат), таким образом,
стационарная точка x1 = 0 – устойчивая. К этому же выводу прихо!
дим, вычисляя производную f 1(x) 2 r 3 2rx. Поскольку f 1(x1) 2 r 3 1, то
условие устойчивости выполняется.
С л у ч а й 2 (1 < r < 3). Как только коэффициент r становится
больше единицы, стационарная точка x1 = 0 теряет устойчивость,
так как f 1(x1) 2 1. Одновременно появляется вторая стационарная
1
точка x2 1 1 2 . Вычисляя для нее f 1(x), получаем f 1(x2 ) 2 2 3 r. Гра!
r
98
фик зависимости f 1(x) от r приведен на рис. 5.6. Область устойчиво!
сти f 1(x2 ) 2 1 на нем выделена штриховкой. Поскольку при 1 < r < 3
модуль производной меньше единицы, то на указанном интервале
стационарная точка х2 – устойчивая. Критическое значение r1 = 1,
при котором появляется стационарная точка x2 и меняется тип ус!
тойчивости стационарной точки x1, представляет собой первую точ!
ку бифуркации.
Построив лестницу Ламерея в этом случае (рис. 5.7), мы обнару!
живаем, что движение по ней идет теперь в направлении точки x2 . В
применении к исходной биологической задаче это означает, что чис!
ленность популяции по прошествии нескольких лет стабилизирует!
ся и перестанет меняться со временем.
2
f1 ( x2 )
1
r
1
x2
f
0,5
0,4
2
3
0,3
0,2
0,1
–1
x
0
0
0,2
Рис. 5.6
0,6
1
Рис. 5.7
По мере приближения ступенек лестницы Ламерея к точке x2 гори!
зонтальные размеры ступенек стремятся к геометрической прогрессии
lim
xn 11 1 xn
3 c 3 f 2(x) x 2x .
2
1 xn 51
n 34 x
n
Отсюда, кстати, вытекает обоснование приведенного выше крите!
рия устойчивости стационарных точек. Отметим, что все стационар!
ные точки лежат на пересечении графика функции f (x) с биссектри!
сой первого квадранта (рис. 5.7).
С л у ч а й 3 (3 < r < 4). Как только r превышает значение 3, ста!
ционарная точка x2 теряет устойчивость, т. е. значения r2 = 3 – это
вторая точка бифуркации. В поведении системы происходят каче!
ственные изменения – в ней возникают устойчивые колебания. Их
период сначала равен двум, затем, по мере увеличения параметра r,
он удваивается и т. д. Значения ri, при которых происходят эти удво!
ения периода, являются точками бифуркации. По мере приближе!
ния r к критическому значению (оно примерно равно 3,5699) число
бифуркаций удвоения периода стремится к бесконечности, а при дос!
99
тижении этого значения всякая регулярность в поведении системы
пропадает и возникает так называемая хаотическая последователь!
ность {xn}. Поведение системы в интервале 3,57 < r < 4 известно, как
перемежающийся хаос, а при r > 4 – как сплошной хаос.
Описанная картина иллюстрируется графиками, приведенными
на рис. 5.8–5.10. На них показано изменение численности популя!
ции во времени при различных значениях коэффициента r и разных
начальных условиях. Для удобства дискретные точки на графиках
соединены прямыми линиями.
С л у ч а й 1 (рис. 5.8): 1 < r < 3. График на рис. 5.8, а построен
для r = 1,2; х1 = 1/7, 1 1 n 1 25. Он получен в пакете MATLAB с помо!
щью команд
x(1) = 1/7; for i = 1:25; x(i+1) = 1.2*x(i)*(1x(i)); end, plot(x)
Видно, что численность популяции монотонно растет, стремясь к
1 1
1 . На рис. 5.8, б величи!
r 6
на r = 2,9 немного меньше бифуркационного значения r2 1 3. После!
довательность {xn} по!прежнему стремится к равновесной точке
1
x2 1 1 2 , однако характер сходимости заметно изменился.
r
а)
б)
устойчивой стационарной точке x2 1 1 2
0,17
0,7
0,16
0,6
0,15
0,14
0,5
0
10
20
30
0,4
0
5
10
15
Рис. 5.8
С л у ч а й 2 (рис. 5.9): 3< r < 4. При r = 3,1 наблюдается периоди!
ческое колебание численности между двумя постоянными уровнями
(рис. 5.9, а). При r = 3,5 число уровней становится равным 4, проис!
ходит бифуркация удвоения периода (рис. 5.9, б).
100
С л у ч а й 3 (рис. 5.10): r = 4. График на рис 5.10 демонстрирует
появление хаотических колебаний.
а)
б)
0,8
0,9
0,7
0,7
0,6
0,5
0,5
0,4
0
5
10
15
20
0,3
0
10
20
30
Рис. 5.9
0,9
0,7
0,5
0,3
0,1
0
20
40
60
80
Рис. 5.10
В рассмотренном примере мы столкнулись с простейшим нелиней!
ным разностным уравнением вида xn+1 = f(xn) и убедились, что даже
при простых функциях f(x) оно демонстрирует сложное поведение.
Заметим, что хаотическая динамика может возникнуть и в непре!
рывных системах. Наиболее известный пример из теории дифферен!
циальных уравнений – так называемая система Лоренца. Она пред!
ставляет собой систему трех нелинейных дифференциальных урав!
101
нений, в которой при определенном сочетании коэффициентов воз!
никает удивительный режим колебаний, получивший название
«странного аттрактора». Траектории системы в трехмерном про!
странстве непредсказуемым образом блуждают вокруг двух особых
точек.
Любопытно, что такое поведение непрерывных систем возможно
только для систем третьего и более высоких порядков, в то время как
у дискретных систем явление детерминированного хаоса возникает
уже в системах первого порядка. Это связано с тем, что траектории
дифференциальных уравнений обязаны удовлетворять некоторым
требованиям гладкости и непрерывности, иначе производные, вхо!
дящие в эти уравнения, могут потерять смысл. В дискретных систе!
мах ограничения такого рода отсутствуют, поэтому они могут демон!
стрировать более «экстравагантное» поведение.
5.4. Задачи и упражнения
1. Найти решение нелинейного уравнения yk11 1
1
, если изве!
2 2 yk
стно значение y0.
z 11
У к а з а н и е . Замена yk 2 k
приводит к линейному уравнению
zk
zk : zk 11 1 zk 2 1.
О т в е т . yk 3
y0 1 k(1 2 y0 )
.
1 1 k(1 2 y0 )
2. Решить нелинейное разностное уравнение
yk 11 1
yk
, y0 1 1.
1 2 yk
1
.
k 21
3. Решить нелинейное разностное уравнение
О т в е т . yk 1
ykyk11 1 yk 2 yk11, y1 1 1.
1
(гармоническая последовательность).
k
4. Решить нелинейное уравнение
О т в е т . yk 1
yk 11 1 2 2
102
1
, y0 1 2.
yk
k12
.
k 11
5. Решить линейное нестационарное разностное уравнение
О т в е т . yk 2
1
yk11 3 ayk 4 0, a 4 16 1 5 27 eh , y0 4 0,
k9
8
в котором коэффициент а зависит от номера такта.
У к а з а н и е . Перейти к линейному стационарному уравнению вто!
рого порядка.
6. Решить нелинейное разностное уравнение
yk11 1 2
ykyk 12
, y1 1 1, y2 1 0,5.
yk 2 yk 12
У к а з а н и е . Замена yk 1
1
приводит к линейному уравнению вто!
zk
рого порядка.
7. Существуют три классических определения среднего значения
двух чисел a и b: среднее арифметическое, среднее геометрическое и
среднее гармоническое
a1b
,
2
a b,
2
.
1 1
1
a b
В любой арифметической прогрессии каждый член представляет
собой среднее арифметическое соседних членов, в геометрической
прогрессии – среднее геометрическое соседних членов. Обе они опи!
сываются линейными разностными уравнениями. Последователь!
ность, начинающаяся с чисел 1 и 1/2, называется гармонической,
если каждый ее член представляет собой среднее гармоническое со!
седних членов. Постройте эту последовательность, получите для нее
разностное уравнение и составьте разностное уравнение для ее част!
ных сумм sk.
О т в е т . Возможны несколько вариантов уравнений, приведем два
из них:
1) y
2
k 11
1
1
1
1
1
2 ,
1
2 1;
yk12 yk sk11 3 sk sk 3 sk21
1
1
1
2) y 1 y 2 1, sk11 3 sk 1 k 2 1.
k 11
k
103
8. Человек держит за конец резиновый жгут длиной 1 м, привя!
занный к дереву. У другого конца жгута сидит жук. Каждую секунду!
жук проползает 1 см по жгуту. Каждую секунду человек, держа ко!
нец жгута, удаляется от дерева на 1 м. Доползет ли жук до человека
и если да, то за какое время?
Р е ш е н и е . Доля жгута в процентах, пройденная жуком за n се!
кунд, описывается разностным уравнением S(n) = S(n–1)+1/n, откуда
S(n) = 1+1/2+1/3+1/4+...+1/n.
Известно, что с ростом n эта сумма
3
неограниченно возрастает. Для ре!
2
шения надо найти такое n, чтобы
величина S(n) достигла (или превы!
1
сила) 100%. Заметим, что попыт!
0
ка прямого компьютерного опреде!
0
2
4
6
8 10 12 14
ления искомого n с помощью
MATLAB!программы типа n = 1;
S=1;while S<100,n=n+1;S=S+1/n;end оказывается безуспешной, по!
скольку ряд возрастает очень медленно и компьютер «зависает».
Воспользуемся формулой Эйлера для оценки частных сумм гармо!
нического ряда S(n) » ln n + 0,577215, которая при n > 20 обеспечи!
вает очень высокую точность. Подставляя в нее S(n) = 100, получа!
ем: ln n = 100 – 0,577215 = 99,42, откуда n » e99,42 с.
О т в е т . Теоретически жук доползет до человека, однако на это
потребуется около 4,8·1032 тысячелетий. И место встречи будет на!
ходиться слишком далеко.
9. Исследовать нелинейное разностное уравнение второго поряд!
x 11
. Построить график решения для 1 < n < 15 при х1 = 1,
ка xn 11 2 n
xn 21
х2 = 2.
О т в е т . График решения при х0 = 1, х1 = 2 приведен на рисунке.
MATLAB – программа для его построения:
clear x;x(1) = 1;x(2) = 2; for i = 2:14; x(i+1) = (x(i)+1)/x(i1);
end, plot(x)
10. Дано нелинейное разностное уравнение второго поряд!
ка xn 11xn 21 1 xn 2 1. Найдите х2007, если х1 = 1799, х2 = 1828 (это
годы рождения наших великих писателей).
У к а з а н и е . При любых положительных начальных условиях
х0 = а, х1 = b, решение – периодическая функция с периодом 5, т. е.
начиная с шестого номера, значение членов повторяются (докажите это).
О т в е т . Поскольку остаток от деления 2007 на 5 равен 2, получа!
ем х2002 = х2 = 1828.
104
11. Рассчитайте первые 50 точек решения линейного неоднород!
ного нестационарного разностного уравнения первого порядка
xn 11 1 xn 1 nxn 2 e, x1 2 1. Постройте график решения. Как оно ведет
себя при больших значениях n?
П р е д у п р е ж д е н и е . Будьте осторожны – попытка рекуррент!
ного расчета точек х2, х3, ... на компьютере может ввести в заблуж!
дение.
О т в е т . Последовательность {xn} положительна и монотонно убы!
вает, стремясь при больших значениях n к нулю по гиперболическо!
e
.
n12
12. Найти стационарные точки нелинейного разностного уравнения
му закону вида
xn 11 1 cos xn
и выяснить их тип. Построить график решения для 1 1 n 1 40 при
х1 = 1.
У к а з а н и е . Постройте графики «итерированного» косинуса
coscos...cos x и посмотрите, к чему они стремятся.
13. Исследовать нелинейное разностное уравнение с «треуголь!
ным» отображением f(x):
xn+1 = r(1 – |1 – 2 xn|) ; x0 = 0,1; 0 < x < 1.
Найти его стационарные точки и точки бифуркации. Построить
графики решения при r = 1 для начальных условий x0 = 1/7, 1/10 и
1/11. Убедитесь, что при этом имеют место периодические колеба!
ния между двумя, четырьмя или восемью постоянными уровнями
соответственно. Проверьте, что при r = 0,9 имеют место хаотические
колебания.
105
Библиографический список
1. Гантмахер Ф. Р. Теория матриц. М.: Наука, 1988.
2. Гноенский Л. С., Каменский Г. А., Эльсгольц Л. Э. Математичес!
кие основы теории управляемых систем. М.: Наука, 1969.
3. Деч Г. Руководство к практическому применению преобразова!
ния Лапласа и z!преобразования. М.: Наука, 1971.
4. Директор С., Рорер Р. Введение в теорию систем. М.: Мир, 1974.
5. Есипов А. А., Сазонов Л. И., Юдович В. И. Практикум по обык!
новенным дифференциальным уравнениям. М.: Вузовская книга,
2001.
6. Малинецкий Г. Г. Хаос: структуры, вычислительный экспери!
мент. Введение в нелинейную динамику. М.: Эдиториал УРСС, 2000.
7. Мальцев А. И. Основы линейной алгебры. М.: Наука, 1970.
8. Маркушевич А. И. Возвратные последовательности. М.: Наука,
1975.
9. Мироновский Л. А., Юдович В. С. Аналоговые и гибридные вы!
числительные машины: Учеб. пособие/ ЛИАП. Л., 1991.
10. Мироновский Л. А. Моделирование динамических систем:
Учеб. пособие/ ГААП. СПб., 1992.
11. Михайлов В. В., Бритов Г. С., Мироновский Л. А. Аналоговые
вычислительные машины: Учеб. пособие/ ЛИАП. Л., 1973.
12. Очков В. Ф. MATHCAD!8 Pro для студентов и инженеров. М.:
Компьютер Пресс, 1999.
13. Романко В. К. Курс дифференциальных уравнений и вариаци!
онного исчисления. М.: Лаб. базовых знаний, 2000.
14. Чен К., Джибмю П., Ирвинг А. MATLAB в математических
исследованиях. М.: Мир, 2001.
106
Оглавление
Предисловие ....................................................................
3
1. Основные свойства разностных уравнений ........................
5
1.1. Примеры разностных уравнений ............................
5
1.2. Решение однородных разностных уравнений ...........
8
1.3. Системы линейных разностных уравнений .............. 17
1.4. Синтез разностных уравнений для получения задан!
ных функций ...................................................... 24
1.5. Задачи и упражнения .......................................... 31
2. Применение жордановых цепочек векторов при решении
дифференциальных и разностных уравнений ..................... 33
2.1. Корневые векторы и циклические пространства ....... 33
2.2. Жордановы цепочки векторов ............................... 35
2.3. Дифференциальные уравнения и жордановы цепочки
векторов ............................................................ 39
2.4. Разностные уравнения и цепочки корневых векторов
47
2.5. Задачи и упражнения .......................................... 55
3. Разностные уравнения и z!преобразование ........................ 58
3.1. Определение z!преобразования .............................. 58
3.2. Решение разностных уравнений с помощью z!преоб!
разования .......................................................... 60
3.3. Цифровые интеграторы и дискретные передаточные
функции ............................................................ 65
3.4. Задачи и упражнения .......................................... 71
4. Структурные модели разностных уравнений ...................... 73
4.1. Модели уравнений с начальными условиями ........... 73
4.2. Пример моделирования в пакете SIMULINK ............ 75
4.3. Модели уравнений с краевыми условиями ............... 77
4.4. Матричное решение разностных уравнений ............. 79
4.5. Задачи и упражнения .......................................... 82
5. Нелинейные разностные уравнения ................................. 83
5.1. Источники нелинейных разностных уравнений ....... 83
5.2. Решение нелинейных разностных уравнений ........... 87
5.3. Анализ стационарных точек и точек бифуркации ..... 95
5.4. Задачи и упражнения .......................................... 102
Библиографический список ................................................ 106
107
Учебное издание
Мироновский Леонид Алексеевич
МОДЕЛИРОВАНИЕ
РАЗНОСТНЫХ
УРАВНЕНИЙ
Учебное пособие
Редактор Г. Д. Бакастова
Компьютерная верстка А. Н. Колешко
Сдано в набор 17.05.04. Подписано к печати 16.11.04. Формат 60´84 1/16. Бумага офсетная.
Печать офсетная. Усл. печ. л. 6,27. Усл. кр.!отт. 6,4. Уч. !изд. л. 6,15. Тираж 100 экз. Заказ №
Редакционно!издательский отдел
Отдел электронных публикаций и библиографии библиотеки
Отдел оперативной полиграфии
СПбГУАП
190000, Санкт!Петербург, ул. Б. Морская, 67
Документ
Категория
Без категории
Просмотров
4
Размер файла
749 Кб
Теги
mironovsky
1/--страниц
Пожаловаться на содержимое документа