Исследование RC-генератора синусоидальных колебаний

СОДЕРЖАНИЕ

1.ПОСТАНОВКАЗАДАЧИ3

2.МАТЕМАТИЧЕСКАЯФОРМУЛИРОВКАЗАДАЧИ4

2.1ДИФФЕРЕНЦИАЛЬНЫЕУРАВНЕНИЯЛИНЕЙНОЙ ЧАСТИ4

2.2Уравнениеусилителя4

2.3Конечно-элементнаямодель усилителя5

3.ФУНКЦИОНАЛЬНОЕПРОЕКТИРОВАНИЕПРОГРАММНОГОКОМПЛЕКСА6

4.МОДУЛИ РЕШЕНИЯСИСТЕМЫ ДИФФЕРЕНЦИАЛЬНЫХУРАВНЕНИЙ7

4.1 Описание методаРунге - Куттачетвертогопорядка7

4.2Описание алгоритмаодного шага8

4.3Блок - схемаалгоритмаодного шагапо методу Рунге- Кутта9

4.4Подпрограммаодного шагапо методуРунге-Кутта.10

4.5Описание алгоритмаметода Рунге- Кутта с автоматическимвыбором шага10

4.6Блок - схемаалгоритмаметода Рунге- Кутта с автоматическимвыбором шага12

4.7Подпрограммаметода Рунге- Кутта с автоматическимвыбором шага13

4.8Тестовая задача15

4.9Результатытестирования16

4.10Квадратичнаяконечно-элементнаямодель усилителя17

4.10.1 Описаниеалгоритма17

4.10.2Блок - схемаалгоритмамодели усилителя18

4.10.3Подпрограмма- модель усилителя18

4.10.4Решение тестовойзадачи19

4.11Подпрограммавычисленияправых частейсистемы уравнений20

4.12Подпрограммавывода20

4.13Главный модульрешения системыуравнений21

5.ЧИСЛЕННОЕМОДЕЛИРОВАНИЕАВТОГЕНЕРАТОРА22

5.1Пробные решения22

5.2Решение дляспектральногоанализа выходногонапряжения24

5.3Решения дляустановлениязависимостейпараметровот 25

6.ПРОГРАММЫОБРАБОТКИРЕЗУЛЬТАТОВМОДЕЛИРОВАНИЯ26

6.1Программачисленногоинтегрированияпо методутрапеций26

6.2Блок - схемаалгоритмавычисленияамплитуд гармоник27

6.3Результатыгармоническогоанализа28

7.ЛИТЕРАТУРА30

1.ПОСТАНОВКАЗАДАЧИ

ВыполнитьисследованиеRC-генераторасинусоидальныхколебаний(Рисунок.1)

Рисунок1

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

Передаточнаяфункция линейнойчасти

где

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

усилителяот его входногонапряжения

приведена втаблице 1

Таблица 1

U₁	-0,125	-0,1	-0,075	-0,05	-0,025	0	0,025	0,05	0,075	0,1	0,125
U₂	3	2,75	2,4	1,73	1	0,02	-1	-1,73	-2,4	-2,75	-3

Численнымиэкспериментамина ЭВМ найтизависимости:

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

,
амплитудыU_2maxвыходногонапряженияU₂(t)от

,
амплитудыA_n n-ойгармоники

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

.

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

Иззависимости

найти значение

,необходимоедля полученияпериода автоколебаний

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

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

2.МАТЕМАТИЧЕСКАЯФОРМУЛИРОВКАЗАДАЧИ

2.1ДИФФЕРЕНЦИАЛЬНЫЕУРАВНЕНИЯЛИНЕЙНОЙ ЧАСТИ

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

(1 )

(2 )

Введемпервую вспомогательнуюпеременную

,определяемуюиз уравнения

(3 )

Подставляя( 3 ) в ( 2 ), получаем

(4 )

Сокращаяна

и группируяв правой частичлены, не содержащие

,получаем

(5 )

Введемвторую вспомогательнуюпеременную

,определяемуюиз уравнения

(6 )

Подставляя( 6 ) в ( 5 ), получаем

(7 )

Сновасокращая на

и группируяв правой частичлены, не содержащие

,получаем

(8 )

Введемтретью вспомогательнуюпеременную

,определяемуюиз уравнения

(9 )

Подставляя( 9 ) в ( 8 ) и сокращаяна

,получаем

(10 )

Переходяв уравнениях( 10 ), ( 9 ), ( 6 ), ( 3 ) от изображенийпеременныхк их оригиналам,получаем системууравнений

(11 )

(12 )

(13 )

(14 )

Здесь

- функция, определяемаянелинейнойхарактеристикойусилителя.

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

2.2Уравнениеусилителя

Уравнение( 11 ) представляетсобой нелинейноеуравнение,которое необходиморешать прикаждом вычисленииправых частейсистемы.

Можнорешать этоуравнениеметодом итераций.Но есть болеепростой путь.

Найдемиз характеристикиусилителяразности

,а затем построимхарактеристику

Значение

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

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

Вычисленнаяхарактеристикапредставленав таблице 2.

Таблица 2

z₃	-3,125	-2,85	-2,475	-1,78	-1,025	-0,02	1,025	1,78	2,475	2,85	3,125
U₁	-0,125	-0,1	-0,075	-0,05	-0,025	0	0,025	0,05	0,075	0,1	0,125

