Лекции 2014 и 2010 годы (Лекц.Мод.сист 2010)

Посмотреть архив целиком

41



Московский Государственный Технический Университет им. Баумана

(МГТУ им. Баумана)

Факультет Информатика и системы управления


Кафедра Информационные системы и коммуникации (ИУ-3)













ЛЕКЦИИ по курсу “МОДЕЛИРОВАНИЕ СИСТЕМ ”

Кафедра ИУ-3, 4-й курс, 8-й семестр.




Автор: Боевкин Виктор Иванович





















Москва 2009



Оглавление

Введение

Раздел 1. Фильтрация помех при моделировании

1.1 Алгоритм разложения сигнала по неортогональному базису

1.2. О сходимости разложений

1.3. Разложение дискретных функций

1.4. Фильтрующие свойства разложений

1.5. Оценка уровня помех в обрабатываемой реализации

1.6. Оценка погрешностей разложения

1.7. Численная оптимизация по параметрам базиса

1.8 Оценка погрешностей разложения методом статистического моделирования



Раздел 2. Моделирование непрерывных систем


2.1. Цифровое моделирование непрерывных динамических систем

2.1.1 Методы цифрового моделирования

2.1.2. Численное решение линейных дифференциальных уравнений методом разложения в ряд Тейлора

2.1.3. Выражение ошибки численного решения через изменения коэффициентов дифференциального уравнения

2.1.4. Выражение ошибки численного решения через изменение корней характеристического уравнения

2.1.5. Устойчивость численного решения

2.1.6. Повышение точности численного решения методом коррекции уравнений движения

2.2. Погрешности аналогового моделирования


2.2.1. Особенности аналоговой вычислительной техники

2.2.2. Взаимосвязь приращений корней и коэффициентов характеристического уравнения

2.2.3. Влияние аналогового интегратора на корни характеристического уравнения

2.3. Погрешности полунатурного моделирования

Упражнения

Литература






































Введение


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

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

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

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

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

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

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

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

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

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

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


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

Слишком простое математическое описание может привести к ошибкам не только количественного, но и качественного характера, например – неправильной оценке факта устойчивости.

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

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

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





Раздел 1. Фильтрация помех при моделировании

1.1 Алгоритм разложения сигнала по неортогональному базису

При обработке экспериментальных функций времени может возникнуть необходимость представления этих функций в виде аналитических выражений. Непрерывный процесс f(t), определенный на интервале T, можно разложить по заранее выбранной системе базисных функций ф1(t)…фn(t), где n – число базисных функций. Разложение S(t) имеет вид:



S(t) = ; f(t) = S(t) + d(t) (1.1) Здесь d(t) = f(t) – S(t) (1.2)

-ошибка разложения

Введем понятие энергии Ey функции (сигнала) y(t):

Ey = (1.3)

В качестве критерия, характеризующего удаление S(t) от f(t) на интервале T, то-есть характеризующего качество разложения, примем энергию ошибки

Ed = (1.4)При заданных f(t) , T и фi(t) энергия ошибки является функцией от коэффициентов Ci разложения. Коэффициенты разложения будем выбирать из условия минимума энергии ошибки, для чего приравняем нулю частные производные

























= 0 ; i = 1…n . (1.5)

Из (4) и (2) можно получить:

Ed = = EfA + Es ; (1.6)

Ef = , A = -2 , Es =

С учетом (6) уравнения (5) станут:

= - + = 0; i = 1…n (1.7)

Используя (1) и (6) , выпишем частные производные:

= 0; = 2 (1.8)

= 2 ; = фi (t)

Подставив (8) в (7), получим:

= Vi; i = 1…n . (1.9)

Здесь введены обозначения:

Ui,k = ; Vi = . (1.10)



Линейная алгебраическая неоднородная система уравнений (9) вместе с обозначениями (10) определяет значения коэффициентов разложения Ci , при которых энергия Ed ошибки разложения минимальна. Используя обозначения (10), энергию ошибки (4) при произвольных значениях коэффициентов Ci можно представить:

Ed = Ef – 2 + . (1.11)

При оптимальных значениях коэффициентов, вычисленных из (9), энергия ошибки (11) имеет вид:

Ed = Ef - = EfEs . (1.12)

Приведенные выше результаты имеют достаточно прозрачное геометрическое представление. Введем в рассмотрение пространство H (бесконечномерное), элементами которого являются функции, определенные на интервале T, и удовлетворяющие условиям Дирихле, то – есть ограниченные и имеющие конечное число экстремумов и разрывов первого рода. Любая такая функция представляется в пространстве H точкой или, что эквивалентно, вектором. Определим скалярное произведение двух векторов в пространстве H следующим образом:

= (1.13)

Пространства, в которых определена операция скалярного произведения, называются Гильбертовыми. Они обладают геометрическими аналогиями, вследствие которых суждения, полученные из геометрического рассмотрения, имеют силу доказательств. Квадрат модуля любого сигнала y(t) в пространстве H равен энергии этого сигнала, что видно из определений энергии (3) и cкалярного произведения (13):

= = = Ey (1.14)

Векторы y(t) и x(t) ортогональны, если их скалярное произведение равно нулю:

= = 0 (1.15)

Базисные функции фi(t) изображаются в пространстве H векторами, которые формируют гиперплоскость Sn размерностью n, то есть набор функций фi (t) является неортогональным базисом гиперплоскости Sn. Разложение S из (1) всегда лежит в гиперплоскости Sn а коэффициенты Ci являются координатами S в неортогональном базисе гиперплоскости Sn. Разлагаемая функция f(t), в общем случае, не лежит в гиперплоскости Sn. Ошибка разложения d(t) изображается в пространстве H вектором, соединяющим концы векторов f(t) и S(t). Квадрат длины вектора ошибки, согласно (14), равен энергии ошибки, т.е. принятому критерию близости (4). Из геометрических соображений ясно, что длина вектора ошибки (а вместе с ней и энергия ошибки) будет минимальной, если вектор ошибки d(t) ортогонален гиперплоскости Sn. Последнее выполняется, когда вектор ошибки d(t) ортогонален ко всем базисным функциям фi(t), формирующим гиперплоскость Sn. Принимая во внимание определение скалярного произведения (13) и условие ортогональности (15), можно записать:

