1. Теоретичні відомості

Лінійне однорідне різницеве рівняння другого порядку з постійними коефіцієнтами має вигляд:

x(t + 2) + a₁ x(t +1) + a ₀ x(t) = f (t). (1)

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

^Xk+2 + ^а1 ^Xk+1 + ^a0 ^Xk = ^fk .

Різницеві рівняння третього й більш високих порядків записуються аналогічно. Відомі два основні методи розв'язку лінійних різницевих рівнянь - за допомогою характеристичного полінома й з використанням z-перетворення, аналогічного перетворенню Лапласа.

1.1. Розв'язок різницевих рівнянь за допомогою характеристичного полінома Будемо шукати розв'язок однорідного різницевого рівняння у вигляді x(t)=z\ де z - деяке число. Підставляючи x(t) у різницеве рівняння (1) при f(t)=0 і скорочуючи на z, одержуємо характеристичне рівняння

z² + a₁ z + a ₀ = 0.

Якщо його коріння z₁, z₂ речовинні й різні, то загальний розв'язок має вигляд:

^x(t) = С ^z1 + ^c 2^z 2. Якщо zz₂,, то в розв'язку з'являється лінійний множник

^x1^(t) = ^(c1 + ^C2^t) z1.

У випадку пари комплексно-сполученого корінь z₁,₂=oc±ifi розв'язок може бути записане в речовинній формі

x(t) = р^t(c₁ sin 9t + c₂ cos9t), р = yja² + —², ф = arctg — .

a

Тут р-р- модуль комплексного числа zi, а р-р- його аргумент.

Формули для рівнянь більш високих порядків виглядають так само, просто збільшується число доданків у розв'язку.

Приклад 1. Розв'яжемо різницеве рівняння x(t+2) -5x(t+1) + 6x(t) = 0. Характеристичне рівняння має вигляд z -5z + 6 = 0.

Його коріння речовинне й різні: z₁ = 3, z₂ = 2. Загальний розв'язок: х(ґ) = з₁ 31 + з₂ 2^Г.

Приклад 2. Розв'яжемо різницеве рівняння х(ґ+2) + 2х(ґ+1) + 4х(ґ) = 0.

Характеристичне рівняння: z + 2z + 4 = 0. У нього комплексне коріння: z₁₂ = -1 ±43 і. Їхнє положення на комплексній площині z = а +і'Р показане на мал. 1.

	І	, В
	2
	Zi
	/ \ \ / і \
	/ 1 \ / І \р ■ / 1 \ / 1 \ / 1 \ / 1 \ А>- і \ і \ і \	1 \ ~Хф=120° \
-2	V В	1 -1	2 а
	-2

Рис. 1. Модуль і аргумент корінь

2л

Модуль і аргумент корінь можна знайти безпосередньо на мал. 1: р = 2, р = .

/ \ „_t ✓ 2л . 2л .

Загальний розв'язок: x(t) = 2 (c cos-1 + c₂ sin 1).

Довільні постійні Cj знаходять, задаючи початкові умови. Нехай, наприклад, у прикладі 2 задані початкові умови x(0)=2; x(1)= -4. Записуємо загальний розв'язок для t=0 і t=1:

2 = 20c₁, 4 = 21(-c₁ -1/2 + c₂ -л/э/2).

2_л

Звідси знаходимо з]=2, з₂=0. Отже, розв'язок має вигляд x(t) = 2^t+1 cos-1.

Загальний розв'язок неоднорідного різницевого рівняння являє собою суму загального розв'язку відповідного однорідного рівняння й будь-який приватного розв'язку неоднорідного рівняння. Приватний розв'язок шукають у тому ж виді, що й права частина, тобто функція f(t) у рівнянні (1):

• якщо f(t) - постійна, то у вигляді константи;

• якщо f(t) - експонента, то у вигляді експоненти з тим же показником;

• еслиf(t) =sin kt або cos kt, те у вигляді Cisinkt+C₂Coskt.