2.3Конечно-элементнаямодель усилителя

Дляпостроенияквадратичногоконечногоэлемента используеминтерполяционнуюформулу Лагранжа

(15 )

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

находить

(16 )

Данныев этом случаенеобходимовыбирать изтаблицы 3, полученнойиз таблиц 1 и2.

Таблица 3

z₃	-3,125	-2,85	-2,475	-1,78	-1,025	-0,02	1,025	1,78	2,475	2,85	3,125
U₂	3	2,75	2,4	1,73	1	0,02	-1	-1,73	-2,4	-2,75	-3

3.ФУНКЦИОНАЛЬНОЕПРОЕКТИРОВАНИЕПРОГРАММНОГОКОМПЛЕКСА

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

программы- модели RC- генератора;
наборапрограмм обработкирезультатовмоделированияавтогенератора.

МодельRC - генераторадолжна, в своюочередь, включать:

модуль,вызывающийподпрограммуметода Рунге- Кутта;
модулиметода Рунге- Кутта;
модуль- модель усилителя;
модульправых частей;
модульвывода результатоводного шагаинтегрирования.

Дляпрограммнойреализацииметода Рунге- Кутта удобноиспользоватьдва модуля:

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

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

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

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

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

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

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

4.МОДУЛИРЕШЕНИЯ СИСТЕМЫДИФФЕРЕНЦИАЛЬНЫХУРАВНЕНИЙ

4.1 Описаниеметода Рунге- Кутта четвертогопорядка

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

Пустьимеется уравнение

или

где

-неизвестнаяфункция отнезависимойпеременной
;
-известнаяфункция.

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

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

вычисленныхв специальноподобранныхчетырёх точках.

Вкачестве первойточки берутначальную точкушага

Производнаяв этой точкеравна

где

-правая частьуравнения .

Вкачестве второйточки на плоскостирешения

выбираютточкус координатами

Производнаяво второйточкеравна

Длятретьей точкиберут координаты

ивычисляютпроизводную

Наконец,для четвёртойточки беруткоординаты

ивычисляютпроизводную

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

Теперь,находят координатыконечной точкишага.

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

4.2Описаниеалгоритмаодного шага

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

Таблица4

Имя	Тип	Назначение
PRAV	Подпрограмма.	Подпрограмма,возвращающаязначенияпроизводных.
N	Целый.	Порядокрешаемой системы.
XN	Вещественный.	Исходныймассив начальнойточки шага.
XK	Вещественный.	Результирующиймассив конечнойточки шага.
F	Вещественный.	МассиввозвращаемыхподпрограммойРRAVпроизводных.
TN	Вещественный.	Начальноена шаге значениенезависимойпеременной.
K	Целый.	Номерпеременной.
J	Целый.	Номерчастичногоприращения.
T	Вещественный.	Независимаяпеременная.
H	Вещественный.	Задаваемаявеличина шага.
P	Вещественный	Массивразмера (4,2),содержащийнеобходимыедля вычисленияи накопленияприращенийконстанты(0,.5,.5,1,6,3,3,6)
R	Вещественный	Рабочиймассив размера(N,3)

Блок-схемаалгоритмаизображенана Рисунке 2.Номер переменнойзаписан какверхний индекс.

Вцикле с номерамиблоков 2, 3, 4, 5 обнуляютсявторой и третийстолбцы рабочегомассива R.

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

Блок7 задает в циклепоследовательнозначения независимойпеременной: TN, TN+0.5H, TN+0.5H, TN+H.Константы0, 0.5, 0.5и 1содержатсяв первом столбцемассива Р.

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

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

Вцикле 13 - 16 во второйстолбец рабочегомассива повычисленнымзначениямпроизводныхзаписываютсяприращения,а в третьемстолбце накапливаютсясуммы четырехприращенийс весовымикоэффициентами1/6, 1/3, 1/3, 1/6 . Константы 6, 3, 3, 6для этого хранятсяво втором столбцемассива Р.

Вцикле 19 - 22 полученныеприращенияприбавляютсяк начальнойточке и результатзаписываетсяв выходноймассив.

Вблоке 23 вычисляютсяпроизводныев конечнойточке шага.

4.3Блок- схема алгоритмаодного шагапо методу Рунге- Кутта

Рисунок2

4.4Подпрограммаодного шагапо методуРунге-Кутта.

SUBROUTINESH(TN,H,XN,XK,F,PRAV,N,R)

DIMENSIONXN(N),XK(N),F(N),P(4,2),R(N,3)

DATAP/0.,.5,.5,1.,6.,3.,3.,6./

DO1 K=1,N

R(K,2)=0.

1R(K,3)=0.

DO4 J=1,4 ! Началовнешнего циклаопределения4-х приращений

T=TN+P(J,1)*H ! Заданиенезависимойпеременной.

DO2 K=1,N ! Циклзадания массивазначений зависимыхпеременных.

2R(K,1)=XN(K)+P(J,1)*R(K,2)

CALLPRAV(T,R,F,N) ! Вычислениепроизводныхв заданнойточке.

DO3 K=1,N ! Циклвычисленияи накоплениячастичныхприращений.

R(K,2)=H*F(K)