= = 0

Легко видеть, что подставив сюда d(t) = f(t) – S(t) из (2) и используя обозначения (10), получим систему уравнений (9). Величины (10) при этом являются скалярными произведениями:

Ui,k = = ; (1.16)

Vi = = .

Матричная запись системы уравнений (9) имеет вид:

U*C = V (1.17)

Здесь U – квадратная матрица размерностью n*n с элементами Ui,k , V – вектор правых частей размерностью n*1 с элементами Vi , C– вектор неизвестных размерностью n*1 с элементами Ck..


1.2. О сходимости разложений



При заданном базисе фi(t) коэффициенты Ui,k системы уравнений (9) могут быть вычислены заранее. В этом случае определение коэффициентов разложения Ci сводится к вычислению правых частей Vi по (16) и решению системы линейных уравнений (9) или, что эквивалентно, матричного уравнения (17). При большой размерности n базиса последняя операция достаточно трудоемка. Если выбрать такие базисные функции, которые в пространстве H ортогональны между собой попарно, то есть

= = 0; (1.18)

i = 1…n; j = 1…n ; i j ,

то матрица U становится диагональной, а система (9) – разрешенной. Разложение по ортогональному базису при возможности увеличения размерности базиса до бесконечности и выполнении некоторых других условий называется обобщенным рядом Фурье [1]. Простота решения системы уравнений (9) в большинстве случаев определяет использование ортогональных базисов (обобщенных рядов Фурье), в том числе и тригонометрических. Однако, в некоторых задачах предпочтение может быть отдано неортогональным разложениям. Степень гладкости разлагаемой функции можно характеризовать порядком R младшей производной, в которой на интервале разложения существуют разрывы первого рода (конечные разрывы). Известно [2], что сходимость тригонометрического ряда Фурье, под которой понимается уменьшение коэффициентов ряда с ростом его номера, имеет порядок 1/nR. Однако, при сколь угодно высокой степени гладкости разлагаемой функции внутри интервала разложения, несовпадение значений разлагаемой функции в начале и конце интервала разложения эквивалентно разрыву в функции, что приводит к весьма слабой сходимости ряда порядка 1/n. Это является следствием периодичности ортогонального тригонометрического базиса на интервале разложения. Возможность улучшения сходимости за счет отказа от ортогональности базиса проиллюстрируем простым примером. Разложим функцию f(t) = exp(-t) на интервале (0…2) по базису:



1, cos(pxt), sin(pxt), cos(2pxt),sin(2pxt), … cos(qpxt),.. (1.19)

Здесь q – число гармоник в базисе. Базис (19) ортогонален при x = 1. Отличие величины x от единицы характеризует удаление базиса от условий ортогональности. На рис.1 представлена зависимость от x среднеквадратической ошибки разложения

для базисов различной размерности: q = 1, q = 2, q = 3.

На Рис.2 представлены графики разлагаемой функции, разложения по ортогональным базисам различных размерностей (кривые q = 1, q =2, q = 3) и по неортогональному базису размерности q = 1 (кривая x = 0.1). Из графиков, представленных на рисунках 4.1 и 4.2, можно видеть, что увеличение размерности ортогонального базиса не приводит к существенному уменьшению ошибки разложения. Применение же неортогонального базиса позволяет существенно уменьшить ошибку разложения уже при малой размерности базиса. Увеличение объема расчетов на решение уравнений (17) может компенсироваться увеличением точности разложения.







1.3. Разложение дискретных функций

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

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

отсчетов:

f(tj) = f(j); j = 1…N.



Для разложения дискретных функций базисные функции и само разложение также должны быть дискретными:

фi(tj) = фi(j), S(j) = , i = 1…n, j = 1…N

Функциональное пространство H, введенное выше, становится N – мерным, а скалярное произведение определяется соотношением:

=

Энергия ошибки (4) (минимизируемый критерий) имеет вид:

Ed = =

Уравнения (9) и (17) для определения оптимальных коэффициентов разложения остаются в силе, но в параметрах Ui,k и Vi интегралы заменяются суммами:

Ui,k = ; Vi = (1.20)



1.4. Фильтрующие свойства разложений



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

f(j) = S0(j) + ah(j); S0(j) = ; j = 1…N (1.21)

Здесь S0(j) – идеальный (незашумленный) сигнал; фk(j), k = 1…n , - известные функции, составляющие идеальный сигнал; h(j) – случайный сигнал с известным законом распределения; а – коэффициент при случайном сигнале, с который условно назовем амплитудой помехи. В качестве базиса для разложения f(j) примем функции фk(j), входящие в идеальный сигнал (21). Целью разложения в этом случае является как можно более точное определение параметров C0k незашумленного сигнала, то – есть фильтрация помех. В результате разложения получим

S(j) = (1.22)

где коэффициенты разложения Ck вычислены по соотношениям (17) и (20).














На Рис. 3 представлено функциональное пространство H (для N = 3), плоскость Sn (для n = 2), а также все величины, входящие в (21) и (22). Вектор aSh(j) является проекцией вектора помехи ah(j) на плоскость Sn, то – есть разложением помехи по базису фk(j). Из треугольника а1а2а3 (рис.3) получим соотношение:

aSh(j) + d(j) = ah(j). (1.23)

Поскольку вектор ошибки d (отрезок a1a2)ортогонален плоскости Sn , то треугольник a1a2a3 является прямоугольным, поэтому, по аналогу теоремы Пифагора, из уравнения (23) можно получить:

a2ESh + Ed = a2Eh (1.24)

Здесь a2Eh– энергия помехи ah(j), a2ESh - энергия разложения aSh помехи ah(j), Ed – энергия ошибки разложения d(j) = f(j) – S(j). Фильтрующие свойства разложения можно характеризовать коэффициентом фильтрации Кф, отражающим уменьшение модуля помехи при проектировании ее на базисную гиперплоскость:

Кф = / = /. (1.25)

Из прямоугольного треугольника а1а2а3 (Рис.3) и (25) можно получить соотношение:

Кф = , (1.26)

где Eh – энергия случайного сигнала h(j), ESh – энергия разложения Sh(j) этого сигнала. При достаточно больших N эти величины обладают статистической устойчивостью и, согласно [4], их осредненные значения могут быть оценены теоретически. По определению, энергия случайного сигнала h(j)

Eh =

Осредним эту величину по ансамблю реализаций h(j):

Eh = M{} = = = NDh = N(Gh)2 (1.27)

Здесь Dh и Gh - дисперсия и среднеквадратическое отклонение случайной величины h(j) соответственно. При известном (или предполагаемом) законе распределения эти величины известны. Так, для равномерного закона распределения с единичным диапазоном

Dh = 1/12 = 0.083; Gh = = 0.289 (1.28)

Разложение случайного сигнала h(j) представим:

Sh(j) = (1.29)

Энергия этого разложения, согласно (12), имеет вид:

Esh = (1.30)

Величины Chk и Vhk связаны соотношениями (9) и (17). Введем матрицу Q, обратную U из (17), тогда

Q = U-1 = ; j = 1…n , i = 1…n . (1.31)

Решение уравнений (9) теперь можно записать:

Chk =

Подставив это в (30), получим:

Esh = * = (1.32)

Осредним (32) по ансамблю реализаций:

Esh = M[Esh] = (1.33)

В соответствии с (20) вычислим

M[Vhr Vhk] = M[*] =

= M[h(j)h(i)] (1.34)

Если случайные величины h(j) и h(i) независимы, то

M[h(j)h(i)]= (1.35)

Подставив (35) в (34) с учетом (20), получим:

M[Vhr Vhk] = Dh = Dh Urk (1.36)

Подставим (36) в (33):

Esh = Dh (1.37)

Внутренняя сумма в (37) является элементом произведения матриц Q и U, которые, согласно (31), взаимно обратны, поэтому

= Qkk Ukk = 1 (1.38)

Используя (38) , из (37) получим:

Esh = Dh n (1.39)

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

Кф = (1.40)



1.5. Оценка уровня помех в обрабатываемой реализации

При обработке экспериментальных реализаций представляет интерес уровень помех, присутствующих в этих данных. Для сигналов, которые можно представить в виде зашумленных сигналов (21), можно оценить амплитуду помехи (а). Из соотношений (24) – (26) можно выразить амплитуду помехи, присутствующую в эксперименте:

a2 = Ed/(EhEsh ) = Ed/Eh(1 – Кф2)

Энергия ошибки разложения Ed вычисляется в процессе разложения обрабатываемой реализации. Величины Eh, Esh, и Кф оцениваются выражениями (27), (39) и (40). Таким образом:

a2 = Ed/Dh(Nn) (1.41)

Напомним, что здесь Dh – дисперсия шума, N – число отсчетов в обрабатываемой реализации, n – количество базисных функций.

1.6. Оценка погрешностей разложения

Для зашумленных сигналов типа (21) в качестве погрешности dS(j) разложения естественно принять разность между разложением S(j) зашумленного сигнала и идеальным (незашумленным) сигналом S0(j).

Учитывая выражения (21) и (22) можно записать:

dS(j) = S(j) – S0(j) == (1.42)

Здесь dCk = Ck – C0k – погрешности разложения в области параметров сигнала.

Из рассмотрения треугольника a0a1a3 (рис.3) можно получить:

dS(j) = S(j) – S0(j) = аSh(j) (1.43)

Вектор аSh(j) был определен выше, как проекция помехи на базисную гиперплоскость Sn, амплитуда помехи а оценивается соотношением (41).

На Рис.4 представлен треугольник а0а1а3 из Рис.3, расположенный в базисной гиперплоскости Sn. При обработке экспериментальных данных вектор S(j) становится известным, вектора же S0(j) и аSh(j) – неизвестны. Однако длина L вектора aSh(j), характеризующего ошибку разложения, может быть оценена с помощью соотношений (39), (40),(41):

L = a= = Кф (1.44)




На Рис.4.4 приведена окружность (в общем случае гиперокружность) с центром в точке а3 и радиусом L (44). Истинное положение точки а3 (конца вектора S0(j)) лежит на

упомянутой окружности, что и является оценкой погрешности разложения.

Сравнивая соотношения (42) и (43) можно получить:

= a (1.45)

Отсюда


dCk = aChk , k = 1…n. (1.46)


Уравнения (46) показывает, что статистические характеристики ошибок коэффициентов разложения Ck могут быть оценены по статистическим характеристикам коэффициентов Chk разложения шума. Дисперсию Chk с учетом (31) можно записать:


D[Chk] = M[Chk 2] = (1.47)


Подставив сюда (36), получим:


D[Chk] = Dh (1.48)

Внутренняя сумма в (48) является элементом единичной матрицы в силу (31) , поэтому:


Отсюда и из (48)

D[Chk]=DhQkk (1.49)

Таким образом, дисперсию ошибок вычисления коэффициентов разложения сигналов типа (21) можно выразить из (46) и (49) следующим образом:

D[dCk]=a2DhQkk (1.50)1.7. Численная оптимизация по параметрам базиса



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

f(j) = S0 (j) + ah(j) , S0 (j) = , (1.51)

i = 1…m , j = 1…N

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

ф1(j,qi), ф2(j,qi) …фm(j,qi) , i = 1…m. (1.52)

Здесь qi – неизвестные параметры базисных функций, m – число этих параметров. Задав произвольную совокупность параметров qi, i = 1…m, произведем разложение сигнала (51) по базису (52), используя соотношения (17) и (20). В результате получим набор коэффициентов разложения и энергию ошибки, которые являются функциями принятой совокупности параметров qi:

Ck = Ck(q1…qm) , Ed = Ed(q1…qm) (1.53)

