35559

МОДЕЛИРОВАНИЕ РАЗНОСТНЫХ УРАВНЕНИЙ

Книга

Математика и математический анализ

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

Русский

2013-09-17

2.4 MB

246 чел.

PAGE  3

Мироновский Л.А.

МОДЕЛИРОВАНИЕ  РАЗНОСТНЫХ

УРАВНЕНИЙ

Санкт-Петербург

2003 г.


МИНИСТЕРСТВО ОБРАЗОВАНИЯ РОССИЙСКОЙ ФЕДЕРАЦИИ

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

профессионального образования

«Санкт-Петербургский государственный университет

аэрокосмического приборостроения»

Л.А. Мироновский

МОДЕЛИРОВАНИЕ

РАЗНОСТНЫХ

УРАВНЕНИЙ

Учебное пособие

Санкт-Петербург

2003 г.


ОГЛАВЛЕНИЕ

  1.  ОСНОВНЫЕ СВОЙСТВА РАЗНОСТНЫХ УРАВНЕНИЙ
    1.  Примеры  разностных уравнений……………………………………………….4
    2.  Решение однородных разностных уравнений …………………………………6
    3.  Системы линейных разностных уравнений…………………………………..13
    4.  Синтез разностных уравнений для получения заданных функций…………18
    5.  Задачи и упражнения…………………………………………………..……….24

2. ПРИМЕНЕНИЕ КОРНЕВЫХ ВЕКТОРОВ ПРИ РЕШЕНИИ

ДИФФЕРЕНЦИАЛЬНЫХ И РАЗНОСТНЫХ УРАВНЕНИЙ

2.1.  Корневые векторы и циклические пространства …………………………….25

2.2.  Жордановы цепочки векторов…………………………………………………26

2.3.  Дифференциальные уравнения и жордановы цепочки векторов……………30

2.4.  Разностные уравнения и цепочки корневых векторов………………………36

2.5.  Задачи и упражнения…………………………………………………………...42

3. РАЗНОСТНЫЕ УРАВНЕНИЯ И Z-ПРЕОБРАЗОВАНИЕ

3.1. Определение  z-преобразования……………………………………………….46

3.2. Решение разностных уравнений с помощью z-преобразования…………… 47

3.3. Цифровые интеграторы и дискретные передаточные функции ……….……52

3.4. Задачи и упражнения………………………………………………………...…56

4. СТРУКТУРНЫЕ МОДЕЛИ РАЗНОСТНЫХ УРАВНЕНИЙ

4.1. Модели уравнений с начальными условиями ………………….……………58

4.2. Пример моделирования в пакете SIMULINK ………………………………..58

4.3. Модели уравнений с краевыми условиями……………………………………61

4.4. Матричное решение разностных уравнений…………………………………..63

4.5. Задачи и упражнения…………………………..………………………………..65

5.  НЕЛИНЕЙНЫЕ РАЗНОСТНЫЕ УРАВНЕНИЯ

5.1.  Источники нелинейных разностных уравнений………………………………66

5.2.  Решение нелинейных разностных уравнений…………………………………68

5.3.  Анализ стационарных точек и точек бифуркации…………………………….74

5.4.  Задачи и упражнения……………………………………………………………79

Библиографический список…………………………………………………………………….81

  1.  
    ОСНОВНЫЕ СВОЙСТВА РАЗНОСТНЫХ УРАВНЕНИЙ

1.1. Примеры разностных уравнений

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

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

,

а разностные уравнения – значения функции в различные моменты времени

.    (1.1)

Перечислим некоторые источники  разностных уравнений: дискретизация обыкновенных дифференциальных уравнений и уравнений в частных производных;  модели объектов с дискретным временем (задача Фибоначчи о размножении кроликов); объекты с дискретным пространством (R-2R-цепь); анализ математических рядов и рекуррентных соотношений.

Пример 1. Непрерывную функцию    можно описать с помощью дифференциального уравнения

.

Если рассматривать дискретные моменты времени  , то эту же функцию можно описать с помощью разностного уравнения

.   (1.2)

Непосредственной подстановкой легко убедиться, что    удовлетворяет обоим уравнениям.

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

Линейное однородное разностное уравнение n-го порядка с постоянными коэффициентами имеет вид:

.   (1.3)

Решением этого уравнения называется решетчатая функция x(t), t = 0, 1, 2, ...,  обращающая уравнение в тождество. Можно говорить и о непрерывном решении – функции x(t), удовлетворяющей уравнению (1.3)  при любом  .

Иногда, чтобы подчеркнуть дискретный характер изменения времени, уравнение (1.3) записывают в форме

.

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

Пример 2. Пусть имеется геометрическая прогрессия   Ее можно описать однородным разностным уравнением первого порядка

.    (1.4)

Действительно, подставляя последовательно  , получаем

.

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

Пример 3. Пусть имеется арифметическая прогрессия . Ее можно описать неоднородным разностным уравнением первого порядка

    (1.5)

или однородным уравнением второго порядка

.     (1.6)

(Убедитесь, что это действительно так).

Пример 4. Предположим, что необходимо вычислить напряжение un  на выходе электрической цепи  (рис. 1.1), зная входное напряжение u0  и полагая R=r . Такие цепи встречаются, например, в цифро-аналоговых преобразователях.

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

Рис. 1.1

Запишем уравнение токов для первого узла цепи:

.

Аналогично, для  k+1-го узла получаем

   (1.7)

или учитывая, что R = r

.

Уравнение (1.7) представляет собой разностное уравнение второго порядка, поскольку связывает значение функции u  в точке  k+2  с ее значениями в двух предыдущих точках  k+1 и k.

Пример 5. Рассмотрим старинную задачу о размножении кроликов (задача Фибоначчи или Леонардо Пизанского; итальянского математика, жившего около 1200г.). В этой задаче требуется определить число пар зрелых кроликов, образовавшихся от одной пары в течение года, если каждая пара кроликов рождает ежемесячно новую пару, и новорожденные достигают зрелости через месяц. Таким образом, получается некоторая последовательность; нас интересует, каким уравнением она описывается.

В первый момент времени число пар равно  x0=1.  Через месяц прибавится пара новорожденных, но число зрелых пар не изменится  x1=1.

Через два месяца крольчата достигнут зрелости, и общее число пар будет  x2=2.

Через k месяцев число зрелых пар будет xk,  а через k+1 месяцев  xk+1 ,  но так как к этому времени появится  xk  пар приплода, то через  k+2  месяцев общее число зрелых пар будет

.       (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)  также является решением.

Общее решение имеет вид

,     (1.10)

где  c1, ..., cn – произвольные коэффициенты.

Известны два основных метода решения линейных разностных уравнений – с помощью характеристического полинома и с использованием  z-преобразования, аналогичного преобразованию Лапласа.

Рассмотрим первый из них. Будем искать решение в виде  ,  где z – некоторое число. Если  t – дискретно, то x(t) – геометрическая прогрессия, если же t – непрерывно, то   x(t) – показательная функция. В обоих случаях ее можно записать в виде  et,  где   = ln z.

Подставляя x(t) в (1.3) и сокращая на  zt,  получаем алгебраическое уравнение

,      (1.11)

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

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

а). Вещественные корни

Если уравнение (1.11) имеет вещественные и различные корни z1, ..., zn, то получаем n линейно независимых частных решений уравнения (1.3):

.

Общее решение в соответствии с (1.10) имеет вид:

.     (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.

Оно имеет вещественные и различные корни   .

Следовательно, общее решение имеет вид

uk =c1 2,62k + c2 0,38k.    (1.13)

Пример 3. Найдем общее решение уравнения Фибоначчи (1.8).

Составляем характеристическое уравнение

z2z – 1 = 0.

Корни этого уравнения также вещественны и различны

.

Общее решение уравнения Фибоначчи имеет вид

.     (1.14)

Пример 4. Решим разностное  уравнение

x(t+2) – 3x (t+1) + 2x(t) = 0.

Составляем характеристическое уравнение:  z2 – 3z + 2 = 0.

Находим его корни:  z1 =1,      z2 = 2.   Общее решение:  x = c1 + c2 2t.

б) Комплексные корни

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

,

где   (сравните с аналогичной формулой для дифференциальных уравнений).

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

Пример 5. Решим разностное уравнение

x(t+2) + 2x(t+1) + 4x(t) = 0.

Характеристическое уравнение:  z2 + 2z + 4 = 0.

Его корни:  .  Модуль и аргумент корней:  .

Общее решение:  

в)  Кратные корни

Если один из корней характеристического уравнения  zi  имеет кратность  k,  то соответствующее ему слагаемое в общем решении умножается на полином степени  k–1

(это полностью с аналогично случаю дифференциальных уравнений).

Пример 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.

Общее решение:

г) Учет начальных и краевых условий

Общее решение, в которое входят произвольные постоянные  ci, задает семейство решений разностного уравнения. Для выделения из этого семейства нужного нам частного решения, необходимо найти эти постоянные. В случае дифференциальных уравнений их находят, задавая начальные или краевые условия, общее число которых равно порядку уравнения. Аналогично для разностных уравнений для нахождения постоянных  c1, ..., cn  надо задать n условий, связывающих значения  x(t)  в различных точках.

Чаще всего встречаются два типа условий:

  •  когда задается начальная последовательность  x(1), x(2), ... , x(n);  
  •  когда часть значений x(t) задается в начале, а часть – в конце интервала решения.

Условия второго типа называются краевыми.

Пример 8. Требуется решить разностное уравнение    для трех вариантов начальной последовательности

Решение. Составляем характеристическое уравнение и находим его корни:

Общее решение имеет вид    

Постоянные  с1  и с2  находим из начальных условий.

В случае а) имеем    с1 + с2 = 1,   с1 + 2с2 = 3,  откуда с1 = – 1,   с2 = 2, т.е.   

В случае б) имеем   с1 + 2с2 = 1,  с1 + 4с2 = 3,  откуда с1 = – 1,   с2 = 1, т.е.  

В случае в) имеем   с1 + 2с2 = 2,    с1 + 4с2 = 0,  откуда с1 = 4,  с2 = –1, т.е.  

Пример 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  + c2 .  Начальные условия:   x0 = 1,  x1 = 1.

Cистема уравнений для определения  с1  и  с2:

              откуда               .

Решение задачи Фибоначчи (общий член последовательности Фибоначчи):

.

Не выполняя округлений, то же решение можно записать в виде

.

Это известная формула Бине. Хотя в нее входят радикалы, она дает при любом  k  целые числа. В частности при k = 0  и  k = 1  получаем  x0 = x1 = 1,  а при  k = 12  (решение задачи о кроликах) имеем  x12 = 233.

Пример 11. Требуется вычислить определитель трехдиагональной  nxn  матрицы  A  вида

Для  n = 1, 2, 3  находим

Чтобы получить общую формулу, раскроем определитель матрицы  An+2  по элементам первой строки:

.

Обозначая  , приходим к разностному уравнению  .

Выписываем характеристическое уравнение  z2z +1 = 0  и находим его корни

,      и   

Общее решение имеет вид:   

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

Ответ.    .

Результат выглядит неожиданно и вряд ли может быть "угадан" путем рассмотрения значений определителей для  n = 1, 2, 3, … .

Пример 12. Вычислить определитель трехдиагональной  nxn  матрицы вида

Для  n = 1, 2  находим

Записываем рекуррентное соотношение, связывающее три последовательных определителя

.

Ему соответствует  разностное уравнение

Корни характеристического уравнения   z1 = z2 =2.  

Общее решение имеет вид:   .

Находим коэффициенты  с1  и  с2:    

Ответ.    .

Пример 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 = (za)(z – 2)2 = 0

и находим его корни z1 = z2 =2,  z3 = a,  

Общее решение имеет вид:   .

Находим коэффициенты  с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.

Ответ.    xn = (b+cn)2n . Удивительно, но решение не зависит от параметра  a.

Пример 14.   Разностные уравнения и многочлены Чебышева. Классические многочлены Чебышева  Tn(x)  определены на интервале  [–1,  1].  Они обладают рядом замечательных свойств. В частности, среди всех полиномов с фиксированным коэффициентом при старшей степени  х  они имеют минимальную амплитуду на интервале  [–1,  1] (полиномы, наименее уклоняющиеся от нуля). Первые четыре многочлена имеют вид

Их графики  показаны на рис. 1.2.

Многочлены могут быть заданы с помощью рекуррентного соотношения

Рис. 1.2

Оно представляет собой линейное однородное разностное уравнение второго порядка. Ему соответствует характеристическое уравнение

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

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

Например, при  n = 2  получаем

Этот многочлен имеет два нуля, расположенных в точках    и минимум в точке 0 (рис. 1.2). Значения многочлена в точке минимума и на краях интервала по абсолютной величине одинаковы и равны единице.