3R(K,3)=R(K,3)+R(K,2)/P(J,2)

4CONTINUE

DO5 K=1,N

5XK(K)=XN(K)+R(K,3) ! Вычислениепеременныхв конечнойточке.

CALLPRAV(TN+H,XK,F,N) ! Вычислениепроизводныхв конечнойточке.

RETURN

END

4.5Описаниеалгоритмаметода Рунге- Кутта с автоматическимвыбором шага

Блок-схема алгоритмаприведена наРисунке 3.

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

Таблица5

Имя	Тип	Назначение
PRAV	Подпрограмма.	Подпрограмма,возвращающаязначенияпроизводных.
OUT	Подпрограмма.	Подпрограмма,составляемаяПользователемдля выводарезультатов.
N	Целый.	Порядокрешаемой системы.
X	Вещественный.	Массивразмера (N,4).
R	Вещественный.	Рабочиймассив размера(N,3).
F	Вещественный.	Массивразмера (N,4).
TN	Вещественный.	Начальноена шаге значениенезависимойпеременной.
TK	Вещественный.	Конецинтервалаинтегрирования.
T	Вещественный.	Независимаяпеременная.
HМ	Вещественный.	Задаваемаявеличинамаксимальногошага.
E	Вещественный.	Задаваемоезначениеабсолютнойпогрешности.
EH	Вещественный.	Погрешность,вычисленнаяна шаге.
IER	Целый.	Выходнойкод ошибки.
H	Вещественный.	Текущийшаг.
HB	Вещественный.	Удвоенныйшаг.
T	Вещественный.	Текущеезначениенезависимойпеременной.
T1	Вещественный.	T1=T+H
T2	Вещественный.	T2=T+2H
KP	Целый.	Признакдостиженияконца интервалаинтегрирования.
KLP	Целый.	Признаквывода последовательнодвух вычисленныхточек.
K	Целый.	Индекс.

Второйстолбец массиваХдолжен содержатьвесовые коэффициентыпогрешности,на которыеумножаютсянайденныепогрешностикаждой интегральнойпеременной,чтобы послесложения этихпроизведенийполучить общийкритерий погрешностисистемы и сравнитьего с заданнойпогрешностью.Во втором столбцеони задаютсяс целью удобстваввода. Первыйстолбец массиваХзаполняетсяначальнымиусловиями, азатем , подряд,вводятся весовыекоэффициенты.Алгоритм начинаетсяс переписываниявесовых коэффициентовв четвертыйстолбец массиваF(блоки 2,3).Номерастолбцов обозначенынижним индексом,а номера строк- верхним. Послезадания начальныхзначений параметров(4) вызываетсяподпрограммавычисленияпроизводных(5) в начальнойточке интервалаинтегрированияиначальная точкас производнымив ней передаетсяподпрограммевывода (40). Затемначинаетсяосновной циклвыполненияшагов интегрирования.Задается шаг,равный максимальному(6), и выполняютсяшаги из точкиТв точку Т1и из точки Т1в точку Т2.Результатызаписываются,соответственно,во второй итретий столбцымассивов XиF.Затем, для проверкиточности выполняетсяудвоенный шагиз точки Тв точку Т2.Результатыэтого шагазаписываютсяв четвертыйстолбец массиваХи в первый столбецмассива F. В цикле (13, 14) накапливаетсякритерий погрешностиЕН,как сумма взятыхс весами погрешностейпо каждому изуравнений.Погрешностькаждой переменнойвычисляетсякак 1/15модуляразности междузначениямиэтой переменной,вычисленнымис разными шагами.Далее выполняетсяанализ критерия(15) и в зависимостиот его значенияшаг увеличивается,уменьшаетсяили остаетсяпрежним. Еслитекущая погрешностьЕНне больше заданнойЕ, то результатышага выводятся(25). При этом, есливыполнялосьдва малых шага(КLР=1),то выводятсяи результатыпредыдущегошага (23). Так случаетсяв начале интервалаинтегрированияи тогда, когдапредыдущийшаг оказалсянеудачным ииз-за большойпогрешностивеличина шагауменьшенавдвое. Послевывода двухшагов признакKLPсбрасываетсяв ноль (24). Выполненныйшаг может бытьпоследним наинтервале(КР=1),тогда осуществляетсявыход из подпрограммы(26,27). В блоке 28 выполняетсяпроверка, можетли быть выполненудвоенный шагбез выхода запределы интервала? Если нет, то в(29) «взводится»признак концаинтервала иустанавливаетсявеличина удвоенногошага равнойоставшейсячасти интервала.В блоке 33 и цикле34-35 последняявычисленнаяточка делаетсяначальной длявыполнениядвух малыхшагов Ни контрольногоудвоенногоНВ. Соответственно,в 36 устанавливаетсяпризнак двухшагов (KLP=1)иосуществляетсявозврат на блок6.Если «дошагивание»не нужно, то в30 проверяется,является литочность расчетовзавышеннойи в 31 можно лиудвоить малыйшаг?Призавышеннойточности шагможно удвоить,если он не превзойдетмаксимальногоНМиудвоенный шагне выведет запределы интервалаинтегрирования.Если увеличениешага допустимо,то блок 32 этовыполняет идалее всепроизводитсякак при дошагивании,но без взводапризнака конца.Если увеличениешага недопустимо,то в цикле 37, 38выполняетсяподготовкак продолжениюрасчетов спрежним шагом.Из трех последнихточек средняяделается начальнойдля выполненияконтрольногошага удвоеннойвеличины НВ,а последняя, - начальнойдля очередногомалого шагаН.

