Файл: Элементы математического моделирования в программных средах MATLAB 5 и Scilab (Андриевский Фрадков).pdf

= |sin(cjfcTo + 7гА:)| = |y[fc]| (для всех к = 0,1,2,...). Таким образом, полученные в результате дискретизации процессы у[к] и 2/i [к] оказались (с точностью до знака, или начальной фазы) неразличимыми, в то время как исходные непрерывные процессы 2/(f) и j/i(*)i конечно, различны.

Поскольку рассматриваемый процесс у(к) в общем случае не является периодическим с частотой П = 2ж/Т, вычисление его спектра с помощью рассматриваемой процедуры являет-

тельность исследуемой реализации должна быть достаточно большой для получения спектра с заданным шагом дискретности по частоте (Т 1/<5и>.)

Задача оценивания спектра по выборке имеет еще одну важную особенность: оценка спектральной плотности, получаемая по ЛПФ, как правило, является смещенной. Это значит, что даже неограниченное уменьшение шага дискретности не приводит к неограниченному убыванию погрешности оценивания. На самом деле средним значением оценки является не точное значение спектральной плотности на заданной частоте, а некоторое сглаженное значение, получаемое усреднением функции - сверткой с некоторым ядром ("окном") [37, 119]. Удачный выбор функции-"окна" позволяет уменьшить погрешность оценивания при заданном конечном

T=max(t) ;

Nt=2~12

T0=T/Nt;

1 Э т а п р о г р а м м а является продолжением программы, приведенной на с. 65.

70

В этом фрагменте программы с помощью функции max определяется длина реализации Т, задается число точек N для вычисления спектров и находится период дискретности Т0. Поскольку заданное число точек не совпадает с полученным в результате моделирования, а также в виду неравномерности моментов времени t} далее производится линейная интерполяция процедурой interpl . Для этого вводится массив моментов отсчета t i с постоянным шагом ТО и формируются подлежащие последующей обработке массивы y l i , y2i: 1

ti=0:T0:(Nt-1)*T0;

y l i = i n t e r p l ( t , y ( : , 3 ) , t i ) ; y 2 i = i n t e r p l ( t , y ( : , 1 ) , t i ) ;

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

Y l = f f t ( y l i , N ) ;

Y2=fft(y2i,N);

Sl=abs(Yl).~2/N;

S2=abs(Y2).~2/N;

Далее вычисляем частоту Найквиста и>п (wn) и шаг по частоте 6ш (dw):

wn=pi/TO;

dw=2*wn/N;

В рассматриваемом примере задано Т = 5 с, N = 4096, откуда следует ~ 2570 с"1, 6ш и 1.26 с"1. Частота Найквиста значительно превышает диапазон частот исследуемых сигналов. Построим графики спектров процессов на входе и выходе фильтра в диапазоне частот [0,100] с"1.

отах=100;

w=0:dw:отах;

lw=length(w); figure subplot(211)

1 Можно использовать и процедуру i n t e r p l q . Тогда соответствующий оператор следует записать с транспонированием массива t i , а именно как

71

stem(w,Sl(1:lw),'к.'),grid t i t l e ( 'S_l(\omega) ') subplot(212) stem(w,S2(l:lw),, k.'), grid

Результат выполнения программы показан на рис. 3.6.

Рис. 3.6. Спектры процессов на входе (Si) и выходе (52 ) ф и л ь т р а (3.10).

3.2.3.Исследование экологической системы

Выше приведена упрощенная математическая модель системы "хищник-жертва". Эта модель представляется уравнением Лотки-Вольтерра (3.7). Получим фазовый портрет системы (3.7) на плоскости (2/1,3/2) Для некоторых безразмерных значений параметров и переменных. 1

На рис 3.7 приведена 5-модель системы (3.7). Начальные условия интегрирующих блоков заданы переменными у1_0,

1 Аналитическое исследование задачи приведено в [52, 90].

72

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

Программа построения фазового портрета экологической системы

а=2; Ъ=1; с=1; d=5; у2.0=2;

Yl=[1.5:0.5:5];

f o r k=l:length(Yl) yl_0=Yl(k);

s i m ( ' s 2 . 2 ' )

p l o t ( y ( : , 1 ) , y ( : , 2 ) ) , h o l d on xlabel{'y _ 1'>,ylabel('y _ 2') axis( [0 12 0 i n f ] )

end hold off

В программе заданы фиксированное начальное значение 2/г(0) = 2 и массив Y1 начальных значений у\(0). Оператором for - end организован цикл по параметру к. Число итераций равно числу элементов массива Y1. После моделирования SIMULINK-программой s2_2 выводится фазовая кривая на плоскости (2/1,3/2)- Чтобы получить фазовый портрет, т.е. изобразить на одном графике множество фазовых кривых, после оператора plot выполняется оператор hold (фиксация) графика (устанавливается опция on). По окончании

73

цикла фиксация снимается (hold off). 1 Оператор axis служит для задания диапазона осей графика. Параметр inf задает автоматический выбор границы по соответствующему направлению. Результат моделирования представлен на рис. 3.8. Стрелочками помечено направление движения по фазовой траектории со временем. Довольно сложно (хотя и возможно) запрограммировать вывод этих стрелок, поэтому на данном графике они добавлены "вручную". Как видно из рисунка, в системе имеют место нелинейные колебания относительно некоторого состояния равновесия, амплитуда которых зависит от начальных условий. Часто эти кривые трактуются таким образом, что экологическая система имеет естественное состояние равновесия и уничтожение "хищников" (точка а на графике) приводит через некоторое время к резкому уменьшению числа "жертв" (точка 6). Если также учесть, что уменьшение численности популяции может привести к ее полному исчезновению, то следует сделать вывод о необходимости соблюдать осторожность при воздействии на экологические системы.

Рис. 3.8. Фазовый портрет системы "хищник - жертва" .

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

74

3.3.Модели состояния динамических систем

3.3.1.Модели общего вида

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

Eo(s) при s > t однозначно определяют (в силу (3.5)) значение Ex(s) при s = t. Для модели динамики популяций (3.7) состоянием является вектор х = со\(yi,y2 ).

В общем случае уравнения состояний - это системы дифференциальных или разностных уравнений первого порядка вместе с уравнениями для выходных величин. Начальное состояние представляет "память" системы о прошлом. Модель состояния непрерывной динамической системы записывается в виде [6, 95, 62]

(3.15)

Vi = Gi{xi,...,xn,ul,...,um,t)

У1 = Gi(xu...,xn,ui,...,um,t)

где иХ)...,ит - входные переменные; y\}...}yi - выходные переменные; хь ...,хп - переменные состояния. Вводя векторные обозначения, можно записать (3.15) в более компактном виде:

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

1 Здесь и далее вектор (столбец) с компонентами ii,i2, ... ,i n обозна-

чается X = со\(хх,хп).

75