д) Решение неоднородных разностных уравнений

Рассмотрим линейное неоднородное разностное уравнение

,

где  fk – известная функция дискретного аргумента.

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

Пример 15. Неоднородное разностное уравнение второго порядка с комплексными корнями

Находим частное решение и корни характеристического полинома

.

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

Коэффициенты  с1,  с2  находим из начальных условий: с1 = , с2 = 2/3.

1.3. Системы линейных разностных уравнений

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

X(t+1)=AX(t), ,     (1.15)

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

Решение этой системы может быть записано в компактной степенной форме:

X(t)=At X0,      (1.16)

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

Для матриц простой структуры существует более удобный способ записи решения. В соответствии с ним сначала надо найти собственные числа и  собственные векторы матрицы  A.  Собственные числа zi получаем как корни характеристического уравнения
det(A–
zE)=0, собственные векторы Hi находим из соотношений АHi  = ziHi . Формула для общего решения имеет вид

X(t)=с1H1  + ... + сnHn.     (1.17)

Постоянные сi  определяются из начальных условий путем решения системы линейных алгебраических уравнений  X(0) = с1H1  + ... + сnHn..

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

              (1.18)

где  Н1, …, Нк  корневые векторы матрицы  А,  отвечающее собственному числу  zi,  причем  Н1 – ее собственный вектор, а  Нк  – циклический вектор, порождающий инвариантные подпространства матрицы для указанного собственного числа.

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

Пример 1.  Найдем решение системы разностных уравнений второго порядка

x(t+1) = x(t) – 3y(t)

y(t+1) = x(t) + y(t).

Выписываем характеристическое уравнение  и находим его корни  .  Собственные векторы матрицы  A:  H1=,   H2 =.

Запишем первую  компоненту решения:

H1 =      или

Выделяем вещественную и мнимую части

и строим из них общее решение

X(t)=.

Пример 2. Задача о размножении бактерий. В колонию бактерий, численность которых удваивается каждую  секунду, попадает вирус. Каждую  секунду он съедает одну бактерию, после чего размножается делением. Какова судьба колонии бактерий, если их начальная численность равна N?

Таким образом, сначала было N бактерий и один вирус. На первой секунде вирус съел одну бактерию и размножился, число бактерий также удвоилось, в результате стало 2(N – 1) бактерий и 2 вируса  и т.д.

Решение. Обозначим через  x и y  текущее число вирусов и бактерий и рассмотрим вектор .  Матричное описание этой задачи имеет вид

X(t+1)=AX(t),  ,  .

Для получения решения воспользуемся сначала степенной формулой (1.16).  Последовательно возводя матрицу  А  в степени 2, 3, ..., n, находим

 An = 2n , X(t)=At X0= 2t 

Переходя к скалярной форме записи, получаем

 x = 2 t,  y = (N–t) 2 t.

Переменная  y  обращается в нуль при  t = N,  следовательно колония вымрет через N секунд.

Второй способ решения опирается на вычисление собственных векторов матрицы  А.  В данном случае мы имеем дело со случаем кратных собственных чисел  z1 = z2 =2,  поэтому решение ищем в форме (1.18):

   (1.19)

где    – единственный собственный вектор матрицы  А.

Для определения корневого вектора    подставляем  t = 0  и используем начальные условия  

Коэффициент  с1 = –1  находим из равенства    которое получается после подстановки формулы (1.19) в исходную систему уравнений.

Следовательно, решение имеет вид

что совпадает с полученным первым способом.

Для численного решения этой задачи можно воспользоваться одним из математических пакетов. График решения, полученный с помощью пакета MATLAB для N = 10, показан на рис. 1.3. Видно, что на девятой секунде численность вирусов и бактерий совпала (точка пересечения кривых), в результате чего на десятой секунде количество бактерий стало равным нулю. Дальнейшая часть кривой y(t)  физического смысла не имеет.

Если начальная численность бактерий равна миллиону, то колония погибнет через миллион секунд (это примерно 11 суток 13 часов и 47 минут).

Рис. 1.3

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

y(t+2) – 4y(t+1) + 4y(t) = 0.

Второй путь – составить  разностное уравнение для переменной  z=y/x. Оно получается очень простым z(t+1) – z(t) = –1 и имеет линейное решение  z(t) = Nt, которое обращается в нуль при  t = N.

Пример 3. Задача Нарайаны (Индия, 14 век). Сколько коров и телок произойдет от одной коровы в течение 20 лет, если корова в начале каждого года приносит телку, а телка дает такое же потомство в начале года, достигнув трех лет.

Решение. Пусть  xn  – число коров, которые принесли потомство в начале n-го года (число новорожденных телок тоже равно  xn);  yn – число телок, которые принесут потомство через год;  zn – число телок, которые принесут потомство через два года.   

Тогда общее число голов в стаде равно  n = 2 xn + yn + zn.  

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

    (1.20)

Перепишем ее в матричной форме:

  (1.21)

Рассмотрим три подхода к решению этой системы разностных уравнений – скалярный, векторный и матричный.

Скалярный способ. Исключая переменные  yn, zn из системы (1.20), перейдем от нее к одному скалярному разностному уравнению третьего порядка

, с начальными условиями  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.

Заметим, что численность стада определяется формулой  

Значит, через 20 лет стадо будет состоять из   =2 872 + 595 + 406 = 2745  голов.

При ручном счете этот способ быстрее всего приводит к цели.

Векторный способ. Используем представление решения через собственные числа  li  и собственные векторы  Нi  матрицы  А  (1.21):

   (1.22)

Характеристическое уравнение  l3l2 –1 = 0  имеет один вещественный корень и пару комплексных корней  

Собственные числа и собственные векторы матрицы А могут быть вычислены в пакете MATLAB  с помощью команды [H, L] = eig(A):

После подстановки этих значений в формулу (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))

Запись этих корней в обычной нотации имеет вид:

  где   

.

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

Кроме того, он позволяет дать ответ на вопрос о структуре стада – процентном соотношении численности коров, телок-двухлеток и телок-однолеток при больших  n. А именно,  их численность будет соотноситься, как компоненты первого собственного вектора матрицы А, т.е. x : y : z = 0,7710 : 0,3589 : 0,5261. Другими словами, коров будет около  46,6%, телок-двухлеток – около 21,6%, телок-однолеток – около 31,8%.  Полученное выше решение
[
x, y, z] = [ 872  406  595]  для  n=19  характеризуется как раз таким соотношением.

Матричный способ. Он опирается на использование степенных формул

При работе в пакете 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:

    

Любопытно, что все они обладают скрытой симметрией, которая становится явной после перестановки второй и третьей строк.

Динамика изменения поголовья стада по годам в обычном и логарифмическом масштабах иллюстрируется рис. 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=x+y+z; t=1:20;

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.4. Синтез разностных уравнений для получения заданных функций

 

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

  1.  получение возмущающих функций при решении неоднородных разностных уравнений;
  2.  воспроизведение переменных коэффициентов при моделировании линейных разностных уравнений;
  3.  получение тестовых воздействий для исследования систем автоматического регулирования.

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

Общий прием генерирования заданной функции  y(t)  предполагает отыскание определяющего линейного разностного уравнения

,    (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) = eh et,  убеждаемся в их линейной зависимости:  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  задает шаг дискретизации аргумента.  Используя формулы тригонометрии, можем записать

 .   (1.24)

Домножим эти равенства на коэффициенты, указанные справа, и сложим:

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

Это – искомое разностное уравнение второго порядка. Начальные условия для него находим из соотношений (1.24): y(0) = 0;   y(h) = sin h.

При  h = 0,1  получаем

Проверка. Чтобы убедиться в правильности построенного уравнения, решим его при начальных условиях  y(0) = 0,   y(h) = sinh.  

Выписываем характеристический полином и находим его корни:

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

Подставляя начальные условия, находим, что  с1 = 0,  с2 = 1,  т.е. при  t = kh  полученное решение совпадает с требуемым.

Пример 3. Найдем методом последовательных сдвигов определяющее разностное уравнение для функции  y(t) = tet.

Решение. Выберем шаг  h  и перейдем к дискретному времени  t = kh:

Раскрывая скобки, выражение для сдвинутой на такт функции  yk+1  можно переписать в форме

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

Это нестационарное уравнение, так как коэффициент  а  зависит от номера такта.

Чтобы получить стационарное уравнение, выпишем дополнительно функцию  :

Cовокупность трех функций  kekh,   (k+1)ekh,   (k+2)ekh  линейно зависима, поэтому дальнейших сдвигов не требуется.

Сложим равенства

домножив их на коэффициенты, указанные справа.

При этом правая часть обращается в нуль, и мы получаем искомое разностное уравнение

В частности, при  h = 0,1  оно принимает вид

Метод характеристического полинома. При использовании этого метода требуется знать соответствие между корнями характеристического уравнения и решениями соответствующего разностного уравнения (табл. 1).

Таблица 1

Вид корня

Характер решения

1

с

2

3

4

5

6

]

По заданному виду функций определяются необходимые корни  z1, …, zn  и составляется характеристическое уравнение  (zz1)(zz2) … (zzn) = 0.  От него легко перейти к разностному уравнению

y(t+nh) + an-1 y(t+(n-1)h)+ … +a1 y(t+h) +a0 y(t) =0,

где коэффициенты  ai  рассчитываются по известным корням  zi  с помощью формул Виета  

Пример 4. Найдем определяющее разностное уравнение для функции  y(t) = tet  методом характеристического полинома.

Решение. Перейдем к дискретному времени  t = kh,  h = 0, 1, 2, … .  Нам требуется получить функцию вида    следовательно характеристический полином должен иметь кратный корень    Восстанавливая полином по корням, получаем

Ему соответствует  разностное уравнение второго порядка

c начальными условиями   

При  h = 0,1  оно принимает вид  

Пример 5. Пусть необходимо воспроизвести функцию  .  

Решение. При  t = kh  первому слагаемому соответствуют корни  , а второму – корни   Коэффициенты  и    находятся из соотношений    откуда  

Характеристическое уравнение имеет вид

.

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

Следовательно, искомое разностное уравнение имеет вид

Начальные условия равны значениям функции  y(t)  при  t = 0, h, 2h, 3h.

Метод списков. В некоторых случаях требуется воспроизвести не одну, а несколько функций времени  y1 = f1(t), …, ym = fm(t). Для этого можно использовать описанные выше способы, что приводит к получению нескольких определяющих разностных уравнений разных порядков. Однако при этом не гарантируется минимальность модели в смысле числа используемых элементов задержки.

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

,    (1.25)

где Х – вектор переменных состояния,  Y – вектор выходных сигналов, реализующих заданные функции времени, Х0 – вектор начальных условий,  А  и  С – постоянные матрицы, подлежащие определению.

Алгоритм получения системы (1.25) минимального порядка, обеспечивающей выполнение равенств  y1 = f1(t), …, ym = fm(t),  состоит в следующем.

Шаг 1. Составляем список  S1  линейно независимых элементарных функций    через которые линейно выражаются заданные функции  f1(t), …, fm(t).  В частности, можно принять 1 = f1(t), …, m = fm(t).

Шаг 2. Сдвигаем функции списка S1  на такт  h,  раскрываем выражения  1 (t+h), …,
q = (t+h)  и составляем список  S2  новых функций   появившихся после такого сдвига и не выражающихся линейно через функции списка  S1.

Шаг 3. Сдвигаем функции списка S2  на такт  h,  раскрываем выражения  q+1 (t+h), …,
r = (t+h)и составляем список  S3  новых функций   появившихся после этого сдвига и не выражающихся линейно через функции списков  S1  и  S2.

Шаги, аналогичные шагам 2 и 3, повторяются до тех пор, пока окажется, что после очередного сдвига на  k-м шаге новых функций, линейно независимых от предыдущих, не появилось. Если такого шага не наступит, то это означает, что системы определяющих разностных уравнений вида (1.25) для заданных функций не существует.

Шаг k+1. Размерность системы (1.25) берется равной общему числу  n  функций    вошедших в списки  S1,  S2, …, Sk.  Все эти функции линейно независимы, чем гарантируется минимальное число уравнений в системе. Cами эти функции принимаются за переменные состояния системы (1.25)

x1=1(t), x2=2(t), …, xn=n(t).     (1.26)

Записывая равенства (1.26) для момента времени  t+h  и выражая их правые части через  x1(t), …, xn(t)  (это возможно, т.к. они должны линейно зависеть от  ), получим систему  разностных уравнений    Матрицу  С  находим, выражая функции fi(t) через  .  Компоненты вектора начальных условий  Х0  получаем из равенств (1.26) при  t = 0.