Коефіцієнти з₁ і з₂ знаходять, підставляючи частки розв'язку в різницеве рівняння й дорівнюючи однойменні функції праворуч і ліворуч.

Приклад 3. Дане неоднорідне різницеве рівняння другого порядку

^xk+2 ^{- xk}+1 + ^xk = ⁶; ^Х0 = ¹; ^Х1 = ¹

Знаходимо коріння характеристичного полінома

₂ -1 ±3i ж 2 + г +¹ = ⁰; = 2 ; ^р = ¹; V = в.

Приватний розв'язок шукаємо у вигляді х_част=с. Підставляючи його у вихідне рівняння, знаходимо, що

^хчаст ².

Загальний розв'язок неоднорідного рівняння одержуємо як суму приватного розв'язку й загального розв'язку відповідного однорідного рівняння

• ж, ж,

х, = 2 + з sm— до + з, з — к.

^{до 1} 3 ² 3

Коефіцієнти з₁, з₂ знаходимо з рівнянь 1 = 2 + з₂, 1 = 2 + з, л/3/2 + з_2/2, звідки З = -В3, з ₂ = 1.

1.2. Розв'язок різницевих рівнянь за допомогою г-перетворення

При описі дискретних систем і розв'язку різницевих рівнянь широко застосовується апарат z-перетворення - це дискретний аналог перетворення Лапласа. Наприклад, множення зображення F(p) на оператор Лапласа р відповідає диференціюванню безперервної функції В(7). Множення зображення F(z) на оператор 2 відповідає зрушенню функції ДО (яка може бути безперервної, дискретної або ґратчастої) на один такт.

Таким чином, якщо оператори р і 1/р - це оператори диференціювання й інтегрування, то оператори г и г"¹ - це оператори зрушення вліво й вправо. З інженерної точки зору оператор 2"¹ являє собою елемент затримки.

Існують також певні паралелі між зображеннями функцій F(p) і F(z). Наприклад, зображенню за Лапласом F(p) = 1 відповідає дельта-функція ) = 5(1), а зображенню F(z) = 1, відповідає одиничний імпульс. У тому й іншому випадку оригіналом є елементарний імпульсний вплив.

Коротка таблиця z-перетворень інших функцій наведена в Додатку. Нехай дане різницеве рівняння п-го порядку