Произведя численную оптимизацию Ed(q1…qm) по параметрам qi, i = 1…m, получим:

Ck опт = Ck опт(q1опт…qmопт), Edопт = Edопт(q1опт…qmопт) (1.54)

Оптимальные параметры из (54)

qiопт , i = 1…m, Ck опт, k = 1…n, (1.55)

являются оценками параметров q0i , i = 1…m , и С0k , k = 1…n незашумленного сигнала (51).



1.8 Оценка погрешностей разложения методом статистического моделирования

Погрешности оценок (55) параметров q0i и C0k незашумленного сигнала S0(j) (50) можно произвести с помощью статистического моделирования. Задав значения параметров q0i и C0k и амплитуду шума а, сформируем ансамбль зашумленных реализаций, отличающихся одна от другой реализациями шума h(j). Разложение каждой реализации из ансамбля даст свою оценку (55) параметров q0i и C0k , что позволяет вычислить погрешности всех параметров для данной реализации:

dqi = q0i – qi опт, i = 1…m, dCk = C0k – Ckопт, k = 1..n (1.56)

Разложение ансамбля зашумленных реализаций позволяет сформировать статистики ошибок по всем параметрам, откуда можно получить оценки законов распределения ошибок и их моменты для заданной амплитуды помехи. При обработке единичной (уникальной) экспериментальной реализации амплитуда помехи определяется по (41), а статистическое моделирование, упомянутое выше, производится вокруг разложения S(j) этой экспериментальной реализации. Это позволяет получить статистические характеристики погрешностей параметров сигнала, представленного единичной реализацией.



Раздел 2. Моделирование непрерывных систем


2.1. Цифровое моделирование непрерывных динамических систем

2.1.1 Методы цифрового моделирования


При решении дифференциальных уравнений на ЦВМ каждая переменная представляется отсчетами, задаваемыми или вычисляемыми в дискретные моменты времени [2.1]. Значения отсчетов также имеют дискретную структуру, определяемую конечностью разрядной сетки машины.

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

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

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

tm+1 = tm + h, где h – шаг дискретизации, необходимо провести ряд математических операций над значениями переменных в предыдущие моменты времени, на что затрачивается машинное время dt. Для решения задачи в реальном масштабе времени необходимо, чтобы dt < h . Выполнение этого условия заставляет искать компромисс между скоростью и точностью решения с учетом, разумеется, допустимых для данной задачи погрешностей.

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

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

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

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

Глобальная ошибка - отличие точного решения (обычно неизвестного) от вычисленного на произвольном шаге интегрирования.

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

2.1.2. Численное решение линейных дифференциальных уравнений методом разложения в ряд Тейлора


Рассмотрим матричное однородное линейное дифференциальное уравнение

(2.1.1)


Одним из методов численного решения дифференциальных уравнений является разложение решения на шаге h в ряд Тейлора [2.4] относительно предыдущего момента времени t


(2.1.2)


Значения входящих в ряд Тейлора производных определим из уравнения (1):

, j = 1…r


Подставив эти значения в (2), получим:


X(t+h) = [I+Ah+A2h2/2…(Ah)r/r! …]X(t) (2.1.3)


Матричный ряд в квадратной скобке представляет собой матричную экспоненту:

eAh = [I+Ah+A2h2/2(Ah)r/r! …]


(2.1.4)


Таким образом, точное решение уравнения (1) на шаге можно представить в виде:

X(t+h) = eAhX(t) (2.1.5)


При численном решении ряды в соотношениях (2), (3) и (4) усекаются. Порядком r метода численного решения называется максимальная степень шага h , оставляемая в разложении. Представим матричный экспоненциальный ряд в виде суммы:


eAh = eAh + deAh , eAh = I+Ah+A2h2…(Ah)r/r! ,

(2.1.6)

deAh = (Ah)r+1/(r+1)! +…


где eAh - усеченная экспонента, deAh - отбрасываемая часть экспоненциального ряда (4).

Результат X(t+h) приближенного численного решения уравнения (1) на шаге h методом r -го порядка можно представить в следующем виде:


X(t+h) = [I+Ah+A2h2…(Ah)r/r!]X(t) = eAhX(t) (2.1.7)


Умножая матрицы (6) справа на X(t) с учетом (7), получим


X(t+h) = eAhX(t) = eAhX(t) + deAhX(t) = X(t+h) + dX(t+h) ,

(2.1.8)


где ошибка численного решения на шаге h или локальная ошибка метода

r -того порядка равна:


dX(t+h) = deAhX(t) = [(Ah)r+1/(r+1)! +…]X(t) (2.1.9)


Для приближенной оценки локальной ошибки можно использовать соотношение

dX(t+h)[ (Ah)r+1/(r+1)! ]X(t)


2.1.3. Выражение ошибки численного решения через изменения коэффициентов дифференциального уравнения


Сравнивая точное решение (5) с приближенным численным решением (7) можно видеть, что процедура вычисления на шаге h для них одинакова: предыдущее значение X(t) умножается слева на некоторую матрицу. Однако в точном решении X(t) умножается на матрицу полной экспоненты eAh , а в приближенном - на матрицу усеченной экспоненты еAh.

Зададимся вопросом: для какого, измененного по отношению к (1), уравнения приближенное численное решение X(t+h) из(7) является точным?

Измененное уравнение (1) запишем следующим образом:


(2.1.10)

Его точное решение на шаге h


X(t+h) = e(A+dA)hX(t) (2.1.11)


Приравняв выражения (2) и (11), можно получить :

e(A+dA)h = eAh (2.1.12)


Это уравнение определяет матрицу dA из уравнения (10).

Матричную экспоненту (12) представим в виде:


e(A+dA)h = eAhedAh


Отсюда и из (6) можно получить:


edAh – I = - e-Ah deAh (2.1.13)


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

dAh- (Ah)r+1 /(r+1)! (2.1.14)


Соотношения (13) и (14) позволяют вычислить матрицу dA, входящую в измененное уравнение (10). Точное решение этого уравнения совпадает с приближенным численным решением уравнения (1) методом r-того порядка. Таким образом, матрицу dA можно использовать в качестве характеристики погрешностей численного решения через изменение коэффициентов дифференциального уравнения (1).