Проиллюстрируем описанный алгоритм на примере.

Пример 6. Пусть требуется найти определяющие разностные уравнения в форме (1.25) для функций

представляющих параметрическое описание спирали Архимеда.

Решение.

Шаг 1. Принимая    получаем список  S1:

Шаг 2. Выполняем сдвиг функций    на такт

.

Новыми функциями здесь являются  sint  и  cost,  поэтому

Шаг 3. Выполняем сдвиг функций    на такт

,    .

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

Шаг 4. Размерность системы (1.25) берем равной 4, равенства (1.26) имеют вид

После сдвига на такт, с учетом формул, полученных на предыдущих шагах, приходим к следующей системе определяющих разностных уравнений

Ей соответствует описание     с матрицами

Моделируя его в любом математическом пакете и строя зависимость y1=f(y2), получим на дисплее график спирали Архимеда.

1.5. Задачи и упражнения

  1.  Показать эквивалентность уравнений (1.5) и (1.6).

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

константа  x(t) = a ; линейная функция x(t) = at + b ;

парабола   x(t) = at2 ;  косинусоида  x(t) = соst.

Аргумент везде изменяется дискретно  t = 0, 1, 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. Найти общее решение разностного уравнения   

Ответ.    

6. Найти общее решение разностного уравнения   

Ответ.    

7. Найти решение разностного уравнения

Ответ.    

8. Найти решение разностного уравнения

Ответ.    

9. Найдите решение системы разностных уравнений второго порядка

x(t+1) = x(t) – 3y(t)

y(t+1) = 3x(t) + y(t).

Указание. Собственные числа и собственные векторы матрицы  A в данном случае имеют вид:  

,   H1=,   H2 =.

10. Найти частное решение неоднородного уравнения   

Ответ.   

11. Найти общее решение неоднородного  уравнения  

Ответ.    

12. Найти общее решение неоднородного уравнения   

Ответ.    

13. Записать аналитические выражения для переменных x, y, z, v в задаче Нарайаны.

14. Довести до конца решение задачи об электрической цепи (пример 4, разд. 1.1 и пример 2, разд. 1.2), полагая  x0 = 100В,  n = 14.  

Ответ:  мВ, где n – числа Фибоначчи: 0=1=1, 2=2, 3=3, 4=5, …

15. Популяция жуков [14]. В популяцию входит три вида жуков: синие, красные и зеленые. В дискретные моменты времени  t=1, 2, 3, … вместо каждого жука появляется три новых, причем:

  •  каждый синий жук делится на двух красных и зеленого;
  •  каждый красный жук делится на синего и двух зеленых;
  •  каждый зеленый жук делится на трех красных.

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

Указание. Решение сводится к исследованию системы  разностных уравнений

Главное собственное число матрицы =3. Ему отвечает собственный вектор  H=[ 3  9  7  ]Т.

Ответ. Динамика относительной численности жуков сначала носит колебательный характер (рис. 1.5). После окончания переходного процесса численности синих, красных и зеленых жуков будут относиться как 3:9:7. 

  1.  ПРИМЕНЕНИЕ ЖОРДАНОВЫХ ЦЕПОЧЕК ВЕКТОРОВ ПРИ РЕШЕНИИ

ДИФФЕРЕНЦИАЛЬНЫХ И РАЗНОСТНЫХ УРАВНЕНИЙ

2.1. Корневые векторы и циклические пространства

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

Процедуры вычисления собственных чисел и собственных векторов входят в состав всех профессиональных математических пакетов, таких как MATHCAD, MATLAB, MAPLE и т.д.

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

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

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

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

где   – собственные числа и собственные векторы матрицы  А,  справедлива только для матриц  А  с попарно различными собственными числами.

Возникает вопрос об обобщении этой формулы на случай кратных собственных чисел. Очевидно, что оно потребует введения корневых векторов и построения жордановых цепочек векторов, зависящих от начальных условий системы. Аналогичный вопрос относится к системам разностных уравнений вида    и формулам (1.17), (1.18).

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

Минимальный многочлен вектора. Пусть дана матрица  А.  Возьмем произвольный вектор  b  и составим ряд векторов  b,  Ab,  A2b, … . Пусть  m – наименьшее число, при котором очередной вектор будет линейно зависим от предыдущих

   (2.1)

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

     (2.2)

Всякий многочлен  ,  для которого имеет место равенство (2.2), называется аннулирующим для вектора  b  (относительно оператора  А). Из всех аннулирующих многочленов для вектора  b  построенный нами многочлен имеет наименьшую степень со старшим коэффициентом 1. Он называется минимальным аннулирующим многочленом вектора  b  или просто минимальным многочленом вектора  b.

В теории управления матрица  известна, как матрица управляемости системы    Число  m  совпадает с ее рангом  rank  R = m  и называется индексом управляемости системы. Если  rank R < n,  где  n – размер матрицы  А,  то система будет неуправляемой – это известный критерий неуправляемости Калмана.

С помощью произвольного входного сигнала  u(t)  такую систему можно вывести из состояния покоя, но она будет двигаться только в  m-мерном подпространстве  n-мерного пространства состояний. Оно называется подпространством управляемости.

Знание минимального многочлена вектора  b  позволяет сформировать такой входной сигнал  u(t)  конечной длительности, реакция на который после его окончания будет равна нулю, т.е. он аннулируется системой. Простой пример – последовательность из  m + 1  прямоугольных импульсов, амплитуды которых равны коэффициентам аннулирующего многочлена. На этом принципе, в частности, основан один из методов тестового диагностирования систем автоматического управления  (метод Шрайбера или метод комплементарного сигнала1).

Циклические пространства. Вновь рассмотрим векторы    где  m – порядок минимального многочлена вектора  b.  Они линейно независимы и образуют базис  m-мерного пространства  Rm,  которое называется циклическим.

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

Минимальный многочлен первого базисного вектора одновременно будет минимальным многочленом пространства  Rm .  Отсюда вытекает следующий критерий цикличности пространства.

Теорема. Пространство циклично относительно оператора  А  тогда и только тогда, когда число его измерений совпадает со степенью его минимального многочлена.

Например, минимальный многочлен единичной матрицы имеет вид    Здесь  m = 1,  т.е. критерий цикличности не выполняется даже для двумерной плоскости.

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

2.2. Жордановы цепочки векторов

Пусть у матрицы  А  размера    имеется собственное число  1  кратности  k  и ему отвечает циклическое подпространство размерности  k  с порождающим вектором  b.  Для этого вектора многочлен  () = (1)k  будет минимальным многочленом.

Рассмотрим векторы

b1 = (A1E)k-1b,     b2 = (A1E)k-2b, …, bk = b.    

Умножая их на  A1E  и учитывая, что  (A1E)kb = 0,  получаем цепочку равенств:

(A – 1E)b1 = 0,     (A – 1E)b2 = = b1        , …,    (A – 1E)bk = bk-1   

или

A b1= 1b1,     A b2= 1b2 + b1        , …,    A bk= 1bk + bk-1.   (2.3)

Линейно независимые векторы  b1,  b2, …, bk = b  образуют цепочку векторов, которая называется жордановой. В совокупности все жордановы цепочки матрицы  А  образуют жорданов базис в  Rn.

Корневые векторы. Векторы жордановой цепочки принадлежат множеству так называемых корневых векторов.

Определение ([7], с.147). Вектор  Н  называется корневым вектором матрицы А,  принадлежащим ее собственному числу  ,  если

Наименьшее число  k,  при котором выполняется это равенство, называется высотой корневого вектора.

Заметим, что при  k = 1  получаем обычное определение собственного вектора матрицы. Таким образом, понятие корневого вектора обобщает понятие собственного вектора в случае кратных собственных чисел.

Векторы жордановой цепочки – это совокупность корневых векторов различной высоты (от 1 до  k  включительно), получаемых на основе одного порождающего вектора путем его последовательного умножения на матрицу  

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

Перечислим некоторые свойства корневых векторов.

  1.  Корневые векторы высоты  k  образуют линейное пространство размерности  k – инвариантное подпространство матрицы  А.
  2.  Умножение корневого вектора высоты  k > 1  на матрицу    уменьшает его высоту на 1.
  3.  Следующие операции не меняют высоты корневого вектора:
  •  умножение на число;
  •  умножение на матрицу  А;
  •  добавление линейной комбинации корневых векторов меньшей высоты.

Отметим, что в некоторых книгах по линейной алгебре вместо названия "корневой вектор высоты  k"  используется термин "обобщенный собственный вектор ранга  k".

Опишем еще один подход к определению корневых векторов  [4].

Пусть у матрицы  А  собственному числу    кратности два отвечает единственный собственный вектор  Н,  так что  АН = Н.  Предположим, что в результате малой вариации матрицы  А  получается матрица  А1  с близкими собственными числами    и   + ,  которым отвечают близкие собственные векторы  Н  и  Н + G  (рис. 2.1). Здесь    и   – величины одного порядка малости. Представляется разумным в качестве обобщенного собственного вектора матрицы  А  взять вектор  G,  вернее тот вектор, к которому он стремится в пределе при  

Чтобы получить уравнение для этого вектора, рассмотрим равенство

A1(H + G) = ( + )(Н + G).

Раскрывая скобки, пренебрегая произведением малых величин и учитывая равенство  A1H = Н,  получим

В пределе при    имеем  ,  так как величины    и    одного порядка малости. Таким образом, вектор  G  является ненулевым решением уравнения

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

Пример 1. Найдем жорданову цепочку векторов матрицы

Она имеет единственный собственный вектор  Н1 =[1  0  0]Т,  отвечающий собственному числу    = 2  кратности k=3. Следовательно, жорданова цепочка будет состоять из трех векторов. Корневые векторы находим из уравнений  А0Н2 = Н1,  А0Н3 = Н2,  где

 

Здесь  с2  и  с3 – произвольные постоянные.

Если собственный вектор  Н1  умножить на произвольный коэффициент  с10,  то проводя аналогичные выкладки, получим

Варьируя коэффициенты  с1, с2 , с3  можно получать различные жордановы цепочки векторов. В частности, при  с1 = 9,  с2 = с3 = 0  имеем

При решении практических задач жорданову цепочку бывает удобнее строить не с начала, а с конца, выбирая в качестве  Н3  конкретный корневой вектор высоты 3. Тогда неопределенность полностью устраняется и векторы  Н2 , Н1  находятся однозначно.

Например, взяв все компоненты вектора  Н3  единичными, получим

Во всех случаях матрица корневых векторов  H = [H1  H2  H3]  оказывается треугольной, что следует из треугольности матрицы  А.

2.3. Дифференциальные уравнения и жордановы цепочки векторов

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

,     (2.4)

где А – постоянная матрица размера  ,  X – вектор переменных,

наряду с формулой  используется формула [5, 13]

   (2.5)

где  i  и  Hi – собственные числа и собственные векторы матрицы  А,  

сi – произвольные постоянные, зависящие от начальных условий.

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

Пусть   – собственное число матрицы  А  кратности  k.  Будем искать соответствующие ему компоненты решения в форме экспоненты  ,  домноженной на полином от  t  степени  k–1:

  (2.6)

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

Подставим это выражение в систему (2.4). После дифференцирования и сокращения на множитель    получим:

Приравняем коэффициенты при одинаковых степенях  t  в правой и левой частях уравнения:

(2.7)

Мы получили уравнения, которые совпадают, с точностью до обозначений, с цепочкой уравнений (2.3). Следовательно, они описывают жорданову цепочку векторов, отвечающих кратному собственному числу матрицы  А.  При этом  Hk-1 – это собственный вектор матрицы  А,  а  Н0 – корневой вектор высоты  k.

Для определения вектора  Н0  используются начальные условия системы (2.4). Особенно просто это делается в случае, если у матрицы  А  кроме собственного числа    кратности  
k = n  нет других собственных чисел (матрица с одной жордановой клеткой). Тогда формула (2.6) описывает общее решение системы, а корневой вектор  Н0  совпадает с вектором начальных условий  Н0 = Х0,  что прямо следует из равенства (2.6) при  t = 0.

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

Требуется найти формулы, описывающие свободное движение этой системы из начальных условий  х1 (0) = а,   х2 (0) = b,   х3 (0) = c.

Решение. Перейдем к описанию в пространстве состояний вида (2.4):

    (2.8)

Матрица  А  этой системы такая же, как в предыдущем примере, она имеет единственный собственный вектор и одно собственное число   = 2  кратности три. Следовательно, решение системы может быть записано в форме (2.6) при  k = 3:

    (2.9)

Полагая  t = 0,  устанавливаем, что  H0 = X0.  Этот вектор порождает жорданову цепочку векторов  Н0,  Н1, Н2,  связанных формулами (2.7):

