Исследование 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
U1 | -0,125 | -0,1 | -0,075 | -0,05 | -0,025 | 0 | 0,025 | 0,05 | 0,075 | 0,1 | 0,125 |
U2 | 3 | 2,75 | 2,4 | 1,73 | 1 | 0,02 | -1 | -1,73 | -2,4 | -2,75 | -3 |
Численнымиэкспериментамина ЭВМ найтизависимости:
-
периодаТустановившихсяавтоколебанийот параметра
, -
амплитудыU2maxвыходногонапряженияU2(t)от
, -
амплитудыAn 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
z3 | -3,125 | -2,85 | -2,475 | -1,78 | -1,025 | -0,02 | 1,025 | 1,78 | 2,475 | 2,85 | 3,125 |
U1 | -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
z3 | -3,125 | -2,85 | -2,475 | -1,78 | -1,025 | -0,02 | 1,025 | 1,78 | 2,475 | 2,85 | 3,125 |
U2 | 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Блок- схема алгоритмаметода Рунге- Кутта с автоматическимвыбором шага
Рисунок34.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
0
.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.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
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
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 | Т | U1MAX | U2MAX | КУС |
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Программачисленногоинтегрированияпо методу трапеций
Длявычисленияамплитуды An n-ойгармоники
выходногонапряжения отее номера nнеобходимонесколько развычислятьопределенныйинтеграл ,Функция
на периоде 
вычислена намии представленав таблице 9.Подынтегральнуюфункцию получим,умножая в каждойточке таблицы 
величину 
на значение 
. Применяя формулутрапеций, интегралзаменим суммой 
(18 )гдеМ=33 ,- количествоточек в таблице9.
Тогдаамплитуду n-ойгармоникиможно вычислить,как
Вычислимв цикле амплитудыдевяти гармоники занесем ихв массив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-генераторасинусоидальныхколебаний.