в(+п)+червоний.1в(+і-1)+ ... +аоу(0=/(0 (2)

с початковими умовами в(0)=в₀; в(1)=в₁; ...;в(п-1)=у_п-1.

Алгоритм його розв'язку за допомогою z-перетворення містить 4 кроку. Крок 1. Застосувати z-перетворення до рівняння (2), заміняючи ) на F(z), в(£) на 7); в(Г+1) на г(7(г)-уо) і т.д.

Крок 2. З отриманого алгебраїчного рівняння виразити 7(г). Крок 3. Виконати розкладання 7(г) на прості дроби. Крок 4. Користуючись таблицею, виконати зворотне z-перетворення. Результатом буде шуканий розв'язок різницевого рівняння.

Приклад 4. Потрібно розв'язати різницеве рівняння другого порядку

Уп+2 -5 Уп+1+ 6уп = ип

с нульовими початковими значеннями в₀, в₁ і вхідним сигналом и_п=1. Крок 1. Застосовуємо до нього z-перетворення

г² 7(г) - 5г7(г) + 67(г) = і(г).

z

z -1

Крок 2. Виражаємо Y(z) і підставляємо U(z) =

1

z

U(z)

Y ( z ) =

z² - 5z + 6 z² - 5z + 6 z - 1

Крок 3. Представляємо праву частину у вигляді суми простих дробів зі змінної г у чисельнику:

Y (z) =

г 0,5г

z 0,5z

- + -

(z - 1)(z - 2)(z - 3) z - 1 z - 2 z - 3

У звіті розкладання робиться методом невизначених коефіцієнтів. У пакеті MATLAB це можна зробити за допомогою команди residue.

(3)

Крок 4. За допомогою таблиці z-перетворень або команди iztrans тулбокса SYMBOLIC пакета MATLAB знаходимо оригінали кожного з доданків і складаємо їх: yn=0,5-2ⁿ +0,5-3 ".

1.3. Система лінійних різницевих рівнянь

Матричний запис системи лінійних однорідних різницевих рівнянь має вигляд

Х(Г+1)=АХ(Г), Х(0)=Хо,

де X е R" - вектор змінних, А - квадратна матриця постійних коефіцієнтів, Х₀ - вектор початкових умов.

Розв'язок цієї системи може бути записане в компактній статечній формі

Х(Г)=А Хо. (4)

При роботі в пакеті МАТЬАВ цей спосіб зручніше за все.

Приклад 5. Знайдемо розв'язок системи різницевих рівнянь другого порядку

х(г+1) = 2х(0, в(Г+1) = -2х( +2в(0,

якщо х(0) = 1; в(0) = 20.

" х(0)"		" 1"
_ В (0)_		20

Матричний опис завдання має вигляд

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

X(0) =

A=

" 2 0 - 2 2

Для одержання розв'язку скористаємося статечною формулою. Послідовно зводячи матрицю А в ступеня 2, 3, ..., п (проробіть це!), знаходимо

10 - n 1

10 - t 1

1

20 - t

X(t)=At X0= 2^l

An = 2

г г

Переходячи до скалярної форми запису, одержуємо остаточну відповідь х = 2, в = (20—) 2 .

Інший шлях розв'язку цього завдання пов'язаний з переходом від системи двох різницевих рівнянь першого порядку до одному рівнянню другого порядку відносно в:

в(1+2) - 4в(Г+1) + 4в(0 = 0.

Воно виходить, якщо виразити із другого рівняння х(), х(+1) і підставити в перше.

1.4. Моделювання дискретних систем вMATLAB і SIMULINK

В MATLAB при моделюванні лінійних системи з дискретним часом використовуються різницеві рівняння, моделі в просторі станів і дискретні передатні функції. У ядрі MATLAB є команда filter, що дозволяє розрахувати вихідний сигнал фільтра по відомому вхідному сигналу, заданому масивом своїх значень. Її вхідними аргументами є вектори а й b коефіцієнтів різницевого рівняння, а також масив значень сигналу й: y=filter(b,a,u). За допомогою четвертого вхідного аргументу можна задавати початкові умови фільтра.

Дискретна система в просторі станів, як і безперервна система, задається четвіркою матриць:

Xk₊₁ = Axk + Buk, yk = Cxk + Du _t. (5)

Тут Xk e Rn - вектор стану; uk, yk - вхідний і вихідний сигнали, A, B, C, D - постійні матриці.

Дискретною передатною функцією лінійної системи (5) називається відношення z- перетворення вихідного сигналу до z-перетворенню вхідного сигналу при нульових початкових умовах Q(z) = Y(z)/U(z). Для систем з одним входом і одним виходом вона являє собою відношення двох поліномів

_{Q( z)} ^bn-1 zn~' +... + ^bi z + ^b0

^{Q(z) =} ~n ni .

	1
X k +	0

-1 -1/2 1/2 -1

^X k+i -

z + an_₁ z +... + a₁ z + a₀

Дискретна передатна функція системи, заданої описом у просторі станів (5), може бути знайдена по формулі Q(z)=C(ze-A)^-1B+D.

Дискретні моделі в MATLAB задаються за допомогою тих же конструкторів ss, tf і zpk, що й безперервні з тою різницею, що останнім параметром задається частота дискретизації ts (sampling time).

Приклад 6. Розглянемо дискретну систему, описувану рівняннями

uk, yk =[0 1]X,

при кроці дискретизації, рівному 1.

Створимо відповідну ss-модель із дискретним часом:

>> a=[-1 -1/2 ;1/2 -1]; b=[1 ;0]; c=[0 1]; d=0; ts=1.0; s1=ss(a,b,c,d, ts).

На дисплей будуть виведені матриці:

a = -1 -0.5 b = 1 c = 0 1 d = 0 0.5 -1 0

і повідомлення: Sampling time: 1.0 Discrete-time model.

Аналогічно створюються tf і zpk-моделі з дискретним часом. Звертання до полів дискретної моделі здійснюється так само, як і в безперервної, наприклад для виводу матриці c треба набрати код s1 .c . Час дискретизації зберігається в поле ts.

Для моделювання дискретних систем використовуються ті ж команди (step, impulse, initial, lsim), що й у безперервному випадку. При одержанні вагової функції дискретних систем замість реакції на дельта-функцію розглядають реакцію на одиничний імпульс виду

. Для цього використовується та ж команда impulse. Реакцію на одиничний стрибок (перехідну функцію) одержують за допомогою команди step.

У тулбоксе SYMBOLIC є команди ztrans і iztrans для виконання прямого й зворотного перетворенн-z-перетворень.

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

Ук+i = ^xk, ₊i = ⁵ " ⁶Ук + ^7uk.

Знайдемо двома способами (символьним і чисельним) передатну функцію від входу u до виходу в.

, C = [0 1].

У цьому випадку маємо наступні матриці опису в просторі станів

	У		" 0 і"		"0"
X =		, A =		, B =
	x		- 6 5		7

Символьний спосіб. Дискретну передатну функцію одержуємо по формулі Q=C(ze-A)^-1B: >>syms z;

>>А=[0 1 ;-6 5];В=[0; 7];З=[0 1]; %уведення матриць

>>Q=C*inv(z*eye(2)-A))*B %передатна функція

Q=7*z/(za2-5*z+6) %результат

Чисельний спосіб. Використовуємо команди тулбокса CONTROL: >> A=[0 1;-6 5];B=[0; 7];C=[0 1];

>>sd=ss(A,B,C,0,1); % Discrete-time model Sampling time: 1, .

>> s=tf(sd)

7 z

Transfer function:

za2 - 5 z + 6

7 _z

В обох випадках одержуємо однаковий результат: Q( z) = — .

z - 5 z + 6

Моделювання дискретних систем в SIMULINK проводиться так само, як і безперервних. На малюнку наведені основні блоки з бібліотеки Discrete.

і		(z-1)		1		y(n)=Cx(n)+Du(n)		T		1
z+0.5		z(z-0.5)		1 +0.5z-1		x(n+1)=Ax(n)+Bu(n)		z-1		z

Discrete Discrete Discrete Filter Discrete State-Space Discrete-Time Unit De|ay

Transfer Fcn Zero-Pole Integrates

Блок затримки Unit delay дозволяє будувати схеми моделювання методом Кельвіна. Блоки Discrete Transfer Fcn (дискретна передатна функція) і Discrete Filter (дискретний фільтр) визначають раціональну функцію відносини поліномів від оператора z і зворотного оператора 1/z відповідно.

2. ЗАВДАННЯ ПО РОБОТІ Й ЗМІСТ ЗВІТУ

Робота складається із двох частин.

1. У першій частині досліджується дискретна система, задана різницевим рівнянням у_п+₂ + Ьу_п+₁ + ау_п = и_п з нульовими початковими условиямиу₀=0; уі=0.

Потрібно для свого варіанта знайти трьома способами її реакцію на вхідний сигнал

и_п=1:

• послідовно розрахувавши крапки y₂, ..., y₅;

• розв'язавши різницеве рівняння;

• використовуючи z-перетворення. Привести всі числові викладення за рішенням різницевого рівняння, дискретну передатну функцію, її розкладання на прості дроби, графік у_п, схему моделювання різницевого рівняння в SIMULINK, побудовану на елементах затримки 1/z і програми для обчислення в МЛТЬЛБ.

2. У другій частині роботи розглядається популяція (риб, тварин), розбитих на три вікові групи - молодшого, середнього й старшого віків. Задані величини р\, р₂ - імовірності дожиття особинами кожної вікової групи до наступного віку, і числа а\, а₂, а₃, що характеризують середню плідність кожної вікової групи.

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

(6)

Вікова структура популяції описується матричним рівнянням

	x		a1	a 2	a 3
X(k + 1) = AX(k), X =	y	, A =	P1	0	0
	z		0	p2	0

де х, в, z - чисельності трьох вікових груп.

Потрібно розрахувати чисельності вікових груп для п=1, 2, 3 при початкових умовах х₀= в₀= z₀=10. Теоретично знайти процентне співвідношення численностей х, в, z для більших п - воно визначається співвідношенням компонентів головного власного вектора матриці А.

Звіт повинен містити всі необхідні розрахунки, програми для обчислення в МЛТЬЛБ і схеми моделювання в SIMULINK. Для визначення власних векторів матриці А можна використовувати символьний варіант команди eig.

3. ПОРЯДОК ВИКОНАННЯ

• Одержати розв'язок різницевого рівняння й побудувати графіки

• за допомогою команд ztrans, iztrans, ezplot тулбокса SYMBOLIC;

• за допомогою команд tf і step тулбокса CONTROL;

• шляхом моделювання в SIMULINK;

• за допомогою команди rsolve пакета MAPLE.

• Розрахувати абсолютні й відносні чисельності вікових груп популяції

• за допомогою статечної формули (4);

• за допомогою команди initial.

Побудувати графіки у звичайному й логарифмічному масштабах.

4. КОНТРОЛЬНІ ПИТАННЯ

Дані різницеві рівняння:

1. x(t+2) -5x(t+1) + 6x(t) = 0, 2. x(t+2) -3x(t+1) + x(t) = 0,

3. x(t+2) - 3x (t+1) + 2x(t) = 0, 4. yk+1 + yk + 5y*-1 + 3yk-2 = 1,

5. yn+₂ = yn+₁ -yn, В1=1, В2=0, 6. xk+3 + 7xk+2 +15xk+1 + 9xk = 0,

7. xk+2-3xk+1 + 2xk = 0, 8. x(t+3)+7x(t+2)+15x(t+1)+9x(t) = 0,

9. х+4) + 2х(ї+2) + х(0 = 0,

^11. Ук₊1 ^- Ук + ² Ук-1 = ^0,

₁₃ Ук+2 + Ук = ^0, Уо = ^2, Уі = ¹

10. В„₊з-5В„₊2+8В„₊1-4В„=0; Уо=0, Уі=2, в2=10,

^12.13 Ук+1 + ⁴ Ук + Ук-І =

₁₄ Ук+1 ^{- 4} Ук + Ук-1 + ⁶ Ук-2 = ⁰ В0 = 6, В1 = 12, в₄ = 276.

Потрібно для рівняння, зазначеного викладачем:

а) Розв'язати його за допомогою характеристичного рівняння;

б) Розв'язати його за допомогою z-перетворення;

в) Намалювати схему моделювання й від неї перейти до матричної форми запису виду (3).

2. Розв'язати двома способами ( по формулі (4) і шляхом перехід до одному рівнянню) систему різницевих рівнянь

^{а) б)} , Г В„+1 - х_п = 0, В0 = 2,

Гх_га = + y_к^, Х0 = 1 Гх£ +¹⁾ = x^{(t) 3}в⁽0 ^в) |х_п+₁ -У_п = 0, х₀ = 0.

IУкн = х_к -3у_к, в0 = 1. Ьа+¹)=^зx(^t)+в(о. ^{п п}

17. Розв'язати систему Х_п+₁ = АХ_п для трьох варіантів матриці А:

	"-1	2	- 4"		"0	2	1"			"20
а) А0 =	- 8	- 3	2	, б) А =	1	0	2	, б) А2 = АТ;	X 0 =	0
	- 2	- 4	6		1	2	0			0