Подставляя эти значения в формулу (2.9), получаем окончательный результат:

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

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

Выделение модальных компонент. Одно из эффективных применений формулы (2.5) связано с определением начальных условий, позволяющих подавлять в решении все модальные компоненты, кроме одной. Это достигается, если в качестве начальных условий брать один из собственных векторов матрицы  А.  Например, полагая  X0 = Hi,  получим    т.е. в решении будет присутствовать только одна модальная компонента.

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

Поясним эту ситуацию на примере матрицы  А  с единственной жордановой клеткой.

Пример 2. Рассмотрим схему из трех последовательно соединенных апериодических звеньев, показанную на рис. 2.3.

Ее описание в пространстве состояний имеет вид

Матрица  А  обладает единственным собственным числом   = – а  кратности три.

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

Требуется указать начальные условия  x0,  y0,  z0,  при которых в этой формуле будет оставаться только одно из слагаемых (первое, второе или третье соответственно).

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

Действительно, эти векторы удовлетворяют равенствам

где   = – а;     т.е. представляют собой  жорданову цепочку с порождающим вектором  X3 (0) = [0  0  1]T.

Находя выходной сигнал для каждого из трех начальных условий

убеждаемся, что в каждом случае он содержит только одно слагаемое.

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

Если в качестве порождающего вектора взять  X3 (0) = [1  1  1]T,  то картина несколько изменится:

При этом выходные сигналы для начальных условий  Х1(0),  Х2(0),  Х3(0)  будут такими:

Видно, что в первом из них отсутствуют две модальные компоненты, а во втором – одна.

В принципе в качестве порождающего вектора жордановой цепочки можно брать любой корневой вектор высоты три – им является всякий вектор с ненулевой третьей компонентой. Однако только вектор  X(0) = [0  0  1]T  удовлетворяет условию задачи. Он порождает цепочку ортогональных векторов, имеющих особенно простой вид и полезные свойства.

Определение жордановой цепочки векторов по начальным условиям. Для того, чтобы применять формулу (2.6) при решении дифференциальных уравнений, нужно уметь определеять корневые векторы  Н0, …, Нk-1,  которые зависят от начальных условий. Выше было показано, что для систем с единственным собственным числом порождающий корневой вектор равен вектору начальных условий  Н0 = Х0,  а остальные векторы  Hi  находятся по формулам (2.7).

В случае, когда у системы есть и простые, и кратные собственные числа, процедура определения жордановой цепочки векторов, входящих в формулу (2.6), усложняется. Теперь вектор  Н0  и другие векторы цепочки будут зависеть не только от начальных условий системы (2.4), но и от всей совокупности собственных векторов матрицы  А.

Опишем эту процедуру подробнее, приняв для определенности, что у матрицы  А  есть собственное число  0  кратности  k   и  nk  простых собственных чисел  1 , …, n-k .  Тогда формула для общего решения системы (2.4) получается объединением формул (2.5) и (2.6):

 (2.10)

Здесь  Gi – собственные векторы, принадлежащие простым собственным числам,

 Hi – корневые векторы, принадлежащие собственному числу 0, причем  Нk-1 – единственный собственный вектор, принадлежащий этому собственному числу.

Постоянные  ci  и векторы  Hi  зависят от начальных условий системы и подлежат определению. Основой для этого служит соотношение

    (2.11)

получаемое из формулы (2.10) при  t = 0.

Изложим процедуру отыскания решения (2.10) системы дифференциальных уравнений (2.4) в виде алгоритма, содержащего 4 шага.

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

Шаг 1. Определяются собственные векторы    матрицы  А,  принадлежащие простым собственным числам.

Шаг 2. Определяются  постоянные  .  Для этого обе части соотношения (2.11) умножают на матрицу  ,  где    и решают полученную систему уравнений:

Вводя обозначения    ее можно переписать в компактной форме

Число  n  уравнений в этой системе превышает число  nk  неизвестных, однако количество линейно независимых уравнений равно  nk,  поэтому система имеет решение, и притом только одно.

Шаг 3. Из соотношения (2.11) определяется неизвестный вектор  Н0 :

Шаг 4.  На основе порождающего вектора  Н0  строится жорданова цепочка векторов

Процедура отыскания решения системы уравнений (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.  

Решение. 

Шаг 1. Описание системы в пространстве состояний характеризуется матрицей

.

Она имеет собственный вектор  G1,  отвечающий простому собственному числу  1 = – 2,  и собственный вектор  Н1,  отвечающий кратному собственному числу   = – 1:

Следовательно, решение ищем в форме

,

где постоянную  с1  и векторы  Н0 ,   Н1  считаем пока неизвестными.

Шаг 2. При  t = 0  имеем равенство    

Умножая его на матрицу    где  ,  

получаем уравнение    

Для варианта  а) оно имеет вид:

откуда  с1 = 1.

Шаг 3. Определяем вектор  

Шаг 4. Определяем вектор    Видим, что он с точностью до постоянного коэффициента совпадает с собственным вектором, определенным на первом шаге.

Подставляя найденные значения в исходную формулу, получаем искомое решение

или в скалярной записи

В случае б) аналогичные выкладки приводят к следующему результату

.

Его скалярная запись имеет вид:

2.4. Разностные уравнения и цепочки корневых векторов

Разностные уравнения с кратными корнями. Перейдем к решению системы разностных уравнений

или, в другой системе обозначений,

.     (2.12)

Здесь  А – постоянная матрица размера  nxn,  X – вектор переменных,  t – дискретное время. Наиболее часто используется показательная форма записи решения  однако она не всегда удобна, например, в тех случаях, когда надо получить аналитические формулы для отдельных компонент вектора  Х.  

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

    (2.13)

где  i , Hi  – собственные числа и собственные векторы матрицы  А,  

ci  – произвольные коэффициенты.

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

Пусть матрица  А  имеет собственное число  0  кратности  m.  Будем искать соответствующие ему компоненты решения в виде

    (2.14)

После подстановки этого выражения в уравнение (2.12) и сокращения на множитель  k,  получим

Раскроем скобки и приравняем коэффициенты при одинаковых степенях  k  справа и слева:

 (2.15)

Вводя обозначение    перепишем эти равенства в форме

Следовательно,  – корневые векторы матрицы  А  высоты  1, 2, …, m  соответственно. Обозначив  H = [H0, …, Hm-1],  запишем равенства (2.15) в матричном виде

где         (2.16)

Треугольная матрица  В  имеет размер  mxm,  в ее строках стоят биномиальные коэффициенты. Вычитая из обеих частей равенства (2.16) матрицу  0Н,  получим эквивалентное матричное равенство

где       (2.17)

Присоединив к нему равенство  Н0 = Х0,  которое следует из выражения (2.14) при  t = 0,  получим систему линейных алгебраических уравнений относительно неизвестных векторов  

Опишем процедуру их определения, не требующую обращения матриц.

Многократно умножая обе части равенства (2.17) слева на матрицу  А0,  получим совокупность соотношений

  (2.18)

Обозначим через    промасштабированные векторы жордановой цепочки, построенные на основе порождающего вектора  Н0 ,  а через  l = [0  1  …  1]Tпервый столбец матрицы  В0.

Приравняем первые столбцы в каждом из матричных равенств (2.18):

  или в матричной форме    h = HB1,

где  

Отбросив в матрице  В1  первую строку (она нулевая), а в матрице  Н – первый столбец (вектор  Н0),  получим систему уравнений

    (2.19)

где  В2 –  нижне треугольная матрица размера  (m–1)x(m–1).

Из системы (2.19) легко последовательно определить все векторы  Hi,  начиная с вектора  Hm-1 (он с точностью до масштабного множителя будет совпадать с собственным вектором матрицы  А).  Тем самым векторы  Hi  окажутся выраженными через векторы  hi  с помощью  формулы вида

   (2.20)

Здесь   D – нижне треугольная матрица, ее элементы постоянны и не зависят ни от матрицы  А,  ни от начальных условий системы.

Приведем матрицы  D  для кратности собственных чисел от  m = 2  до  m = 5:

Видно, что каждая следующая матрица получается окаймлением предыдущей. Элементы первого столбца матрицы  D  описываются формулой  ,  а диагональные элементы – формулой    Это позволяет выписать в явном виде выражения для векторов  Н1  и  
Hm-1  при произвольном  m:

  (2.21)

Пример 1.  Найдем решение системы разностных уравнений

В данном случае  1 = 2 = 3 = 2,   m = 3,  поэтому решение будет иметь вид

Находим векторы  Н0,  h1,  h2:

Векторы  Н1  и  Н2  вычисляем по формулам (2.21):

Следовательно, решение системы разностных уравнений имеет вид:

Определение цепочки корневых векторов по начальным условиям. Корневые векторы  Н1, …,  Hm-1,   входящие в формулу (2.20), однозначно определяются вектором  Н0,  порождающим жорданову цепочку  h1, …, hm-1.  Вектор  Н0  в  свою очередь, зависит от начальных условий системы. Если у матрицы  А  системы  разностных уравнений (2.20) нет других собственных чисел, кроме  0,  то эта зависимость крайне проста:  Н0 = Х0.

В более общем случае, когда у матрицы  А  кроме собственного числа  0  кратности m есть простые собственные числа 1, …, n-m,  эта зависимость становится более сложной. Формула для решения системы (2.12) принимает вид

 (2.22)

где  Gi – собственные векторы матрицы  А,  отвечающие простым собственным числам,  

Hi – корневые векторы матрицы  А,  принадлежащие собственному числу  0.

Постоянные  сi  и векторы  Нi  зависят от начальных условий и подлежат определению.

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

   (2.23)

и исключение вектора  Н0  путем умножения на матрицу  :  

  (2.24)

Из этой системы линейных уравнений определяются коэффициенты  сi,  после чего из соотношения (2.23) находится вектор  Н0.  Далее на его основе строится жорданова цепочка векторов  h1, …, hm-1  и из системы уравнений (2.19) находят векторы  Нi.

Опишем эту процедуру в виде алгоритма, содержащего 5 шагов.

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

Шаг 1. Определяются собственные векторы    матрицы  А,  принадлежащие ее простым собственным числам.

Шаг 2. Из системы (2.24) определяются  постоянные  .  Вводя обозначения    ее можно записать в более компактной форме

Число  n  уравнений в этой системе на  m  превышает число неизвестных из-за того, что часть уравнений является линейно зависимой (и при  решении может не учитываться).  

Шаг 3. Из соотношения (2.23) определяется вектор  Н0:

Шаг 4. На основе порождающего вектора  Н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). Тогда динамику популяции можно описать системой трех разностных уравнений:

Переходя к матричной форме записи, получаем

Примем для определенности  а1 = 0,  а2 = 9,  а3 = 12,  р1 = 1/3,  р2 = 1/2.

Здесь компоненты вектора  Х  отражают возрастную структуру популяции, в которой выделено три возрастных группы. Числа 9 и 12 характеризуют среднюю плодовитость поколений (младшая группа не дает приплода), величины 1/3  и  1/2  – доля рыб, выживающих в младшей и средней группе и переходящих в следующую возрастную группу.

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

,

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

Первая строка характеризует динамику изменения численности мальков, вторая и третья – рыб среднего и старшего возрастов. Из приведенного массива трудно усмотреть закон изменения численности популяции.

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

Шаг 1. Для определения собственных чисел матрицы  А  выпишем ее характеристическое уравнение

Матрица  А  системы имеет собственное число  1 = 2,  которому отвечает собственный вектор  G1 = [2 4  4  1]T,  и собственное число  0 = –1  кратности два. Следовательно, решение должно иметь вид:

Шаг 2. Постоянную  с1  находим из соотношения (2.24)  .  Подставляя в него значения

получаем  с1 = 1/9.

Шаг 3. Вектор  Н0  определяем из соотношения (2.23):

Шаг 4. Определяем вектор  h1:

Шаг 5. Вектор  Н1  находим из системы (2.20). Учитывая, что в нашем случае  D2 = 1,  получаем  Н1 = h1.

Следовательно, решение имеет вид:

или в скалярной форме записи

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

Графики относительной численности рыб каждого возраста, определяемые формулами    приведены на рис. 2.5. Их установившиеся значения относятся, как  24:4:1.

Рис. 2.5

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

Ответ прост: для поддержания популяции на постоянном уровне главное собственное число    матрицы  А  должно быть равно 1, если же оно больше единицы, то доля изымаемой части должна составлять  100 ( – 1)/ %.  Например, при   = 2  следует регулярно изымать половину популяции, а при   = 1,1 – около 9%.

2.5. Задачи и упражнения

1. Построить жорданову цепочку векторов для матриц