Используя (14) можно выразить локальную ошибку численного решения (9) через dA:


dX(t+h) = - dAhX(t)


Если по условиям решаемой задачи можно сформулировать допустимые изменения коэффициентов dA, то соотношения (13) и (14) позволят осуществить предварительный выбор порядка метода r и шага дискретизации h.

2.1.4. Выражение ошибки численного решения через изменение корней характеристического уравнения


Приведем матрицу A к диагональному виду [2.6]:

A = S L S-1 (2.1.15)


Здесь L –диагональная матрица, составленная из собственных значений Zk матрицы А, являющихся корнями характеристического уравнения


det(IZA ) = 0 , (2.1.16)


S - диагонализирующая матрица, составленная из собственных векторов матрицы А.

Известно, что диагонализирующая матрица S сохраняет свой вид для любой степени А и для любого матричного полинома по степеням А, включая сходящиеся степенные ряды. Диагонализируя матричную экспоненту (4) и умножая ее слева на S-1 и справа на S, получим:


eLh = I + Lh + (Lh)2/2 + …+ (Lh)r/r! + (Lh)r+1 /(r+1)! +…

(2.1.17)


Здесь eLh – диагональная матрица со скалярными экспонентами

eZkh по диагонали.

Диагонализируя соотношения (6) и умножая их слева на S-1 и справа на S, получим:

eLh = eLh + deLh , eLh = I+Lh+L2h2…(Lh)r/r! ,

(2.1.18)

deLh = (Lh)r+1)/(r+1)! + …

Все матрицы в (18) диагональные.

Диагонализируем соотношение (12):


e(L+dL)h = eLh (2.1.19)

Напомним, что матрица dA в (12) характеризует ошибку численного решения в области коэффициентов решаемого уравнения (1).

Аналогично, диагональная матрица dL из (19) характеризует ошибку численного решения в области корней характеристического уравнения (16). Точное решение (19) относительно dL имеет вид:


dL = ln(eLh )/hL (2.1.20)

Диагонализируя (13), получим другю форму точного решения (19):


edLh – I = - e-Lh deLh (2.1.21)


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

dZk h- (Zkh)r+1 [1 – Zkh(r+1)/(r+2) +…]/(r+1) dZk h- (Zkh)r+1 /(r+1)! (2.1.22)


Как видно из выражений (20) - (22), смещение dZk зависит только от

k -того корня Zk и не зависит от других корней характеристического уравнения.

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

Z1 , ..., Zn эквивалентно точному решению другой системы уравнений со смещенными корнями Zk + dZk , k = 1,…n, где смещения корней определяется соотношениями (2.1.20)…(2.1.22).

Если по условиям решаемой задачи можно сформулировать

допустимые изменения корней dZk характеристического уравнения, то соотношения (2.1.20) …(2.1.22) позволят осуществить предварительный выбор порядка метода r и шага дискретизации h.

2.1.5. Устойчивость численного решения


Выше было показано, что численное решение системы дифференциальных уравнений (1) с матрицей L собственных значений эквивалентно точному решению системы уравнений со смещенной матрицей собственных значений L+dL, определяемой выражением (19). Так как входящие в него матрицы диагональны, то легко написать соответствующие скалярные соотношения:


eZkh = e(Zk+dZk)h = ebkh , bk = Zk+dZk , k=1…n (2.1.23)


где bk – смещенные собственные значения.

Устойчивость исходной системы уравнений (1) или, что то же самое, ее точного решения, определяется следующими условиями:


Re(Zk) < 0 , k = 1…n.


При численном решении дифференциального уравнения корни смещаются на dZk . Это может привести к тому, что какой-то смещенный корень bk переместится в правую полуплоскость, что будет соответствовать неустойчивости численного решения.

Условием устойчивости численного решения является, очевидно, соотношение:

Re(bkh) = Re[(Zk +dZk)h] < 0 , k=1…n. (2.1.24)


Из (23) с учетом (18) найдем значение смещенного корня:


bkh = ln(eZkh) = ln[1+Zkh+…+(Zkh)r/r!] (2.1.25)


Представим eZkh в виде комплексного числа:


eZkh = u + iv = Neiф , (2.1.26)

N ==, ф= arctg(v/u)


Учитывая (25) и (26) , условия устойчивости численного решения (24) можно записать:


Re(bk ) = ln < 0 (2.1.27)


В координатах u и v из (26) областью устойчивости является круг единичного радиуса с центром в начале координат:


u2 + v2 < 1 (2.1.28)


Соотношения (27) и (28) позволяют избежать грубых ошибок при выборе шага дискретизации h и порядка

метода r.


2.1.6. Повышение точности численного решения методом коррекции уравнений движения


Приближенное численное решение X(t+h) = eAh X(t) дифференциального уравнения (1) эквивалентно точному решению Х(t+h) = e(A+dA)h X(t) измененного уравнения


,


где dA определяется из соотношения


eAh = e(A+dA)h


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

X(t+h) = eAh X(t) (2.1.29)


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


(2.1.30)


Численное решение уравнения (30) на шаге можно записать так:


X(t+h) = e(A+dA)h X(t) (2.1.31)


Здесь

e(A+dA)h = [I + (A+dA)h + …(A+dA)rhr /r! ] (2.1.32)


- усеченная экспонента матрицы (A+dA) .

Приравнивая решения (29) и (31), можно получить соотношения:


eAh = e(A+dA)h ; e(A+dA)heAh = deAh (2.1.33)


Отсюда, используя соотношения (6) и (32)и пренебрегая старшими степенями dA , получим :


dAh [I + Ah + …+ (Ah)(r-1) /(r-1)! ] = deAh (2.1.34)


В первом приближении оценку dA получим из (5.34) с учетом (6):


dAh(Ah)(r+1) /(r+1)! ( 2.1.35)


Сравнивая (35) с (14) можно видеть, что в первом приближении


dA-dA