4.6Блок- схема алгоритмаметода Рунге- Кутта с автоматическимвыбором шага

Рисунок3

4.7Подпрограммаметода Рунге- Кутта с автоматическимвыбором шага

ПодпрограммаARK позволяет решать произвольнуюсистему N-гопорядка савтоматическимвыбором шагаинтегрирования.Эта подпрограммаобращается:

к подпрограмме одного шага-SH,
кподпрограммевычисленияправых частейсистемы,
к подпрограммевывода.

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

Вглавном модулеоператоромEXTERNALдолжны бытьобъявлены именаподпрограммправых частейи вывода. ОператоромDIMENSIONдолжны бытьобъявленымассивы - фактическиепараметрыподпрограммыARK.Эти массивы,по желанию,могут объявлятьсякак одномерныеили как двухмерные.Размеры массивов(N,4),(N,3),(N,4),где N-порядоксистемы. Формальныеимена этихмассивов вподпрограммеARK,соответственно,X,R,F.В главном модулепервые Nэлементовмассива, соответствующегоX,заполняютсяначальнымиусловиями, аследующие Nэлементовзаполняютсявесовымикоэффициентамипогрешности. В подпрограммахправых частейи вывода впервыхNэлементахмассива, соответствующегоX,при входе содержатсятекущие значениявсех Nпеременныхсистемы и недолжны тампереопределяться.Первые Nэлементовмассива, соответствующегоF,должнызаполнятьсяв подпрограммеправых частейвычисляемымитамзначениямиправых частейсистемы.

Формальнымипараметрамив подпрограммеправых частейдолжныбытьпараметры(T,X,F,N),где T-независимаяпеременнаясистемы.

Формальнымипараметрамиподпрограммывывода должныбыть параметры(T,X,F,N,IER),где IER-код ошибки,определяемыйв подпрограммеARK:

IER=0,-ошибки нет;
IER=1,-знак заданногоначальногошага не соответствуетдвижению отначала интервалаинтегрированияк его концу;
IER=2,-начальный шагили/и длинаинтервалаинтегрированияошибочно заданыравными нулю;
IER=3,-шаг в процессесчёта сталболее чем в1000разменьше начального.

МассивыX и Fв подпрограммахправых частейи вывода можнообъявлять какодномерные,с регулируемымразмером X(N),F(N).

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

ПодпрограммаARKвычисляетрешение системыи в каждой точке,удовлетворяющейусловиям точности,обращаетсяк подпрограммевывода, передаваяей значенияпараметровT,X,F,IER.Пользовательможет запрограммироватьздесь печатьнеобходимыхпеременныхили накоплениеих в дополнительныхмассивах дляпоследующейобработки. (Впоследнемслучае дополнительныемассивы следуетпередаватьв главный модульчерез общуюобласть памятис помощью оператораCOMMON). Послевозвратаиз подпрограммывывода,ARKпродолжаетвычислениеследующей точкирешения.

SUBROUTINEARK(HM,TN,TK,X,R,F,N,E,PRAV,OUT,IER)

C Подпрограммаавтоматическоговыбора шага.

C HM-Задаваемыймаксимальныйшаг.

C TN,TK-Начало и конецотрезка интегрирования.

C N-Порядок системы.

C E-Задаваемоезначение абсолютнойпогрешности.

EXTERNALPRAV,OUT

C PRAVи OUT имена составляемыхПользователемподпрограммправых частейи вывода.

C IER-Выходной кодошибки.

DIMENSIONX(N,4),R(N,3),F(N,4)

C Первый столбецмассива Xдолжен привходе содержатьначальныеусловия,

С на выходе внем содержитсярешение.

C Второй столбецмассива Xдолжен привходе содержатьвесовые коэффициентыпогрешности.

C Первый столбецмассива Fдолжен заполнятьсявычисляемыми

C в подпрограммеPRAVзначениямиправых частейсистемы уравнений.

C Остальныеэлементы массивовX,R,F-рабочие.

DO3 K=1,N

3 F(K,4)=X(K,2)

T=TN

HB=2*HM

IER=0

KP=0

KLP=1

CALLPRAV(T,X,F,N)

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

C T-Независимаяпеременнаясистемы.

IF((TK-TN)*HM)4,5,60

4 IER=1

GOTO 60

5 IER=2

60 CALL OUT(T,X,F,N,IER)

C Вызов составленной Пользователем подпрограммывывода результатовшага

IF(IER.NE.0)RETURN

6 H=HB/2

CALLSH(T,H,X,X(1,2),F(1,2),PRAV,N,R)

8 T1=T+H

CALLSH(T1,H,X(1,2),X(1,3),F(1,3),PRAV,N,R)

T2=T+HB

CALLSH(T,HB,X,X(1,4),F,PRAV,N,R)

EH=0

DO14K=1,N

14 EH=EH+ABS((X(K,3)-X(K,4))/15*F(K,4))

IF(EH-E)21,21,16