взяв в качестве порождающего вектора  b  вектор с единичными компонентами.

Указание.  Собственные числа матрицы  А1:  1 = 2 = 7;  матрицы  А2:  1 = 2 = 3 = 0;  матрицы  А3:  1 = 2 = 3 = 4 = 2.

2. Построить корневые подпространства для матриц

Ответ.    А1:  базис корневого подпространства для   = 0 –  вектор  [0  1  –1]T,  для   = 1  векторы  [0  1  1]T,  [0  1  0]T.

А2:  базис корневого подпространства для   = –1 –  векторы  [1  0  0]T,  [0  1  0]T,  [0  0  1]T.

А3: базис корневого подпространства для   = 2 – векторы  [2  –1  0  0]T,  [1 0 1 0]T, [2 0 0 1]T.

А4:  базис корневого подпространства для   = –2 –  вектор  [0  1  0  –1].

3. Обобщить формулу (2.10) на случай двух кратных корней. Как в этом случае получить систему линейных уравнений для определения постоянных  сi?

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

Указание. Базис корневых подпространств для   = –1 –  векторы  [1  1  0  0]T,  
[0  0  1  1]
T. Базис корневых подпространств для   = 1 –  векторы  [3  1  0  0]T,  
[0  -2  3  1]
T.

5. Если порядок минимального полинома матрицы  А  меньше порядка ее характеристического полинома, то один из корней  1, …, n-m  может совпадать с корнем  0,.  Как будет выглядеть формула (2.10) в этом случае?

6. Решить систему дифференциальных уравнений    С матрицей  ,  у которой имеются две жордановы клетки, отвечающие собственному числу  0 = 3.  Принять  Х0 =[a  b  c]T.

7. Обобщить формулу (2.12) на случай двух кратных корней. Как в этом случае получить систему линейных уравнений для определения постоянных  сi?

8. Решить систему разностных уравнений  Xk+1 = AXk  с матрицей  .

Указание. Все собственные числа матрицы  А  равны  единице, им отвечают две жордановы клетки второго порядка.

9. Если порядок минимального полинома матрицы  А  меньше порядка ее характеристического полинома, то один из корней  1, …, n-m  может совпадать с корнем  0,.  Как будет выглядеть формула (2.24) в этом случае?

10. Решить систему разностных уравнений Xk+1 = AXk с матрицей  .

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

  1.  РАЗНОСТНЫЕ УРАВНЕНИЯ И Z-ПРЕОБРАЗОВАНИЕ

3.1. Определение  z-преобразования

При исследовании непрерывных динамических систем с дискретным временем (они называются также импульсными системами) и при решении разностных уравнений широко применяется математический аппарат z-преобразования (другие названия: дискретное преобразование Лапласа, преобразование Лорана).

Определение. Пусть f(t) – функция целочисленного аргумента t = 0, 1, 2 ... .
Ее
z-преобразованием называется функция комплексного аргумента z, определяемая формулой

   (3.1)

Условная запись z-преобразования: .

Свойства z-преобразования

  1.  Линейность  
  2.  Сдвиг по времени  

.

3. Умножение на время  .

Приведем примеры z-преобразований простейших функций:  

f(t)  =          F(z)=1

При выводе последних формул удобно использовать теорему о дифференцировании изображений (свойство 3). Она доказывается через дифференцирование соотношения  (3.1) по переменной z.

Заметим, что в тулбоксе Symbolic пакета MATLAB имеются команды ztrans и
i
ztrans для выполнения прямого и обратного z–преобразований. Например, в результате выполнения команд  syms t, ztrans(t^2)   получим  ans =z*(z+1)/(z–1)^3.

Таблица  z-преобразований

Оригинал   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)

z (F(z)–f0)

f(t+2)

z2(F(z)–f0z–1f1)

3.2. Решение разностных уравнений с помощью z-преобразования

Пусть дано  разностное уравнение  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-преобразования содержит 4 шага.

Шаг 1. Применить z-преобразование к исходному  разностному уравнению, заменяя  y(t)  на  Y(z);  y(t+1)  на  z(Y(z)–y0)  и т.д.

Шаг 2. Из полученного алгебраического уравнения выразить  Y(z).

Шаг 3. Выполнить разложение на простые дроби.

Шаг 4. Пользуясь таблицей, выполнить обратное z-преобразование.

Рассмотрим примеры применения этого алгоритма.

Пример 1. Требуется решить разностное уравнение первого порядка

y(t+1)–0,9y(t)=0,1;      y(0)=5.       (3.2)

Шаг 1.    Применяем  z-преобразование:  z(Y(z)–5)–0,9Y(z)=.

Шаг 2.    Выражаем  Y(z):    .

Шаг 3.     Раскладываем на простые дроби

.

Шаг 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.

График полученного решения показан на рис. 3.1.

Пример 2.  Пусть дано  разностное уравнение второго порядка

yn+2 +a1 yn+1+ a0 yn=un;

с заданными начальными значениями  y0  и y1.

Шаг 1. Применяем к нему z-преобразование

Шаг 2. Выражаем Y(z):

Шаг 3. Чтобы найти  yn,  надо разложить правую часть на простые дроби.

Положим  а1 =–1;  а0=0,5;  un=0;   y0=0;  y1=1.  Тогда получаем

.

Корни знаменателя – комплексные, следовательно, эта дробь – простая.

Шаг 4. С помощью таблицы преобразований находим оригинал  yn=can sinwn.  Постоянные с,  а,  w  определяем из соотношений:    откуда  

Окончательный вид решения:    

Пример 3. Задача о банковском вкладе. Мы открываем счет в банке, наш начальный взнос составляет  В  долларов. Каждый год добавляем по b долларов, годовой процент равен  k. Какая сумма окажется на нашем счете через  n  лет? Какова будет прибыль?

Обозначая искомую величину через  yn,  можно записать следующее разностное уравнение:

yn+1=yn(1+k)+b;    y0=B.

Заметим, что этому уравнению  можно сопоставить схему моделирования, содержащую один элемент задержки ЭЗ (рис. 3.2,а)

а)       б)

   

     

Рис. 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)nb/k.

Неизвестный коэффициент  С  находим из начальных условий

y0=B=Cb/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 = ynxn,    где  xn – общий взнос. Она определяется формулой

zn = yn(B+nb).  

Пусть  B =100$,   b =50$,    k =5%,    n =10 лет.  Подставим эти данные в полученное решение:

.

При  n=10  получаем:   y10 =791,78,   z10 =y10 –600=191,78.

Таким образом, всего за 10 лет мы внесли 600$, получили в итоге  792$,  прибыль составила 192$.

Второй способ (z–преобразование)

Исходное разностное уравнение имеет вид   yn+1 =(k+1)yn +b;    y0 =B.

Применяя к нему z–преобразование  

получаем алгебраическое уравнение:

.

Выражаем из него Y(z):

Раскладываем это выражение на две простые дроби

.

Возвращаемся к оригиналам

.

Ответ получили тот же, что и первым способом.

Пример 4.  Уравнение второго порядка с вещественными корнями

.

Находим изображение

В данном случае имеет место сокращение нуля и полюса.

Рассмотрим еще несколько примеров.

Пример 5. Уравнение третьего порядка с кратными корнями

 yn+3 – 5yn+2 + 8yn+1 – 4yn= 0;           y0= 0,   y1= 2,   y2= 1.  

Находим изображение и выполняем разложение на простые дроби

.

Выполняем обратное  z-преобразование

.

Пример 6. Найдем решение задачи о электрической цепи (см. рис. 1.1) с помощью z-преобразования. Примем для определенности  R=2r,  u0=1.  Эта цепь описывается  разностным уравнением (1.7), которое при R=2r  принимает вид

Применяем к нему z–преобразование:

Отсюда, учитывая что  u0 =1,  находим

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

Переходя к оригиналам, получаем формулу для расчета напряжения  в произвольной точке  R-цепи:

uk  = ch k + C sh k .     (3.4)

Константа    находится из условия  ch =2,  откуда    1,32;  С – постоянная, зависящая от  u1.  

Проблема заключается в том, что значение u1  неизвестно и для определения постоянной  С  нужно дополнительное условие. Обойти эту проблему удается, рассмотрев последнее звено схемы рис. 1.1 (при  R = 2r  оно представляет собой делитель напряжения с коэффициентом  1/3) и выписав краевое условие  un = un-1 /3.

С учетом этого можно записать

3(ch n + C sh n) =ch (n – 1) + C sh (n – 1),

откуда находим  

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

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)  и  sh (n – 1)  в числителе и знаменателе и учитывая, что ch =2,    получаем  

Подставляя это выражение в формулу (3.4) и выполняя преобразования, получаем окончательный результат:

В частности для  k = 0  (вход схемы) имеем  u0 = 1,  а для  k = n  (выход схемы) получим

Если схема содержит только одно звено (простой делитель напряжения), то

3.3. Цифровые интеграторы и дискретные передаточные функции

Для описания цифровых линейных систем используются разностные уравнения различных порядков и дискретные передаточные функции Q(z).

Определение. Дискретной передаточной функцией линейной системы с входом  x(t)  и выходом  y(t)  называется отношение z-преобразования выходного сигнала к z-преобразованию входного сигнала при нулевых начальных условиях  Q(z) = Y(z)/X(z).

Аналоговый интегратор описывается дифференциальным уравнением первого порядка   и передаточной функцией Q(p)=1/p.  Простейший цифровой интегратор описывается  разностным уравнением первого порядка, которое реализует численное интегрирование по методу Эйлера (методу прямоугольников). Несколько большую точность обеспечивает цифровой интегратор, использующий формулу трапеций.

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

Метод Эйлера. Согласно этому методу интегрирование приближенно выполняется с помощью схемы, приведенной на рис. 3.3.

Рис. 3.3

Первый блок схемы осуществляет усиление входного сигнала в  Т  раз, второй – задержку на время  Т,  где  Т – шаг дискретного времени  t = nT,  n = 0, 1, 2, … . Схеме соответствует разностное уравнение:

y((n+1)T)–y(nT)=Тx(nT)     или     yn+1yn xn.

Применяя  z-преобразование,  находим дискретную передаточную функцию

     Q(z) = Y(z) / X(z) = T / (z–1).  

Ее можно было получить из непрерывной передаточной функции интегратора  Q(p)=1/p   с помощью замены переменных  

p = (z–1) / T.      (3.5)

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

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

Это стандартная формула замены при переходе от  Q(p) к Q(z)  по методу трапеций. При этом обе передаточные функции Q(p) и  Q(z)  будут дробно-рациональными передаточными функциями  одного порядка.

Операторный вывод этой формулы  основан на приближенном соотношении:

,  

где  e-pT – известная из теории преобразования Лапласа передаточная функция блока запаздывания на время  Т.

Выразим отсюда  р  через  z:  2+pT=2zpTz;       Tp(1+z)=2(z–1);   откуда

         (3.6)

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

 

На рис. 3.4, а в качестве иллюстрации показан результат интегрирования синусоидального сигнала с помощью трех цифровых интеграторов, использующих методы Эйлера (кусочно-постоянный сигнал), трапеций (кусочно-линейный сигнал) и Рунге-Кутта (непрерывная линия).

a)

б)

Рис. 3.4

Графики получены путем моделирования в пакете SIMULINK с помощью схемы, приведенной на рис. 3.4, б  (она соответствует значению  Т=1).

Формулы (3.5) и (3.6) позволяют осуществлять приближенный переход от непрерывной передаточной функции любого объекта  Q(p)  к его дискретной передаточной функции  Q(z).  Шаг квантования времени  Т  должен выбираться с учетом динамики объекта и быть по крайней мере на порядок меньше его постоянных времени или периода собственных колебаний.

Пример 7. Рассмотрим колебательное звено с передаточной функцией

Его импульсная весовая функция имеет вид синусоиды  q(t) = sin kt  с периодом  2p/k.  Найдем дискретную передаточную функцию этого колебательного звена.

Шаг  Т  выберем таким образом, чтобы на период собственных колебаний приходилось 20 точек дискретного времени:  T = 0,1p/k.  Используя, согласно формуле (3.5), замену  
p = (z–1)/T,  получаем

Применение формулы трапеций (3.6) дает более точный результат:

Учитывая, что  кТ = 0,314;    к2 Т2 = 0,0987,  получаем

 

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

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

Рассмотрим случай матричного описания дискретной системы

   (3.7)

где  – вектор состояния;  – входной и выходной сигналы,  A, B, C, D – матрицы соответствующих размеров.

Для того, чтобы найти дискретную передаточную функцию такой системы, применяем z-преобразование:

Выразим из первого уравнения  X(z)  и подставляя его во второе:

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

    (3.8)