Уравнения (5.33)…(5.35) с разной степенью точности определяют корректирующую матрицу dA , которую нужно прибавить к матрице исходного уравнения (5.1), чтобы численное решение скорректированного уравнения (5.30) совпало с точным решением на шаге исходного уравнения (1). Диагонализируя соотношения (33)…(35), получим:


eLh = e(L+dL)h ; e(L+dL)heLh = deLh


dLh [I + Lh + …+ (Lh)(r-1) /(r-1)! ] = deLh (2.1.36)


dLh = (Lh)(r+1) /(r+1)!


Все выражения в (36) являются диагональными матрицами, поэтому могут быть переписаны в скалярном виде:


eZkh = e(Zk+dZk)h ; e(Zk+dZk)heZkh = deZkh


dZkh [I + Zkh + …+ (Zkh)(r-1) /(r-1)! ] = deZkh (2.1.37)


dZkh = (Zkh)(r+1) /(r+1)! k= 1…n


Уравнения (37) с разной степенью точности определяют

приращения dZk , на которые нужно скорректировать

собственные значений Zk матрицы А, чтобы численное решение измененного таким образом уравнения совпало с точным решением на шаге исходного уравнения (1).

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

- неоднородные уравнения [2.5]

- линеаризуемые уравнения

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



2.2. Погрешности аналогового моделирования


2.2.1. Особенности аналоговой вычислительной техники


В аналоговых вычислительных машинах (АВМ) каждая переменная величина из решаемой системы дифференциальных уравнений представляется (моделируется) напряжением на соотвтствующем блоке [2.7]. Каждый блок –это решающее устройство, выполняющее свою математическую операцию одновременно с другими блоками, то-есть параллельно во времени. Основой решающего блока является операционный усилитель постоянного тока с большим коэффициентом усиления и глубокой отрицательной обратной связью. Основными решающими блоками являются сумматор и интегратор. Отрицательная обратная связь сумматора осуществляется через активное сопротивление, интегратора – через конденсатор. Стабильность коэффициентов передачи решающих блоков, а вместе с тем и точность выполнения математических операций, зависит от качества используемых сопротивлений и конденсаторов, и не может быть высокой. Изучение влияния погрешностей выполнения операций суммирования и интегрирования на точность решения дифференциальных уравнений рассмотрим на примере уравнения n – ного порядка:


xn + an-1xn-1 + an-2xn-2 +…+ a1x^ + a0x = f(t) (2.2.1)



2.2.2. Взаимосвязь приращений корней и коэффициентов характеристического уравнения


Динамические свойства системы, описываемой уравнением (1), в большой степени определяются корнями характеристического уравнения:


P(p,ak) = pn + an-1pn-1 + an-2pn-2 +…+ a1p + a0 = 0 (2.2.2)


При решении уравнения (1) на аналоговой модели погрешности коэффициентов передачи решающих блоков формируют погрешности выставки коэффициентов уравнения (1), а вместе с тем и погрешности в коэффициентах характеристического уравнения (2). Это, в свою очередь, приводит к изменению корней характеристического уравнения, то – есть к изменению динамических свойств моделируемой системы. Для выяснения влияния погрешности выставки параметров аналоговой модели на изменение динамических свойств найдем зависимость малых приращений корней от малых приращений коэффициентов [2.8]. При известных корнях, характеристическое уравнение (2) можно представить ввиде произведения:

P(p,pj) = (p-p1)(p-p2)…(p-pn) = 0 (2.2.3)


Приравняв коэффициенты при одинаковых степенях p в (2) и (3), получим формулы Виета, выражающие коэффициенты характеристического уравнения через его корни :


an-1 = -(p1+p2…+pn)


an-2 = (p1p2+p1p3+…pn-1pn) (2.2.4)

……………………………

a0 = (-1)np1p2pn


Обратных прямых выражений корней через коэффициенты в общем случае не существует, однако для малых приращений они могут быть получены. Для этого возьмем полные дифференциалы по параметрам ak и pj от выражений (2) и (3) соответственно:



dP(p,ak) ===

= da0+pda1+…+pn-1 dan-1 (2.2.5)



dP(p,pj) == - (2.2.6)


Для обеспечения равенства полиномов (5) и (6) при любых p


dP(p,ak) = dP(p,pj) (2.2.7)


необходимо и достаточно, чтобы оно выполнялось при n различных значений p. В качестве этих значений удобно выбрать p = pk , k = 1…n. Подставив эти значения в выражение (7) с учетом (5) и (6) получим :


dpk= - (da0+da1 pk+…+pkn-1 dan-1 )/(pk-p1)

(pk-pk-1)(pk-pk+1)…(pk-pn), k =1,2,…n. (2.2.8)

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



2.2.3. Влияние аналогового интегратора на корни характеристического уравнения



Структурная схема аналогового интегратора представлена на рис.1.


Примерные значения параметров схемы [2.9]: коэффициент усиления операционного усилителя К = 108, входное сопротивление R = 1 Мом, емкость конденсатора С = 1мф, сопротивление утечки конденсатора Rk = 10000 Мом. Передаточную функцию Wa(p) схемы можно привести к следующему виду:


Wa(p) = Ka /(p + a) ; Ka = RC = 1сек ; a = 10-4 1/сек


С учетом инерционных свойств операционного усилителя, передаточную функцию Wa(p) аналогового интегратора примем:

Wa(p) = 1 /(p + a)(1+Tp) ; T = 10-4 сек (2.2.9)


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

Напряжения Uk , k = 0…n моделируют производные x(k) . При идеальных интеграторах с передаточной фикцией


Wи(p) = 1/p (2.2.10)


передаточная функция структурной схемы рис.5.2 имеет вид


W(p) = 1/( pn + an-1pn-1 + an-2pn-2 +…+ a1p + a0) , (2.2.11)


что соответствует уравнению (1).

Характеристическое уравнение имеет форму (2), а при известных корнях – форму (3). При наличии в уравнении (3) m вещественных корней и (n-m)/2 комплексных корней, его можно привести к виду:


P(p) = Пк=1(n-m)/2 (p2+2gkp+dk2i=1m(p+qi) = 0 (2.2.12)


Определение корней уравнения (12) сводится к приравниванию нулю каждого из сомножителей. Апериодическим составляющим переходного процесса соответствуют сомножители (моды) с вещественными корнями:


(p+qi) = 0 , pi = - qi , i = 1…m (2.2.13)


Колебательным составляющим переходного процесса соответствуют сомножители (моды) с комплексными корнями:


(p2+2gkp+dk2) = 0 , p1k,2k = -gk,

k = 1…(nm)/2 (2.2.14)


Рассмотрим теперь работу схемы на рис.2. при аналоговых интеграторах с передаточной функцией (9), которую представим в следующем виде:


Wa(p) = 1/D(p), D(p) = (p+a)(Tp+1) (2.2.15)


Здесь D(p) – оператор аналогового интегратора.

Сравнивая передаточные функции идеального (10) и аналогового (15) интеграторов, можно видеть, что передаточную функцию схемы моделирования рис.2.2. при аналоговых интеграторах можно получить, подставив в (11) вместо оператора p из (10) оператор D из (15):


W(p) = 1/( Dn + an-1Dn-1 + an-2Dn-2 +…+ a1D + a0)

(2.2.16)


Характеристическое уравнение (12), соотношения (13) и (14) для аналоговых интеграторов также формируются путем замены p на D. Для вещественных корней, из (13)


(D + qi) = (p + a )(Tp + 1) + qi = 0 (2.2.17)


Отсюда, при малых a и T, получим:


pi1 = - qi – a = pi – a ; pi2 = - 1/T + qi ; i = 1…m (2.2.18)


Из выражения (18) видно, что каждому вещественному корню рi характеристического уравнения (2) исходного уравнения соответствуют два корня схемы моделирования (рис.2.2) с аналоговыми интеграторами. Один из них – pi1 – смещенный на малую величину а корень pi . Второй, дополнительный, корень pi2 – располагается вблизи большой по модулю величины – 1/T .

Для комплексных корней (14), используя (15), получим уравнение четвертой степени:


(D2+2gkD+dk2) = p4 + 2p3T(1+aT) + p2 (1+ 4aT+a2T2+2gkT) + 2p(1+aT)(a+gk) + (dk2 + 2gka + a2) = 0 (2.2.19)


Проектируем приближенное (для малых a и T) решение уравнения (19) ввиде:


(p2 + 2gkp + dk2)(T2p2 + 2Tzkp +1) = 0 (2.2.20)


Приравняв коэффициенты при одинаковых степенях p в уравнениях (19) и (20), можно получить приближенные соотношения:


gk 2= gk2 + 2dka ; dk = dk +(a – gk2T)

zk = 1 – T(dk +a) (2.2.21)


Из выражений (20) и (21) можно видеть, что аналоговые интеграторы приводят к смещению основных комплексных корней (параметры gk и dk ) , и к появлению пары дополнительных комплексных корней с вещественной частью, близкой к –1/T и затуханием zk ,близким к единице.

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

Мерой пригодности аналоговых интеграторов для решения конкретной задачи является малость смещения основных корней и достаточно большая удаленность дополнительных корней от мнимой оси комплексной плоскости в сторону отрицательных значений. Для определенного типа аналоговой модели (с известными a и T) класс решаемых задач задается ограничениями на корни характеристического уравнения:

1/T>>[Re(pk)] >>a ; 1/T>>[Im(pk)] >>a




2.3. Погрешности полунатурного моделирования


При полунатурном моделировании система управления замыкается через преобразующее устройство ,как показано на рис .2.3, где обозначено: Wo(p) – передаточная функция разомкнутой системы управления (вместе с объектом управления), Wp(p) – передаточная функция преобразующего устройства [Л5.1].?

Пусть

Wo(p) = B(p)/A(p), Wp(p) = (1+G(p))/(1+Q(p))


B(p) =, A(p) =, an = 1 (2.3.1)


G(p) =, Q(p) =,


При идеальном преобразующем устройстве , когда Wp(p) = 1, передаточная функция замкнутой системы имеет вид:


WI(p) = B(p)/[A(p) + B(p)] ,


а характеристическое уравнение


RI(p) = A(p) +B(p) =0 . (2.3.2)


Корни характеристического уравнения (23) будем считать известными и обозначим


pI1 , pI2 , …pIn (2.3.3)


Характеристическое уравнение передаточной функции преобразующего устройства, и его корни, в соответствии с (22), запишем:


Rp(p) = 1 + Q(p) = 0 ,

pp1 , pp2…ppr (2.3.4)


Коэффициенты полинома Q(p) будем считать малыми, так, что модули корней (25) много больше модулей корней (24).В пределе, при Q(p) => 0,


limQ=>0 [pIj] =


Передаточная функция замкнутой системы (см.рис. 2.3.) с учетом преобразующего устройства и ее характеристическое уравнение имеют вид:


W(p) = B(p)(1 + Q(p))/[A(p)(1+Q(p))+B(p)(1+G(p))]

R(p) = [A(p)(1+Q(p))+B(p)(1+G(p))] = 0 (2.3.5)


Характеристическое уравнение (5) имеет n + r корней:


p1 , p2pn ,pn+1pn+r


Представляется очевидным, что при Q(p)0 и G(p)0 первые n корней уравнения (5) стремятся к корням уравнения (2), а последние r корней стремятся в бесконечность по отрицательной полуоси, как и корни уравнения (4).

Это позволяет предположить, что характеристический полином R(p) (5) при малых Q(p) и G(p) незначительно отличается от произведения:


R1(p) = RI(p)Rp(p) (2.3.6)


Введем мало измененные полиномы R, входящие в (6):


RI(p) = RI(p) + d RI(p), Rp(p) = Rp(p) + d Rp(p), (2.3.7)

R1(p) = RI(p)Rp(p), dRI(p) = dA(p) + dB(p), dRp(p) = dQ(p)


Здесь

dA(p) =, dB(p) =,

dQ(p) =. (2.3.8)


Величины dakdqk являются малыми приращениями соответствующих коэффициентов.

Приравняем теперь характеристический полином (5) измененному полиному (6):


R(p) = RI(p) , (2.3.9)


Подставив сюда выражения (7) для приращений полиномов, сокращая одинаковые слагаемые и пренебрегая произведениями малых приращений, получим:


dA +dB + QdB + AdQ +BdQ = - B(Q – G) (2.3.10)


Левая и правая части этого уравнения являются полиномами по p с наибольшей степенью pn+r . Для выполнения равенства (10) при любых p нужно приравнять коэффициенты при одинаковых степенях p в левой и правой частях уравнения, что даст n+r соотношений между приращениями dak , dbk , dqk . Количество этих приращений, в соответствии с выражениями (8), равно n+r+m+1. Следовательно, для выполнения равенства (10) можно задать произвольно m+1 приращение, что позволяет принять

dB = 0. Тогда соотношение (10) приобретет следующий вид:


(1 + Q)dA + (A + B)dQ = - B(Q – G) (2.3.11)


Система алгебраических уравнений , эквивалентная полиномиальному уравнению (11), - линейная и, в общем случае, неоднородная. Однако, при выполнении условий


Q(p) = G(p) , (2.3.12)


она станет однородной, т.е. будет иметь тривиальное нулевое решение dA(p) = 0, dQ(p) = 0. Этот результат становится очевидным, если вспомнить, что при выполнении условий (12) передаточная функция преобразующего устройства

Wp(p) = 1, что видно из (1).

Другой способ решения уравнения (11) исходит из того, что для обеспечения равенства двух полиномов степени (n+r) достаточно обеспечить их равенство в (n+r) точках по аргументу p. Для p , равных корням (3) уравнения(2), соотношение (11) приобретает вид:


dA(pIj)(1 + Q(pIj)) = -B(pIj)(Q(pIj) – G(pIj)); j = 1…n (2.3.13)


С учетом обозначений (7), систему уравнений (13) можно записать:


da0 + da1pIj +…dan-1pIjn-1 = (2.3.14)

= -{ B(pIj)(Q(pIj) – G(pIj))}/ (1 + Q(pIj)); j = 1…n


Решение системы уравнений (14) дает приращения коэффициентов полинома A(p) и , в силу соотношения dB = 0, приращения коэффициентов полинома RI(p) = 0 (2).

Зная приращения коэффициентов полинома (2), можно вычислить приращения его корней по соотношениям (2.2.8).

Сравнивая (2.2.8) с (14) можно видеть, что в них входят одинаковые комбинации приращений коэффициентов, поэтому не нужно решать (14) относительно dak , а можно сразу вычислять приращения основных корней системы по формуле


dpIj = {B(pIj)(Q(pIj) – G(pIj))}/ (1 + Q(pIj))(pIj - pI1)*


*(pIj - pI2)… (pIj - pIn) , j = 1…n (2.3.15)


Подставив в (11) корни (4), получим:


dQ(ppj)[A(ppj) + B(ppj)] = -B(ppj)(Q(ppj) – G(ppj));

j = 1…r (2.3.16)


Проведя аналогичные предыдущим рассуждения и выкладки, получим приращения остальных корней:


dppj = B(ppj)(Q(ppj) – G(ppj))/ [A(ppj) + B(ppj)]*


*(ppj – pp1) … (pIj - pIn); j = 1…r (2.3.17)

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

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


Упражнения


Колебательное звено без затухания описывается дифференциальным уравнением:

(3.1)


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

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

Решить уравнение (3.1) методами Рунге – Кутта первых четырех порядков. Определить фактическую частоту и декремент затухания (раскачки), сравнить с результатами п.1).