16 HB=HB/2

IF(HB.LT.HM/512)IER=3

IF(IER.EQ.3)RETURN

KP=0

GOTO 6

21 T=T1

IF(KLP.EQ.1)CALLOUT(T1,X(1,2),F(1,2),N,IER)

KLP=0

CALLOUT(T2,X(1,3),F(1,3),N,IER)

IF(KP.EQ.1)RETURN

IF(ABS(TK-T2)-ABS(HB))29,29,30

29 KP=1

HB=TK-T2

GOTO 33

30 IF(EH-E/50)31,37,37

31 IF(2*H.LE.-HM.AND.ABS(TK-T2).GE.ABS(2*HB))GOTO 32

37 DO38K=1,N

X(K,1)=X(K,2)

38 X(K,2)=X(K,3)

GOTO 8

32HB=2*HB

33 T=T2

DO35K=1,N

35 X(K,1)=X(K,3)

KLP=1

GOTO 6

END

4.8Тестовая задача

Решимдифференциальноеуравнение

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

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

. Вычислим решениена интервале

, что составитпочти

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

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

SUBROUTINEPRAV(T,X,F,N)

DIMENSIONX(N),F(N)

F(1)=X(2)

F(2)=-X(1)

RETURN

END

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

SUBROUTINEOUT(T,X,F,N,IER)

DIMENSIONX(N),F(N),D(3,1000),C(300)

COMMONK,L,KP,D,C

IF(T.LT.31.4)THEN

K=K+1

D(1,K)=T

D(2,K)=X(1)

D(3,K)=X(2)

ENDIF

IF(X(1).LT.0.AND.KP.EQ.1)THEN

L=L+1

C(1)=X(1)

KP=0

ENDIF

IF(X(1).GT.C(L).AND.X(1).GT.0)THEN

C(L)=X(1)

KP=1

ENDIF

IF(T.EQ.270)PRINT*,’T=270’,’ X(270)=’,X(1)

RETURN

END

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

EXTERNALPRAV,OUT

DIMENSIONX(2,4),F(8),R(2,3),D(3,1000),C(300)

COMMONK,L,KP,D,C

READ*,N,TN,TK,HM,((X(K,J),K=1,N),J=1,2),E

K=0

L=1

C(1)=1

KP=1

CALLARK(HM,TN,TK,X,R,F,N,E,PRAV,OUT,IER)

PRINT1, (C(J),J=1,L)

1 FORMAT(I4/(5E15.7))

CALLKRIS(D,3,K,2,0,0.,0.)

END

4.9Результатытестирования

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

Рисунок4

Амплитудыколебаний равныединице, период

Выходнойфайл решенияприведен ниже.

T=270 X(270)= 9.810482E-01

.1000000E+01 .9994009E+00 .9976879E+00 .9948635E+00 .9930583E+00

.9963406E+00 .9985125E+00 .9995713E+00 .9995162E+00 .9983473E+00

.9960660E+00 .9926749E+00 .9945613E+00 .9972748E+00 .9988768E+00

.9993657E+00 .9987408E+00 .9970031E+00 .9941545E+00 .9925186E+00

.9957730E+00 .9979174E+00 .9989495E+00 .9988685E+00 .9976745E+00

.9953687E+00 .9919540E+00 .9940073E+00 .9966935E+00 .9982686E+00

.9987311E+00 .9980807E+00 .9963180E+00 .9934454E+00 .9919787E+00

.9952052E+00 .9973223E+00 .9983279E+00 .9982209E+00 .9970015E+00

.9946712E+00 .9912329E+00 .9934532E+00 .9961117E+00 .1015252E+00

Значениефункции в точкеТ=270 отличаетсяот точногопримерно на0,4% , а положительныемаксимумыотличаютсяот единицы неболее , чем на0,9% . При этом следуетучесть, что вэту погрешностьвошла и погрешностьпроцедурынахождениямаксимума сшагом, равнымшагу интегрирования.Тенденции кзатуханию илираскачиваниюколебаний нет.Все это доказываетработоспособностьалгоритма ипрограммы.

4.10Квадратичнаяконечно-элементнаямодель усилителя

4.10.1 Описаниеалгоритма

Функцияэтого модулязаключаетсяв том, что позаданной входнойвеличине (обозначимее Z3) выдается иливеличина U1,определяемаяиз таблицы 2,или величинаU2,определяемаяиз таблицы 3.Эти таблицыпредставимв виде одногодвухмерногомассива W,в первой строкекоторого запишемтабличныезначения входнойпеременнойZ3, а вовторой и третьейстроках - имсоответствующиетабличныезначения переменныхU1 и U2. Значение ещеодного входногопараметра L,- номерастроки, будетопределять,какую выходнуюпеременнуювычисляетмодель (L=2или L=3). Выходнуюпеременнуюмодуля обозначимU , а длямодуля назначимимя US.Блок - схемаалгоритмаприведена нарисунке 5.

Вцикле с индексомJ определяетсятот конечныйэлемент, в областикоторого находитсявходная величинаZ3 , азатем вычисляетсявыходная величинапо формулеЛагранжа сиспользованиемL-той строкимассива W.

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

4.10.2Блок- схема алгоритмамодели усилителя

Рисунок5

4.10.3Подпрограмма- модель усилителя