Для систем с одним входом и одним выходом она представляет собой отношение двух полиномов

   (3.9)

причем знаменатель совпадает с характеристическим полиномом матрицы  А.

Формула (3.8) справедлива как для скалярных систем, так и для систем с несколькими входами и выходами. В общем случае  С  и  В – прямоугольные матрицы,  Q(z) – матричная передаточная функция, ее элементы – рациональные дроби типа (3.9).

Пример 8. Дискретная система задана уравнениями

Найти передаточную функцию от входа  u  до выхода  у  и импульсную весовую функцию.

В данном случае имеем

Выписываем матрицу  zE – А  и преобразуем ее

Дискретную передаточную функцию получаем по формуле (3.8):

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

Все эти выкладки можно сделать в тулбоксе Symbolic пакета MATLAB. Соответствующая программа имеет вид:

syms z;

А=[0 1;-6 -5];В=[0; 7];С=[0 1];  %описание в пространстве состояний

Q=С*inv(z*eye(2)-А))*В;   %передаточная функция

q=iztrans(Q);     %весовая функция

Q,q     

Q=7*z/z^2+5*z+6),q=7*(-2)^n-7*(-3)^n %результат

3.4. Задачи и упражнения

1. Найти функцию  f(k),  если ее z-преобразование имеет вид   

Решение. Выполняем разложение на простые дроби      

Переходя к оригиналам, получаем ответ:  

2. Решить  разностное уравнение    используя  z-преобразование.

3. Найти передаточную функцию и реакцию дискретной системы

yk+2 + yk+1 –  6yk = 3uk+1 –  3uk

с нулевыми начальными условиями на линейно возрастающее входное воздействие  uk = k  при  k0,  (uk = 0  при  k<0).

Решение. Применяем  z-преобразование:

где  Q(z) – дискретная передаточная функция.

Так как     то   

Далее следует выполнить разложение на простые дроби и перейти к оригиналам.

Ответ.    

4. Цепная схема представляет собой соединение  n идентичных Т-образных резистивных звеньев (рис. 3.5). Найти ее выходное напряжение  un,  если  u0 = 1В  и все резисторы одинаковы.

Указание. Эта задача отличается от примера 6 только краевым условием  un = un-1 /2,  откуда  C = –thn.

Ответ. Формула  для выходного напряжения имеет вид   

5. Найти дискретную передаточную функцию системы, заданной матрицами описания в пространстве состояний:

Ответ.  Дискретная передаточная функция имеет вид   


  1.  СТРУКТУРНЫЕ МОДЕЛИ РАЗНОСТНЫХ УРАВНЕНИЙ

4.1. Модели уравнений с начальными условиями

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

Альтернатива состоит в использовании структурного моделирования, когда разностному уравнению сопоставляется некоторая структурная схема на сумматорах, элементах задержки и других блоках. Возможности для такого моделирования имеются в пакетах VISSIM, SIMULINK и ряде других. Остановимся подробнее на двух принципах структурного моделирования разностных уравнений.

Первый из них хорошо известен. Согласно нему схема моделирования составляется так же, как это делается при моделировании дифференциальных уравнений методом понижения производной (методом Кельвина) [9]. Если, например, дано разностное уравнение второго порядка

,    (4.1)

то его сначала разрешают относительно "старшего" члена, в нашем случае это  x(t + 2):

.    (4.2)

Затем, предполагая,  что член  x(t + 2)  известен, пропускают соответствующий сигнал через два элемента задержки ЭЗ1  и ЭЗ2 (рис. 4.1), каждый из которых задерживает свой входной сигнал на один такт.  Выходные сигналы элементов задержки складывают в соответствии с уравнением (4.2) и получают сигнал  x(t + 2),  подавая который на вход первого элемента задержки, замыкают схему моделирования.

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

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

.   (4.3)

Вводя индексацию переменных, перепишем его в форме

.     (4.4)

Если нам известны начальные значения функции  x(t)

,

то, подставляя эти значения в формулу (4.4) при  k = 0,  можно найти  x(n).  Далее, полагая в формуле (4.4)  k = 1  и зная  x(1), ..., x(n),  можно найти  x(n+1);  полагая  k = 2,  можно найти  x(n + 2)  и т.д.

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

Рис. 4.2

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

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

1) Геометрическая прогрессия  (пример 2, разд.1) описывается разностным уравнением первого порядка    xk+1 = axk ,             x0 = 1.

Оно моделируется с помощью цепочки последовательно соединенных усилителей (рис. 4.3).  

Арифметическая прогрессия, уравнение электрической цепи и ряд Фибоначчи (примеры 3, 4, 5 разд. 1) описываются однородным разностным уравнением второго порядка вида

.     (4.5)

Схема моделирования такого уравнения при заданных начальных условиях  x0,  x1  приведена на рис. 4.4.

Рис. 4.4

4.2. Пример моделирования в пакете SIMULINK    

Проиллюстрируем методику моделирования методом шагов на примере разностного уравнения второго порядка

1,9x(t+2) – 0,81x(t+1)+ 0,9x(t)=0,   x0 =1,   x1 =–1.   (4.6)

Требуется получить его решение на интервале  0t12.

Переписываем исходное уравнение в виде

x(t+2)= 1,426x(t+1) – 0,474x(t).

 

Используя библиотеку блоков SIMULINK, создадим схему, показанную на рис. 4.5, и оформим ее в виде отдельного блока  с помощью  процедуры маскирования. Его входные сигналы  x(t)  и  x(t+1),  выходные сигналы  x(t+1)  и  x(t+2).

Рис.4.5

Чтобы получить решение в 12 точках, построим схему из 12 таких блоков, соединенных последовательно (рис. 4.6).  Блок "To Workspace" обеспечивает передачу результатов в рабочую область MATLAB для их последующей обработки.

На выходах схемы получаем искомые значения:    x =   [ 1    -1    -0.9     0.09     0.891    0.7209   -0.1531   -0.7866   -.5702          0.1948     0.6885    0.4443   -0.2197].

Рис. 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

Рис. 4.7

Аналогичным образом выполняется структурное моделирование нелинейных разностных уравнений и краевых задач.

4.3. Модели уравнений с краевыми условиями

Аналитическое решение краевых разностных уравнений практически ничем не отличается от решения уравнений с заданной начальной последовательностью. Однако с точки зрения компьютерного моделирования или численного решения эти задачи существенно различаются. В частности, для краевых задач не удается непосредственно применить описанный выше метод шагов (подумайте, почему?) и, следовательно, связанные с ним итерационный алгоритм решения и схему моделирования по методу Кельвина (рис. 4.1).

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

Проиллюстрируем это на примере.

Пусть требуется найти решение уравнения (4.5) с краевыми условиями

x(0) = x0,          x(k) = xk.

Схема моделирования этого уравнения показана на рис. 4.4. Для получения искомого решения достаточно на вход первого сумматора подать известную величину  x0  и плавно изменять величину  x1  до тех пор, пока напряжение на выходе схемы не станет равным заданной величине  xk . На этом моделирование задачи завершается.

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

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

В рассматриваемом случае решение дает метод, который можно назвать методом краевой структурной схемы: согласно этому методу структурная схема моделирования должна отражать не только само уравнение, но и вид краевых условий.

Структурная схема (рис. 4.2) была построена по уравнению (4.3), разрешенному относительно старшего слагаемого. Это позволило легко учитывать заданную начальную последовательность. Разрешая то же уравнение относительно других слагаемых, будем получать другие варианты структурной схемы, позволяющие учитывать различные виды краевых условий.

Например, уравнение (4.5) можно разрешить относительно среднего члена

и использовать это выражение для составления схемы моделирования.

Мы по-прежнему получим цепочку последовательно включенных сумматоров, однако теперь каждая точка решения формируется не на основе двух предыдущих, как это было на рис. 4.4, а на основе двух соседних точек – предыдущей и последующей (рис. 4.8). Благодаря

Рис. 4.8

 

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

Увеличение порядка моделируемого уравнения не вызывает принципиальных трудностей. Пусть требуется решить уравнение (4.3) с краевыми условиями, часть из которых задана в начале интервала решения  (n1),  а часть – в конце  (n2),  причем  n1 + n2 = n:

Разрешим уравнение (4.3) относительно члена  x(t + n1)  и по полученному выражению составим схему моделирования. Она будет отличаться от схемы, показанной на рис. 4.3, тем, что у каждого сумматора будет  n1 входов, соединенных с предшествующими сумматорами и n2  входов, соединенных с последующими сумматорами. У крайнего левого сумматора на первую группу входов будут поданы напряжения  ,  а у крайнего правого сумматора на вторую группу входов будут поданы напряжения  .  Благодаря наличию в схеме прямых и обратных связей, на выходах всех сумматоров мгновенно установится решение, удовлетворяющее краевым условиям.

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

4.4. Матричное решение разностных уравнений

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

    (4.7)

с краевыми условиями  х0 = а,  хn+1 = b.

Для того, чтобы найти решение этого уравнения на интервале    выпишем систему линейных алгебраических уравнений

   (4.8)

Перейдем к матричной форме записи:

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

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

Можно сказать, что структурные схемы, приведенные на рис. 4.1, 4.4, 4.8 отвечали различным численным алгоритмам решения системы (4.8). Например, схема рис. 4.4 реализует рекуррентный способ решения этой системы при известных значениях  х0, х1,  а схема рис. 4.8 – процедуру решения краевой задачи методом простой итерации. На матричном языке итерационный процесс описывается формулой

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

Пример. Решим матричным способом задачу Фибоначчи о кроликах. Она сводится к  разностному уравнению (1.8)

Его матричное представление имеет вид  АХ = F,  где    А – треугольная матрица 11 порядка с единицами на главной диагонали и с элементами –1  на двух нижележащих диагоналях. Два первые элемента вектора F равны  2 и 1, а остальные равны нулю.

Сформируем указанные матрицы в пакете MATLAB:

F=[2 1 0 0 0 0 0 0 0 0 0]’;  с =[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);

.

Решение получаем с помощью команды  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. Сформулируйте краевые условия в задаче об электрической цепи (пример 4 разд. 1.1).

5. В чем состоит метод краевой структурной схемы? Составьте с помощью этого метода схему моделирования электрической цепи (пример 4 разд. 1.1, пример 2 разд. 1.2, пример 6 разд. 3.2).

6. Составьте три варианта структурных схем для моделирования уравнения (4.5).

7. Найдите аналитическое решение разностного уравнения (4.6) и сравните его с графиком рис. 4.7, полученным путем моделирования в SIMULINK.

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

9. Выполните задание п.8  для электрической цепи, приведенной на рис. 3.5.

МРУ_5.doc   11.07.03

  1.  5. НЕЛИНЕЙНЫЕ РАЗНОСТНЫЕ УРАВНЕНИЯ

5.1. Источники нелинейных разностных уравнений

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

Отметим две формы записи нелинейных разностных уравнений:  в виде одного уравнения  порядка к

F(xn+k, ..., xn)=0;     (5.1)

и в виде системы уравнений первого порядка (описание в пространстве состояний)

Xn+1 = Ф(Xn),       (5.2)

Где   – вектор состояния, Ф – нелинейная векторная функция.

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

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

,     (5.3)

где  r – некоторый параметр.

Уравнение (5.3) можно рассматривать как явную формулу, которая позволяет по числу xn находить следующее число xn+1, и таким образом определяет всю последовательность {xn}. Можно сказать, что соотношение (5.3) соответствует некоторому итерационному процессу. Поэтому последовательность чисел {xn} часто называют итерациями начальной точки

Простейшие уравнения вида (5.3) были рассмотрены в разд. 1. При   соотношение (5.3) определяет геометрическую прогрессию со знаменателем  r.  При  – арифметическую прогрессию.

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

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

  1.  1. Моделирование процессов с дискретным временем.

Допустим, исследуется динамика численности популяции некоторых животных и существует возможность оценивать эту численность только раз в год. Тогда естественно временную переменную считать дискретной, которая может принимать только целые действительные значения п = 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),                  (5.4)

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

yn+1 =qyn(Nyn),                  (5.5)

Отображения (5.4) и (5.5) часто называют логистическим отображением.

Нас интересует вопрос о том, что произойдет с различными видами по прошествии достаточно долгого времени, то есть каковы аттракторы изучаемого отображения. Для ответа на него в этой простейшей модели достаточно выяснить, какой будет последовательность {xn},   при различных значениях  r.

2. Дискретизация дифференциальных уравнений.

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

Единицы измерения выбраны так, чтобы предельная численность популяции равнялась 1. Если
0 <
х(0) < 1,  то  х(t) ® 1 при  t ® ¥.  В этом нетрудно убедиться, проводя качественное исследование этой модели. Неподвижная точка  х = 0 является неустойчивой, точка  х = 1 устойчивой.