Произвести предварительную коррекцию уравнения (3.1), определить теоретические значения частоты и затухания (раскачки).

Решить скорректированное уравнение (3.1) методами Рунге – Кутта первых четырех порядков, сравнить с результатами п.3).

Вычислить смещение корней уравнения (3.1) при аналоговом моделировании с передаточной функцией интегратора (2.2.15). Найти границу по моделируемой частоте w , отделяющую затухание от раскачки.

Промоделировать уравнение (3.1) на аналоговой модели, сравить с результатами п.5).

Включить в схему моделирования уравнения (3.1) преобразующее устройство с передаточной функцией


W(p) = 1/(1+Tp)

Вычислить изменение частоты и декремент раскачки. Предъявить требования к допустимому значению T в зависимости от w.



Литература

Л5.1. Боевкин В.И. Оценка точности математического моделирования динамических систем. – М.: Изд. МГТУ, 1990. – 54 с.

Л5.2. Калиткин Н.Н. Численные методы. – М.: Наука, 1978. – 512 с.

Л5.3. Форсайт Дж., Малькольм М., Моулер К. Машинные методы математических вычислений. - М .: Мир, 1980. – 279 с.

Л5.4. Беки Дж., Карплюс У. Теория и применение гибридных вычислительных систем. – М.: Мир, 1970. – 483 с.

Л5.5 Боевкин В.И., Багдамян О. Н. Оценка точности решения задачи Коши по смещению корней характеристического уравнения. – Труды №474. – М.: МВТУ, 1987. – 82 с.

Л5.6 Стренг Г. Линейная алгебра и ее применения. – М. : Мир, 1980. - 454 с.


Л5.7. Анисимов В.В., Голубкин В. Н. Аналоговые вычислительные машины. – М.: Высш. школа, 1971. – 447 с.

Л5.8. Боевкин В.И., Айнутдинова И. Н., Кучминская А.И. Чувствительность динамических свойств линейной системы к изменению ее параметров. – Труды № 284 . – М.: МВТУ, 1979. – 170 с.

Л5.9. Боевкин В.И. Моделирование при испытаниях автоматических информационных устройств. – М.: МВТУ, 1985. – 72 с.


Случайные файлы

Файл
138118.rtf
27927.rtf
69626.rtf
33840.rtf
141247.rtf