SUBROUTINE US(L,Z,U)

С Подпрограмма- модель усилителя.

DIMENSIONW(3,11)

C характеристикиусилителя изтаблиц 2 и 3 постолбцам

DATAW /-3.125 ,-0.125 , 3. ,

= -2.85 , -0.1 , 2.75 ,

= -2.475 , -0.075 , 2.4 ,

= -1.78 , -0.05 , 1.73 ,

= -1.025 ,-0.025 , 1. ,

= -0.02 , 0. , 0.02

= 1.025 , 0.025 , -1. ,

= 1.78 , 0.05 , -1.73 ,

= 2.475 , 0.075 , -2.4 ,

= 2.85 , 0.1 , -2.75 ,

= 3.125 , 0.125 , -3. /

C Поиск интервала,заключающегоZ3.

DO J=2, 10, 2

IF(Z3.GE.W(1,J-1).AND.Z3.LT.W(1,J+1)) GO TO 8

ENDDO

PRINT*,‘ ОШИБКА ‘

STOP

C ФормулаЛагранжа.

8 U=W(L,J-1)*(Z3-W(1,J))*(Z3-W(1,J+1))/

= ((W(1,J-1)- W(1,J))*(W(1,J-1)-W(1,J+1)))+

= W(L,J)*(Z3-W(1,J-1))*(Z3-W(1,J+1))/

= ((W(1,J)-W(1,J-1))*(W(1,J)-W(1,J+1)))+

= W(L,J+1)*(Z3-W(1,J-1))*(Z3-W(1,J))/

= ((W(1,J+1)-W(1,J-1))*(W(1,J+1)-W(1,J)))

RETURN

END

4.10.4Решениетестовой задачи

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

DIMENSIOND(3,1000)

READ*,XN,XK,DX

K=0

DOX=XN,XK,DX

K=K+1

С Вычислениезначения входнойпеременнойU1

CALLUS(2,X,U1)

С Заполнениестроки аргументаU1

D(1,K)=U1

С Вычислениезначения выходнойпеременнойусилителя U2.

CALLUS(3,X,U2)

С Заполнениестроки переменнойU2.

D(2,K)=U2

ENDDO

CALLKRIS(D,3,K, 1, 1,0.,0.)

END

Рисунок6

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

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

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

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

SUBROUTINEFUN(T,Z,F,N)

С Подпрограммавычисленияправых частейсистемы уравнениймодели автогенератора.

DIMENSIONZ(N*4),F(N*4),D(4,15000)

COMMONK,TZ,TAU,D

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

CALLUS(2,Z(3),U1)

A=1/TAU

F(1)=- A*U1

F(2)=A*(Z(1)-5*U1)

F(3)=A*(Z(2)-6*U1)

RETURN

END

4.12Подпрограммавывода

Вподпрограммесохраненынаименованияпеременныхмодели.

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

SUBROUTINEPRIN(T,Z,F,N,IER)

С Подпрограммавывода результатовшага.

DIMENSIONZ(N*4),F(N*4),D(4,15000)

COMMONK,TZ,TAU,D

IF(T.GE.TZ)THEN

K=K+1

С Вычислениезначения переменнойвхода U1.

CALLUS(2,Z(3),U1)

C Вычислениезначения переменнойвыхода U2.

CALLUS(3,Z(3),U2)

С Заполнениемассива.

D(1,K)=T

С Выход усилителябудет изображатьсяна графикахкривой номер1.

D(2,K)=U2

С Вход усилителябудет изображатьсяна графикахкривой номер2.

D(3,K)=U1

ENDIF

RETURN

END

4.13Главныймодуль решениясистемы уравнений

Вглавном модулев соответствиес требованиямиподпрограммыметода Рунге- Кутта ARK объявиммассивы длярешения системытретьего порядка.Имена массивовсохраним такимиже, как именаформальныхпараметровподпрограммыARK. Зададимнулевые начальныеусловия и равныедля всех интегральныхпеременныхвесовые коэффициентыпогрешности.Из исходногофайла будемвводить:

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

,
времяначала интегрированияТN,
времяконца интегрированияТК,
максимальныйшаг интегрированияНМ
задаваемуюпогрешностьЕР.

DIMENSIONZ(12),RAB(9),F(12),D(4,15000)

С Главный модульрешения системыуравнений

EXTERNALFUN,PRIN

COMMONK,TZ,TAU,D

С Задание начальныхусловий и весовыхкоэффициентовпогрешности.

DO1 K=1,3

Z(K)=0.

1 Z(K+3)=0.33333

READ*,TZ,TAU,TN,TK,HM,EP

K=0

С Решение системы.

CALLARK(HM,TN,TK,Z,RAB,F,3,EP,FUN,PRIN,IER)

С Вывод результатовв форме графикови таблиц.

CALLKRIS(D,4,K,2,1,0.,0.)

END

5.ЧИСЛЕННОЕМОДЕЛИРОВАНИЕАВТОГЕНЕРАТОРА

5.1Пробныерешения

Пробноерешение выполнимс параметрами,указаннымив таблице 6

Таблица6

TZ		TN	TK	HM	EP
0	1	0	370	1	0.0001

Рисунок7

Изрисунка видно,что возбуждениеавтогенераторадлится примерно20 периодовколебаний,период колебанияпримерно равен16с., что составляет

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