Применим к этому уравнению метод Эйлера, заменив производную приращением функции

где  h –шаг дискретного времени.

Обозначая    получаем

Мы вновь получили одномерное разностное уравнение вида (5.4). Анализ показывает, что его свойства похожи на свойства дифференциального уравнения только при небольших  h. Когда шаг по времени  h  превышает некоторое критическое значение, поведение этих двух объектов начинает качественно отличаться.

3. Решение нелинейных алгебраических уравнений.

Стандартным методом численного решения нелинейных алгебраических уравнений вида F(x) = 0 является метод простой итерации. В соответствии с ним для нахождения корней исходное алгебраическое уравнение переписывается в виде x = f(x),   и затем ему сопоставляется разностное уравнение   

xn+1 = f(xn),    n = 0, 1, 2, ... .

При определенных условиях решение этого уравнения  образует последовательности {xn}, сходящуюся к корню исходного уравнения. Соответствующий  итерационный процесс иллюстрируется  рис. 5.1,а.

Другой стандартный метод численного решения алгебраического уравнения вида  x =f(x) – метод Ньютона (метод касательных). В нем последовательность {xn} строится на основе формулы

,

т.е. для получения следующей точки проводится касательная к кривой в предыдущей точке.

Графическое изображение итерационного процесса по методу Ньютона показано на рис. 5.1,б.

Рис 5.1

В обоих случаях возникают разностные уравнения первого порядка вида (5.3).

4. Обработка экспериментальных данных.

Пусть некоторый сложный процесс, развивающийся во времени, характеризуется функцией  х(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).  

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

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

5.2. Решение нелинейных разностных уравнений

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

Остановимся на одном из таких случаев, когда  разностное уравнение имеет вид

    (5.6)

где  a,  b,  c,  d постоянные коэффициенты.

Это разностное уравнение вида (5.3) с дробно-линейной функцией  Его нелинейность связана с наличием операции деления в правой части.

Уравнение (5.6) можно переписать в эквивалентной форме

   (5.7)

Теперь нелинейность связана с наличием члена, содержащего произведение переменных.

Изложим два способа решения таких уравнений.

Способ 1 (замена переменной).

Выполняя в уравнении (5.6) замену переменных по формуле   получим

,       

Избавимся от знаменателей, выполнив переход к уравнению типа (5.7) и потребуем, чтобы коэффициент при нелинейном члене  zk zk+1 равнялся нулю:

Положим с = 1 (это всегда можно сделать, поделив числитель и знаменатель правой части формулы (5.6) на  с):

Если дискриминант этого квадратного уравнения положителен

то существует вещественное число ,  для которого указанная замена переменных приводит к линейному уравнению для  z.  

Найдя    и выполнив обратную замену переменных, получим решение исходного нелинейного разностного уравнения (5.6).

Пример 1. Требуется решить нелинейное разностное уравнение

Здесь  a = 2,  b = -1,  c = 1,  d = 0.  Выписываем квадратное уравнение для определения  :

Оно имеет единственный корень   = 1.  Делаем замену переменных    

Приводим к общему знаменателю

Раскрывая скобки, получаем линейное разностное уравнение     

Его общее решение имеет вид    откуда находим решение исходного нелинейного разностного уравнения

где  С – произвольная постоянная.

Способ 2 (переход к системе линейных уравнений).

Представим переменную    в виде отношения двух новых переменных    

Тогда уравнение (5.6) принимает вид

Приравнивая отдельно числители и знаменатели, получаем систему двух линейных разностных уравнений

или в матричной форме

Ее решение может быть записано с помощью одной из формул

где  X0 –вектор начальных условий,  zi  и  Hiсобственные числа и собственные векторы матрицы  А.

Окончательное решение получаем в виде отношения компонент вектора  Хk.

Пример 2. Найдем указанным способом решение предыдущего примера. Выписываем эквивалентную систему линейных разностных уравнений

Матричная запись этой системы:

Вычислим последовательные степени матрицы  А:

Следовательно, решение линейной системы имеет вид

Решение исходного уравнения получаем как отношение uk  и  vk:

где  – произвольная постоянная, зависящая от начального условия.

Видим, что решения, полученные двумя способами, совпадают.

Пример 3. Рассмотрим задачу об определении входного сопротивления цепной схемы, содержащей  n одинаковых Г-образных звеньев (рис. 5.2,а), если все сопротивления в схеме имеют величину 1 Ом.

Непосредственное вычисление входного сопротивления для случаев n =1, 2, 3, …  дает:

R1 = 2;  …

Мы видим, что в полученных дробях фигурируют отношения чисел Фибоначчи 1, 1, 2, 3, 5, 8, 13, 21, … . Можно предположить, что справедлива общая формула   где  – n-ое число Фибоначчи.  Тогда в пределе при n  ∞  получаем, что входное сопротивление схемы стремится к золотому сечению

Приведем аналитический вывод формулы для входного сопротивления. Пусть найдено входное сопротивление  Rn-1  схемы, содержащей  n–1  звено. Добавляя одно звено на входе, получим схему рис. 5.2,б. Находим ее входное сопротивление  Rn:

Мы получили разностное уравнение с дробно-линейной функцией в правой части.

Для его решения воспользуемся способом 2:

Отсюда получаем систему двух линейных разностных уравнений

Находим характеристическое уравнение этой системы, его корни и собственные векторы матрицы  А:

      

Отметим, что корни связаны с золотым сечением соотношениями

Общее решение системы имеет вид:

Коэффициенты  с1  и  с2  находим из начальных условий  и0 = 2,  v0 = 1:

Отсюда

Формулу для  Rn  получаем как отношение

Упростим эту формулу, используя уравнение золотого сечения    вытекающие из него тождества  

а также приведенные выше соотношения  

Последнее равенство следует из формулы Бине для чисел Фибоначчи, приведенной в разд. 1 (пример 10).

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

Пример 4. [12] Динамика эпидемии. В городе с 20 000 жителей появляются 50 инфекционных больных, что вызывает эпидемию. Предположим, что прирост больных за день пропорционален числу контактов больных и здоровых, т.е. произведению числа здоровых (еще не переболевших и не приобретших иммунитет) на число больных. Коэффициент пропорциональности (он характеризует скорость распространения эпидемии) примем равным 10-4.

Спрашивается: как развивается эпидемия – как изо дня в день меняется число больных? На какой день будет максимальное число заболевших?.

Решение.

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

В соответствии  с условиями задачи имеем два уравнения

Первое из них характеризует число заболевших в очередной  (n + 1)-й день, второе – число здоровых (еще не болевших) к  (n + 1)-му дню.

Система разностных уравнений нелинейная из-за наличия произведения  x(n)y(n), это не позволяет решить их аналитически.

Воспользуемся методом численного моделирования. Мы имеем задачу с начальными условиями  х(0) = 50,  у(0) = 20000, это позволяют шаг за шагом вычислять изменение числа больных по дням, полагая  n = 1, 2, 3 … .

Например, при  n = 1  получаем  

Программа для моделирования этой задачи   в пакете 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 больных.

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

х(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)

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

Требуется составить программу для вычисления членов этой последовательности при  ,  если известно, что  x1=1/e.  

Очевидный путь состоит в построении программы с простым циклом, вычисляющей члены  x2, x3 и т.д., непосредственно по формуле (5.8). Приведем ее вариант, написанный на языке пакета MATLAB:

x=zeros(1,50); x(1)=exp(-1);  

for i=2:50,                 

   x(i)=1-i*x(i-1);

end;

Вычисление первых четырех членов дает:

x1=1/e=0,3678794, x2=1–2 x1=0,2642411,  x3=1–3 x2=0,2072766,  x4=1–4 x3=0,1708934.

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

Попытка проверить качество вычислений путем подстановки вычисленных членов последовательности в формулу (5.8) только вводит в заблуждение: равенство будет выполнено с высокой степенью точности при всех  n, в том числе и  при  n=20.  Таким образом "очевидный" алгоритм оказывается практически неработоспособным. Это лишний раз говорит о необходимости критически относиться к результатам компьютерных вычислений. Более подробный разбор этого примера содержится в учебном пособии  [10].

5.3. Анализ стационарных точек и точек бифуркации

В тех случаях, когда нелинейные разностные уравнения не имеют аналитического решения, наряду с численными методами применяют методы качественного анализа. Один из важных вопросов, с которого обычно начинают исследование разностных уравнений (5.2) или (5.3) – это анализ стационарных точек. Так называются точки, удовлетворяющие уравнению  X = Ф(X)  в первом случае и   – во втором. Они характеризуют положения равновесия системы, когда ее следующее состояние совпадает с предыдущим. Отметим, что стационарные точки разностного уравнения совпадают с так называемыми неподвижными точками отображений  Ф(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.

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

Это разностное уравнение первого порядка. Для определения его стационарных точек (неподвижных точек функции  f) построим график в плоскости  (y, х),  где    (рис. 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 – полуось    Следовательно, описанный вариант метода простой итерации позволит найти только один из трех корней.

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

Критерий устойчивости стационарной точки. Пусть  хi – стационарная точка уравнения  xn+1 = f(xn).  Если    то эта точка – устойчивая (притягивающая), если же    то стационарная точка   xi – неустойчивая (отталкивающая). Случай, когда модуль производной функции  f(x)  в стационарной точке равен нулю, требует дополнительного исследования (анализа производных высших порядков).

Применим этот критерий к рассмотренному примеру 6. В этом примере функция f(x) и ее производная имеют вид:  

Вычисляя значение производной в стационарных точках  х = 0;  0,5;  2,  получаем:

В соответствии с критерием заключаем, что стационарная точка  х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  зависит от плодовитости популяции, условий питания и реальной площади.

Найдем стационарные точки этой модели, рассмотрев алгебраическое уравнение
x = rx(1–x). У него есть очевидный корень  x1 = 0.  Если численность популяции в начальный момент равна 0, то и она будет равняться 0 и далее. Второй корень равен  .  Таким образом, x1 и x2  – стационарные точки данной задачи. Проанализируем их устойчивость и найдем точки бифуркации (критические значения параметра r).  Это позволит ответить на вопрос, как будет себя вести популяция при разных значениях  r.

Рассмотрим отдельно случаи  0<r<1,   1<r<3,   3<r<4,  полагая везде, что  

Случай 1 (0<r<1). В этом случае имеется одна стационарная точка  х1 = 0 (корень х2 отрицателен). Чтобы проанализировать ее устойчивость, нарисуем кривую  ,  где  f(x)=rx(1–х)  

и прямую  y=x  (рис. 5.5).  Отложим  х1  по оси абсцисс, проведем вертикаль до пересечения с кривой    (точка  А),  затем из нее горизонталь до пересечения с осью  х.  Полученную точку пересечения обозначим через  х2.  Легко проверить, что    Взяв точку  х2  за начальную и повторив все те же операции, получим х3,  затем х4  и т.д.  Эта процедура называется построением лестницы Ламерея.  Она позволяет графически строить члены последовательности {xn}.

В данном случае эта последовательность стремится к нулю (спуск по лестнице Ламерея приводит в начало координат), таким образом, стационарная точка x1 = 0 – устойчивая. К этому же выводу приходим, вычисляя производную    Поскольку    то условие устойчивости выполняется.

Случай 2 (1<r<3). Как только коэффициент  r  становится больше единицы, стационарная точка  x1 = 0  теряет устойчивость, так как   Одновременно появляется вторая стационарная точка    Вычисляя для нее    получаем    График зависимости   от  r  приведен на рис. 5.6. Область устойчивости   на нем выделена штриховкой. Поскольку при  1<r<3  модуль производной меньше единицы, то на указанном интервале стационарная точка  х2 – устойчивая. Критическое значение  r1 = 1,  при котором появляется стационарная точка    и меняется тип устойчивости стационарная точка  x1,  представляет собой первую точку бифуркации.

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

По мере приближения ступенек лестницы Ламерея к точке    горизонтальные размеры ступенек стремятся к геометрической прогрессии

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

Случай 3 (3<r<4).  Как только  r  превышает значение 3, стационарная точка    теряет устойчивость, т.е. значения  r2 = 3 – это вторая точка бифуркации. В поведении системы происходят качественные изменения – в ней возникают устойчивые колебания. Их период сначала равен двум, затем, по мере увеличения параметра  r,  он удваивается и т.д. Значения  ri,  при которых происходят эти удвоения периода, являются точками бифуркации. По мере приближения  r  к критическому значению (оно примерно равно 3,5699) число бифуркаций удвоения периода стремится к бесконечности, а при достижении этого значения всякая регулярность в поведении системы пропадает и возникает так называемая хаотическая последовательность  {xn}.  Поведение системы в интервале  3,57<r<4  известно, как перемежающийся хаос, а при  r>4 – как сплошной хаос.

Описанная картина иллюстрируется графиками, приведенными на рис. 5.8–5.10. На них показано изменение численности популяции во времени при различных значениях коэффициента  r  и разных начальных условиях. Для удобства дискретные точки на графиках соединены прямыми линиями.

Рис. 5.8 относится к случаю  1<r<3.  График на рис. 5.8, а построен для  r =1,2;  х1 = 1/7,     Он получен в пакете MATLAB с помощью команд

x(1)=1/7; for i=1:25; x(i+1)=1.2*x(i)*(1-x(i)); end, plot(x)

Видно, что численность популяции монотонно растет, стремясь к устойчивой стационарной точке    На рис. 5.8, б  величина  r = 2,9  немного меньше бифуркционного значения    Последовательность  {xn}  по-прежнему стремится к равновесной точке   однако характер сходимости заметно изменился.

   

Рис. 5.9 относится к случаю  r>3.  При  r=3,1  наблюдается периодическое колебание численности между двумя постоянными уровнями (рис. 5.9, а).  При r=3,5  число уровней становится равным 4, происходит бифуркация удвоения периода (рис. 5.9, б). Наконец, график рис. 5.10, построенный для r=4, демонстрирует появление хаотических колебаний.

В рассмотренном примере мы столкнулись с простейшим нелинейным разностным уравнением вида   xn+1=f(xn)  и убедились, что даже при простых функциях  f(x)  оно демонстрирует сложное поведение.

Заметим, что хаотическая динамика может возникнуть и в непрерывных системах. Наиболее известный пример из теории дифференциальных уравнений – так называемая система Лоренца. Она представляет собой систему трех нелинейных дифференциальных уравнений, в которой при определенном сочетании коэффициентов возникает удивительный режим колебаний, получивший название «странного аттрактора». Траектории системы в трехмерном пространстве непредсказуемым образом блуждают вокруг двух особых точек.

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

5.4. Задачи и упражнения

1. Найти решение нелинейного уравнения   если известно значение y0.

Указание. Замена    приводит к линейному уравнению   

Ответ.     

2. Решить нелинейное разностное уравнение  

Ответ.      

3. Решить нелинейное разностное уравнение  

Ответ.      (гармоническая последовательность).

4. Решить нелинейное уравнение  

Ответ.     

5. Решить линейное нестационарное разностное уравнение

в котором коэффициент  а  зависит от номера такта.

Указание. Перейти к линейному стационарному уравнению второго порядка.

6. Решить нелинейное разностное уравнение  

Указание. Замена    приводит к линейному уравнению второго порядка.

7. Существуют три классических определения среднего значения двух чисел a и b:

среднее арифметическое, среднее геометрическое и среднее гармоническое

,   ,  .

В любой арифметической прогрессии каждый член представляет собой среднее арифметическое соседних членов, в геометрической прогрессии – среднее геометрическое соседних членов. Обе они описываются линейными разностными уравнениями. Последовательность, начинающаяся с чисел 1 и ½, называется гармонической, если каждый ее член представляет собой среднее гармоническое соседних членов. Постройте эту последовательность, получите для нее разностное уравнение и составьте разностное уравнение для ее частных сумм sk.

Ответ. Возможны несколько вариантов уравнений, приведем два из них:

 1)        2)  