Таблица7

TZ		TN	TK	HM	EP
370	1	0	400	1	0.0001

Графикирешения приведенына Рисунке 8, ачисленныезначения втаблице 8. Рисунокпоказывает,что выходноенапряжениеавтогенератора(кривая 1) достаточноблизко к синусоидальному,чего нельзясказать о входномнапряженииусилителя(кривая 2).

Таблица8

АРГУМЕНТ ФУНКЦИЯ 1 ФУНКЦИЯ2 ФУНКЦИЯ 3 ФУНКЦИЯ4 ФУНКЦИЯ 5

370.0 -1.753 .5084E-01 .0000

370.5 -1.291 .3469E-01 .0000

371.0 -.7804 .1970E-01 .0000

371.5 -.2281 .6177E-02 .0000

372.0 .3466 -.8225E-02 .0000

372.5 .9243 -.2303E-01 .0000

373.0 1.476 -.4105E-01 .0000

373.5 1.974 -.5888E-01 .0000

374.0 2.395 -.7481E-01 .0000

374.5 2.699 -.9564E-01 .0000

375.0 2.860 -.1103 .0000

375.5 2.885 -.1127 .0000

376.0 2.792 -.1037 .0000

376.5 2.600 -.8794E-01 .0000

377.0 2.324 -.7205E-01 .0000

377.5 1.961 -.5838E-01 .0000

378.0 1.527 -.4280E-01 .0000

378.5 1.038 -.2625E-01 .0000

379.0 .5052 -.1226E-01 .0000

379.5 -.5797E-01 .1948E-02 .0000

380.0 -.6338 .1614E-01 .0000

380.5 -1.202 .3169E-01 .0000

381.0 -1.729 .4996E-01 .0000

381.5 -2.190 .6695E-01 .0000

382.0 -2.559 .8495E-01 .0000

382.5 -2.793 .1038 .0000

383.0 -2.885 .1127 .0000

383.5 -2.849 .1092 .0000

384.0 -2.706 .9619E-01 .0000

384.5 -2.472 .7926E-01 .0000

385.0 -2.152 .6553E-01 .0000

385.5 -1.753 .5082E-01 .0000

386.0 -1.290 .3467E-01 .0000

386.5 -.7795 .1968E-01 .0000

387.0 -.2272 .6154E-02 .0000

387.5 .3476 -.8250E-02 .0000

388.0 .9253 -.2306E-01 .0000

388.5 1.477 -.4108E-01 .0000

389.0 1.975 -.5892E-01 .0000

389.5 2.396 -.7484E-01 .0000

390.0 2.699 -.9568E-01 .0000

390.5 2.861 -.1103 .0000

391.0 2.885 -.1127 .0000

391.5 2.791 -.1037 .0000

392.0 2.600 -.8792E-01 .0000

392.5 2.323 -.7203E-01 .0000

393.0 1.960 -.5836E-01 .0000

393.5 1.526 -.4277E-01 .0000

394.0 1.037 -.2622E-01 .0000

394.5 .5042 -.1223E-01 .0000

395.0 -.5907E-01 .1975E-02 .0000

395.5 -.6350 .1617E-01 .0000

396.0 -1.203 .3172E-01 .0000

396.5 -1.730 .4999E-01 .0000

397.0 -2.191 .6699E-01 .0000

397.5 -2.560 .8500E-01 .0000

398.0 -2.793 .1039 .0000

398.5 -2.885 .1127 .0000

399.0 -2.849 .1091 .0000

399.5 -2.705 .9616E-01 .0000

400.0 -2.472 .7922E-01 .0000

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

Рисунок8

5.2Решениедля спектральногоанализа выходногонапряжения

Выделимодин периодколебаний исделаем третьерешение.

Таблица9

TZ		TN	TK	HM	EP
379,5	1	0	395	1	0.0001

Рисунок9

Таблица 9

АРГУМЕНТ ФУНКЦИЯ 1 ФУНКЦИЯ2 ФУНКЦИЯ 3 ФУНКЦИЯ4 ФУНКЦИЯ 5

379.5 -.5797E-01 .1948E-02 .0000

380.0 -.6338 .1614E-01 .0000

380.5 -1.202 .3169E-01 .0000

381.0 -1.729 .4996E-01 .0000

381.5 -2.190 .6695E-01 .0000

382.0 -2.559 .8495E-01 .0000

382.5 -2.793 .1038 .0000

383.0 -2.885 .1127 .0000

383.5 -2.849 .1092 .0000

384.0 -2.706 .9619E-01 .0000

384.5 -2.472 .7926E-01 .0000

385.0 -2.152 .6553E-01 .0000

385.5 -1.753 .5082E-01 .0000

386.0 -1.290 .3467E-01 .0000

386.5 -.7795 .1968E-01 .0000

387.0 -.2272 .6154E-02 .0000

387.5 .3476 -.8250E-02 .0000

388.0 .9253 -.2306E-01 .0000

388.5 1.477 -.4108E-01 .0000

389.0 1.975 -.5892E-01 .0000

389.5 2.396 -.7484E-01 .0000

390.0 2.699 -.9568E-01 .0000

390.5 2.861 -.1103 .0000

391.0 2.885 -.1127 .0000