8. Человек держит за конец резиновый жгут длиной 1м, привязанный к дереву. У другого конца жгута сидит жук. Каждую секунду жук проползает 1 см по жгуту. Каждую секунду человек, держа конец жгута, удаляется от дерева на 1м. Доползет ли жук до человека и если да, то за какое время?

Решение. Доля жгута в процентах, пройденная  жуком за n секунд, описывается разностным уравнением  S(n)= S(n–1)+1/n, откуда  S(n)=1+1/2+1/3+1/4+…+1/n.  Известно, что с ростом n эта сумма неограниченно возрастает. Для решения надо найти такое n, чтобы величина  S(n) достигла (или превысила) 100%. Заметим, что попытка прямого компьютерного определения искомого 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,  ne99,42 секунд.

Ответ. Теоретически жук доползет до человека, однако на это потребуется около 4.8*1032 тысячелетий. И место встречи будет находиться слишком далеко.

9. Исследовать нелинейное разностное уравнение  второго порядка    Построить график решения для 1< n <15 при  х1 = 1,  х2 = 2.

Ответ.  График решения при  х0=1,  х1=2  приведен на рисунке.  

10. Дано нелинейное разностное уравнение второго порядка  Найдите  х2007,  если  х1 = 1799,  х2 = 1828  (это годы рождения наших великих писателей).

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

Ответ.  Поскольку остаток от деления 2007 на 5 равен 2, получаем  х2002 = х2 = 1828.

clear x;x(1)=1;x(2)=2; for i=2:14; x(i+1)=(x(i)+1)/x(i-1); end, plot(x)

11. Рассчитайте первые 50 точек решения линейного неоднородного нестационарного  разностного уравнения первого порядка    Постройте график решения. Как оно ведет себя при больших  значениях  n?

Предупреждение. Будьте осторожны – попытка рекуррентного расчета точек  х2,  х3, …  на компьютере может ввести в заблуждение.

Ответ. Последовательность  {xn}  положительна и монотонно убывает, стремясь при больших значениях  n  к нулю по гиперболическому закону вида

12. Найти стационарные точки нелинейного разностного уравнения   

и выяснить их тип. Построить график решения для   при  х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 имеют место хаотические колебания.

Библиографический список

  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.

1 Мироновский Л.А. Диагностирование линейных систем методом комплементарного сигнала. //Приборы и системы. Управление, контроль, диагностика. 2002, N5, с.52-57.  


R

R

R

R

r

r

r

u2

u1

un

u0

r

un-1

-1

-.5

0

.5

1

-1

-.5

0

.5

1

Тi(x)

x

EMBED Word.Picture.8  

EMBED Word.Picture.8  

РРис. 1.5

G

H+G

H

Рис. 2.1

EMBED Equation.3  

3

EMBED Equation.3  

3

EMBED Equation.3  

x3

x2

x1

a

b

Рис. 2.2

c

Рис. 2.3

z

0

EMBED Equation.3  

y0

EMBED Equation.3  

x

x0

EMBED Equation.3  

z

z0

y

Рис. 2.4

EMBED Equation.3  

х3

с

EMBED Equation.3  

b

EMBED Equation.3  

x1

а

х2

EMBED Word.Document.6  

0

5

10

15

20

25

30

0

1

2

3

4

5

 

Рис.3.1

0

2

4

6

8

10

0

200

400

600

800

В

B

yn+1

yn

b

ЭЗ

k+1

y

x

Т

Э З

z

1

Unit Delay

Sine Wave

Scope4

р

1

Integrator

z

1

Unit Delay

0.5

Gain

Рис. 3.5

un

u0

u1

un-1

. . .

x(t)

u(t)

-a1

b

-a0

EMBED Equation.3

Э З2

Э З1

x(t+2)

x(t+1)

Рис. 4.1

EMBED Equation.3

x0

x1

xn-1

x2

x3

xk

xk+1

x1

x2

a1

...

an-1

EMBED Equation.3

EMBED Equation.3

a0

a0

a0

a1

...

a1

...

an-1

an-1

xn+2

xn+1

xn

EMBED Equation.3

a0

a1

...

an-1

xn+k

a

a

a

a

x0

x1

x2

xk

x3

Рис. 4.3

x1

x0

x2

a0

a0

a0

a0

a1

a1

a1

a1

xk

x4

x3

-

-

-

-

1.426

-0.474

Out2

2

2

1

Sum

1

In1

In2

Out1

Gain2

Gain1

Блок1

Блок 12

-1

X1

1

X0

x

To Workspace

Mux

X2

X12

xk

b

a

b

b

b

x2

x1

x0

xk--1

xk-2

x4

x3

a

a

a

a

b

u2

u1

un

u0

Rn-1

1

1

u0

Рис. 5.2,а

Рис. 5.2,б

EMBED Word.Picture.8  

Рис. 5.3

-1

0

1

2

3

4

5

 

.5

1

1.5

2

2.5

x

y

Рис.5.4

-1

1

2

1

2

3

r

EMBED Equation.3  

Рис. 5.6

 x

0

.2

.4

.6

.8

 1

0

.1

.2

.3

.4

 x

 y

A

 B

 C

r=1

Рис.5.5

0

.2

.4

.6

.8

 1

0

.1

.2

.3

.4

.5

 f

 x2

Рис.5.7

б)

a)

Рис. 5.8

Рис. 5.9

б)

a)

Рис. 5.9

Рис. 5.10


 

А также другие работы, которые могут Вас заинтересовать

22131. Осадка. Распределение накопленной деформации (εi) по объему осаженной заготовки 182 KB
  Расчет силы деформирования при осадке и построение графика технологических нагрузок. Мощность и работа пластической деформации при продольной осадке цилиндра. Работа деформирования при продольной осадке.Схема осадки:1 – нижняя плита; 2 – верхняя подвижная плита; 3 – цилиндрическая заготовка при продольной осадке; 4 – цилиндрическая заготовка при поперечной осадке.
22132. Метод баланса работ 36 KB
  В основу метода положено следующее положение: при пластической деформации работа внешних сил на соответствующих им перемещениях равна работе внутренних сил работе пластической деформации. Работа пластической деформации 2 Если упрочнение отсутствует то Чаще принимают равным выбранному по АВ – работа внешних сил: активной силы силы деформирования; сил трения. Работа сил трения берется со знаком минус. 3 где X Y Z – проекции силы действующей по участку поверхности dF на оси координат а UX UY UZ –...
22133. Феноменологическая теория разрушения металлов при холодной пластической деформации 98 KB
  Феноменологическая теория базируется на сложившихся в настоящее время физических представлениях о закономерностях разрушения металла при пластической деформации. Различными экспериментальными методами было показано что величина пластического разрыхления возрастает пропорционально степени деформации сдвига. Авторами данной теории была выдвинута следующая гипотеза: 1 где степень разрыхления частицы накопленная частицей деформация сдвига ab коэффициенты...
22134. Выдавливание. Расчет силы деформирования и построение графика технологических нагрузок 617.5 KB
  Основы теории штамповки выдавливанием на прессах М. Прямое выдавливание – технологическая операция в процессе которой происходит истечение металла 2 заключенного в замкнутой полости контейнер 3 в направлении движения рабочего инструмента 1 через отверстие поперечное сечение которого определяет поперечное сечение выдавливаемой части деформируемой заготовки. Обратное выдавливание – технологическая операция в процессе которой происходит истечение металла из замкнутой полости в направлении обратном встречном движению рабочего...
22135. Вытяжка без утонения 314 KB
  Схема операции – вытяжка из осесимметричной полой заготовки. При этом величина зазора между матрицей и пуансоном составляет не менее толщины исходной листовой заготовки Рис. Пример заготовки и детали.
22136. Вытяжка с утонением стенки 165 KB
  Механическая схема деформации и распределение деформации по очагу пластической деформации. Степень деформации при вытяжке оценивают коэффициентом вытяжки: или см. Частицы расположенные у нижней границы очага пластической деформации получают максимальную деформацию: . Частицы расположенные у верхней границы очага пластической деформации получают минимальную деформацию.
22137. Волочение 197 KB
  а б Рис. Рис. Допущения: напряжённое состояние плоское; продольные скорости металла одинаковы в пределах поперечного сечения ОПД очаг пластической деформации; и считаем главными напряжениями Сечениями z и zdz выделим элемент ОПД Рис. Рис.
22138. Метод верхней оценки 162.5 KB
  Сущность метода верхней оценки заключается в разбиении заготовки на жесткие блоки наделённые возможностью относительного скольжения и составлении баланса мощностей внешних и внутренних сил. При этом мощность пластической деформации рассчитывается как сумма мощностей сил трения по всем поверхностям скольжения жестких блоков относительно друг друга и инструмента. Скорости скольжения рассчитываются путём построения годографа скоростей. Строят годограф скоростей и определяют все скорости относительного скольжения всех блоков.
22139. Вырубка и пробивка 183 KB
  В верхнем небольшом по толщине слое металла примыкающем к пуансону. В нижнем небольшом по толщине слое металла прилегающем к матрице. 4 В срединном слое металла наибольшом по толщине двухосная схема напряжений и схема деформации сдвига. Местное поверхностное смятие развивается по толщине пока вся толщина металла не будет охвачена пластической деформацией; на третьей стадии происходит пластическая деформация в узкой по толщине кольцевой зоне пластический сдвиг.