391.5 2.791 -.1037 .0000

392.0 2.600 -.8792E-01 .0000

392.5 2.323 -.7203E-01 .0000

393.0 1.960 -.5836E-01 .0000

393.5 1.526 -.4277E-01 .0000

394.0 1.037 -.2622E-01 .0000

394.5 .5042 -.1223E-01 .0000

395.0 -.5907E-01 .1975E-02 .0000

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

Изменяя величину

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

Таблица10

TZ		TN	TK	HM	EP	Т	U_1MAX	U_2MAX	К_УС
370	1	0	400	1	0,0001	15,5	0,1127	2,885	25,6
3200	10	0	3700	10	0,0001	155	0,1127	2,884	25,59
16000	50	0	20000	40	0,0001	780	0,1128	2,886	25,85
32000	100	0	36000	80	0,0001	1560	0,1129	2,886	25,62

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

(17 )

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

6.ПРОГРАММЫОБРАБОТКИРЕЗУЛЬТАТОВМОДЕЛИРОВАНИЯ

6.1Программачисленногоинтегрированияпо методу трапеций

Длявычисленияамплитуды A_n n-ойгармоники

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

Функция

на периоде

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

величину

на значение

. Применяя формулутрапеций, интегралзаменим суммой

(18 )

гдеМ=33 ,- количествоточек в таблице9.

Тогдаамплитуду n-ойгармоникиможно вычислить,как

(19 )

Вычислимв цикле амплитудыдевяти гармоники занесем ихв массивD дляпостроенияграфика с помощьюподпрограммыKRIS.

Блок- схема и программавычисленияамплитуд гармоникприведены ниже.

DIMENSIONT(200),U2(200),F(200),A(9),D(2,9)

READ*,M,L,(T(K),U2(K),X,Y,K=1,M)

DON=1,9

DOK=1,M,L

F(K)=U2(K)*SIN(N*0.405366*T(K))

ENDDO

S=0

DOK=1,M-1,L

S=S+(T(K+1)-T(K))*(F(K)+F(K+1))

ENDDO

A(N)=S/15.5

D(1,N)=N

D(2,N)=A(N)

ENDDO

CALLKRIS(D,2,9,1,0,0.,0.)

PRINT16,(N,A(N),N=1,9)

16 FORMAT(I4,E14.6)

END

Изменениешага L позволяетоценить погрешностьинтегрирования.ПеременныеX и Yнужны в спискеввода для считыванияданных прямоиз выходногофайла третьегорешения.

6.2Блок -схема алгоритмавычисленияамплитуд гармоник

Рисунок10

6.3Результатыгармоническогоанализа

Зависимостьамплитудыгармоники отее номера приведеныв таблицах 11,12 и на рисунке11.

Таблица 11

1 .284373E+01

2 .222451E-02

3 .103735E-01

4 .498333E-03

5 -.751302E-02

6 .191248E-03

7 .318412E-02

8 -.107523E-04

9 .145544E-03

Рисунок11

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

Таблица 12

1 .284373E+01

2 .222451E-02

3 .103735E-01

4 .498333E-03

5 -.751302E-02

6 .191248E-03

7 .318412E-02

8 -.107523E-04

9 .145544E-03

Интегрированиепроведено свысокой точностью,так как обарешения совпадают.

Четныегармоникипрактическиравны нулю, анаибольшаяиз нечетных,- третья составляетвсего 0,36%от первой. Втаких условияхаппроксимацияэтой характеристикине имеет смысла.

7.ЛИТЕРАТУРА

Б.П.ДЕМИДОВИЧ, И.А. МАРОН, Основывычислительнойматематики,«Наука», М., 1966.
Б.П.ДЕМИДОВИЧ, И.А. МАРОН, Э.З.ШУВАЛОВА, Численныеметоды анализа,«Наука», М., 1967.
И.С.БЕРЕЗИН, Н.П.ЖИДКОВ, Методывычислений,Физматгиз,1961.
Н.Н.КАЛИТКИН, Численныеметоды, «Наука»,М., 1978.
Н.С.БАХВАЛОВ, Численныеметоды, «Наука»,М., 1975.
Д.ХИММЕЛЬБЛАУ,Прикладноенелинейноепрограммирование,«Мир», М., 1975.
А.А.ФЕЛЬДБАУМ,А.Д. ДУДЫКИН,А.П. МАНОВЦЕВ,Н.Н. МИРОЛЮБОВ,Теоретическиеосновы связии управления,Физматгиз, М.,1963.
З.С.БРИЧ, Д.В. КАПИЛЕВИЧ,Н.А. КЛЕЦКОВА,ФОРТРАН 77 дляПЭВМ ЕС, «Финансыи статистика»,М., 1991.
П.В.СОЛОВЬЕВ, FORTRANдляперсональногокомпьютера,«ARIST»,М., 1991.
Г.Н.РЫБАЛЬЧЕНКО,Численныеметоды решениязадач строительствана ЭВМ, КиевУМК ВО, 1989.
Г.Н. РЫБАЛЬЧЕНКО,Методическиеуказания ккурсовой работепо дисциплине«Основы вычислительнойматематики»,Кривой Рог,КТУ, 1997.

ИсследованиеRC-генераторасинусоидальныхколебаний.