Дослідження RC-генератора синусоїдальних коливань

[pic].

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 Рішення задля встановлення залежностей параметрів від [pic] 25.

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

6.1 Програма чисельного інтегрування методом трапецій 26 6.2 Блок — схема алгоритму обчислення амплітуд гармонік 27 6.3 Результати гармонійного аналізу 28.

7. ЛІТЕРАТУРА 29.

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

Виконати дослідження RC-генератора синусоидальных коливань (Малюнок. 1).

[pic].

Малюнок 1.

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

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

[pic], де [pic].

Нелінійна залежність вихідного напруги [pic]усилителя з його вхідного напруги [pic] приведено в таблиці 1.

Таблиця 1.

|U1 |-0,12|-0,1|-0,07|-0,0|-0,02|0 |0,02|0,05|0,07|0,1 |0,12| | |5 | |5 |5 |5 | |5 | |5 | |5 | |U2 |3 |2,75|2,4 |1,73|1 |0,02|-1 |-1,7|-2,4|-2,7|-3 | | | | | | | | | |3 | |5 | |.

Численными експериментами на ЕОМ знайти залежності: періоду Т встановлених автоколебаний від параметра [pic], амплітуди U2max вихідного напруги U2(t) від [pic], амплітуди An n-ой гармоніки [pic] вихідного напруги від неї номери n.

[pic], коефіцієнта посилення електронного підсилювача як встановлених автоколебаний від [pic]. Знайдені експериментально залежності апроксимувати степенными многочленами.

З залежності [pic] знайти значення [pic], необхідне отримання періоду автоколебаний [pic], і розрахунком коливань перевірити правильність отриманого значення параметра [pic].

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

МАТЕМАТИЧНА ФОРМУЛЮВАННЯ ЗАДАЧИ.

1 ДИФЕРЕНЦІАЛЬНІ УРАВНЕНИЯ ЛІНІЙНОЇ ЧАСТИ.

Запишемо систему диференційних рівнянь лінійної частини RCгенератора. І тому перетворимо її передатну функцию.

[pic] (1).

[pic] (2).

Введемо першу допоміжну зміну [pic], котру визначаємо з уравнения.

[pic].

(3).

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

[pic] (4).

Скорочуючи на [pic] і групуючи у правій частині члени, які містять [pic], получаем.

[pic] (5).

Введемо другу допоміжну зміну [pic], котру визначаємо з уравнения.

[pic].

(6).

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

[pic] (7).

Знову скорочуючи на [pic] і групуючи у правій частині члени, не містять [pic], получаем.

[pic] (8).

Введемо третю допоміжну зміну [pic], котру визначаємо з уравнения.

[pic].

(9).

Підставляючи (9) в (8) і скорочуючи на [pic], получаем.

[pic].

(10).

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

[pic].

(11).

[pic].

(12).

[pic].

(13).

[pic].

(14).

Тут [pic] - функція, обумовлена нелінійної характеристикою усилителя.

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

2 Рівняння усилителя.

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

Можна вирішувати це рівняння методом ітерацій. Але є простіший путь.

Знайдемо з характеристики підсилювача різниці [pic], та був побудуємо характеристику [pic] Значення [pic] відомо спочатку з початкових умов, та був при кожному зверненні до вирахування правих частин системи та з побудованої нами характеристики можна обчислити [pic] для підстановки в праві частини інших уравнений.

Розрахований характеристика представленій у таблиці 2.

Таблиця 2.

|z3 |-3,12|-2,8|-2,47|-1,7|-1,02|-0,02|1,025|1,78|2,475|2,85|3,12| | |5 |5 |5 |8 |5 | | | | | |5 | |U1 |-0,12|-0,1|-0,07|-0,0|-0,02|0 |0,025|0,05|0,075|0,1 |0,12| | |5 | |5 |5 |5 | | | | | |5 |.

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

Для побудови квадратичного кінцевого елемента використовуємо интерполяционную формулу Лагранжа.

[pic] (15).

Для обчислення вихідний величини автогенератора слід також по формулі Лагранжа по заданому значенням [pic] знаходити [pic]. [pic] (16) Дані у разі необхідно вибирати з таблиці 3, отриманою з таблиць 1 і 2.

Таблиця 3.

|z3 |-3,12|-2,8|-2,47|-1,7|-1,02|-0,02|1,025|1,78 |2,47|2,85|3,125| | |5 |5 |5 |8 |5 | | | |5 | | | |U2 |3 |2,75|2,4 |1,73|1 |0,02 |-1 |-1,73|-2,4|-2,7|-3 | | | | | | | | | | | |5 | |.

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

Функціонально програмний комплекс має складатися з двох незалежних частин: програми — моделі RC — генератора; набору програм обробки результатів моделювання автогенератора.

Модель RC — генератора повинна, своєю чергою, включати: модуль, викликає підпрограму методу Рунге — Кутта; модулі методу Рунге — Кутта; модуль — модель підсилювача; модуль правих частин; модуль виведення результатів один крок интегрирования.

Для програмної реалізації методу Рунге — Кутта зручно використовувати два модуля: модуль, виконує один поставлене крок методу; модуль, управляючий величиною кроку залежність від одержуваної похибки решения.

Взаємодія цих модулів таке. Зухвалий модуль вводить значення параметра [pic], початок і поклала край інтервалу інтегрування, максимальний крок, початкові умови і задану похибка. Далі ця модуль звертається до модуля управління методу Рунге — Кутта. Останній задає величину кроку подпрограмме один крок і призводить процес інтегрування системи рівнянь, утримуючи похибка в заданих межах. При виконання кроку, в відповідність до методом Рунге — Кутта, модуль кроку чотири рази звертається до модулю правих частин, а той, своєю чергою, — до моделі підсилювача як функції [pic]. По виконанні кроку, задовольняючого умовам точності, модуль управління викликає підпрограму виведення результатів кроку, а вона, в своє чергу звертається до моделі підсилювача як функції [pic]. Модуль управління завершує своє роботу після досягнення кінця інтервалу інтегрування. Тоді викликає модуль звертається до подпрограмме виведення таблиць і графіків KRIS.

У набір підпрограм обробки результатів моделювання необхідно включити дві незалежні програми: програму чисельного інтегрування методом трапецій; програму апроксимації експериментальних залежностей степенными многочленами методом найменших квадратов.

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

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

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

Нехай є уравнение.

[pic] чи [pic] де [pic] - невідома функція від незалежної перемінної [pic]; [pic] - відома функція. Усі чисельні на методи вирішення завдання Коші засновані на наближеною заміні шуканої функції степенными многочленами.

У методі Рунге-Кутта четвертого порядку відшукується прирощення, яке дає що наближує багаточлен на кроці інтегрування. Прирощення шуканої функції обчислюється як твори довжини кроку на значення похідною від цього функції. Як похідною береться середньозважене від значень похідних [pic] вирахуваних в спеціально підібраних чотирьох точках.

Як першу точки беруть початкову точку кроку [pic]. Похідна у цій точці равна.

[pic], де [pic] - права частина рівняння .

Як другу крапки над площині рішення [pic] вибирають точку з координатами [pic]. Похідна на другий точці равна.

[pic] Для третьої точки беруть координати [pic] і обчислюють производную.

[pic] Нарешті, для четвертої точки беруть координати [pic] і обчислюють производную.

[pic] За отриманими чотирьох значенням похідною знаходять середньозважене значение.

[pic] Тепер, знаходять координати кінцевої точки шага.

[pic].

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

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 і одну зберігають у першому стовпці масиву Р.

Цикл 8 — 11 записує у стовпець робочого масиву значення змінних для обчислення похідних. І тому до початковому значенням перемінної додається спочатку нульовий прирощення, потім половина збільшення, одержуваного на кроці багатозначно похідною у початковій точці, потім половина збільшення, одержуваного на кроці з значенням похідною у другий точці, і, нарешті, повне прирощення, одержуване на кроці зі значенням похідною у третій точке.

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

У циклі 13 — 16 на другий стовпець робочого масиву по вичисленим значенням похідних записуються збільшення, а третьому стовпці накопичуються суми чотирьох збільшень з ваговими коефіцієнтами 1/6, 1/3, 1/3, 1/6. Константи 6, 3, 3, 6 при цьому зберігаються в другому стовпці масиву Р.

У циклі 19 — 22 отримані збільшення додаються до початковій точці і результати записується в вихідний массив.

У блоці 23 обчислюються похідні в кінцевої точки шага.

3 Блок — схема алгоритму один крок методом Рунге — Кутта.

[pic].

Малюнок 2.

4 Підпрограма один крок методом Рунге-Кутта.

SUBROUTINE SH (TN, H, XN, XK, F, PRAV, N, R).

DIMENSION XN (N), XK (N), F (N), P (4,2), R (N, 3).

DATA P/0.,.5,.5,1., 6., 3., 3., 6./.

DO 1 K=1,N.

R (K, 2)=0.

1 R (K, 3)=0.

DO 4 J=1,4 ! Початок зовнішнього циклу визначення 4-х приращений.

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

DO 2 K=1,N ! Цикл завдання масиву значень залежних переменных.

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

CALL PRAV (T, R, F, N) ! Обчислення похідних в заданої точке.

DO 3 K=1,N ! Цикл обчислення та накопичення грошових часткових приращений.

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

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

4 CONTINUE.

DO 5 K=1,N.

5 XK (K)=XN (K)+R (K, 3) ! Обчислення змінних в кінцевої точке.

CALL PRAV (TN+H, XK, F, N) ! Обчислення похідних в кінцевої точке.

RETURN.

END.

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 виконується підготовка продовження розрахунків із колишнім кроком. Із трьох останніх точок середня робиться початковій до виконання контрольного кроку подвоєною величини НВ, а остання , — початковій чергового малого кроку Н.

6 Блок — схема алгоритму методу Рунге — Кутта з автоматичним вибором шага.

[pic] Малюнок 3.

7 Підпрограма методу Рунге — Кутта з автоматичним вибором шага.

Підпрограма ARK дозволяє вирішувати довільну систему N-го порядку з автоматичним вибором кроку інтегрування. Ця підпрограма звертається: до подпрограмме один крокSH, до подпрограмме обчислення правих частин системи, до подпрограмме вывода.

Підпрограма SH записана в універсальному вигляді й наведена вище. головний викликає модуль, і навіть підпрограми правих частин 17-ї та виведення Користувач має скласти самостоятельно.

У головному модулі оператором EXTERNAL повинні прагнути бути оголошено імена підпрограм правих частин 17-ї та виведення. Оператором DIMENSION мали бути зацікавленими оголошено масиви — фактичні параметри підпрограми ARK. Ці масиви, по бажанню, можуть оголошуватися як одномірні чи як двухмерные. Розміри масивів (N, 4); INSERT INTO `ref` (`id_predmet`, `name_predmet`, `id_ref`, `name_ref`, `text_ref`) VALUES (N, 3); INSERT INTO `ref` (`id_predmet`, `name_predmet`, `id_ref`, `name_ref`, `text_ref`) VALUES (N, 4), де N-порядок системи. Формальні імена цих масивів в подпрограмме ARK, відповідно, X, R, F. У головному модулі перші N елементів масиву, відповідного X, заповнюються початковими умовами, а такі N елементів заповнюються ваговими коефіцієнтами похибки. У підпрограмах правих частин 17-ї та виведення у перших 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 в підпрограмах правих частин 17-ї та виведення можна оголошувати як одномірні, із регульованим розміром X (N), F (N).

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

Підпрограма ARK обчислює рішення системи та у кожному точці, задовольняє умовам точності, звертається до подпрограмме виведення, передаючи їй значення параметрів T, X, F, IER. Користувач може запрограмувати тут печатку необхідних змінних й нагромадження в додаткових масивах для наступної обробки. (У разі додаткові масиви слід передавати їх у головний модуль через загальну область пам’яті з допомогою оператора COMMON). Після повернення з підпрограми виведення, ARK продовжує обчислення наступній точки решения.

SUBROUTINE ARK (HM, TN, TK, X, R, F, N, E, PRAV, OUT, IER) З Підпрограма автоматичного вибору кроку. З HMЗадаваемый максимальний крок. З TN, TKПочаток і поклала край відрізка інтегрування. З NПорядок системи. З EЗадаваемое значення абсолютної погрешности.

EXTERNAL PRAV, OUT З PRAV і OUT імена який складають Користувачем підпрограм правих частин 17-ї та виведення. З IERВихідний код ошибки.

DIMENSION X (N, 4), R (N, 3), F (N, 4) З Перший стовпець масиву X повинен перед входом утримувати початкові умови, З не вдома у ньому є рішення. З Другий стовпець масиву X повинен перед входом утримувати вагові коефіцієнти похибки. З Перший стовпець масиву F повинен заповнюватися вычисляемыми З в подпрограмме PRAV значеннями правих частин системи рівнянь. З Інші елементи масивів X, R, Fрабочие.

DO 3 K=1,N.

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

T=TN.

HB=2*HM.

IER=0.

KP=0.

KLP=1.

CALL PRAV (T, X, F, N) З Виклик складеної Користувачем підпрограми правих частин системи рівнянь. З TНезалежна змінна системы.

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

4 IER=1.

GO TO 60.

5 IER=2 60 CALL OUT (T, X, F, N, IER) З Виклик складеної Користувачем підпрограми виведення результатів шага.

IF (IER.NE.0)RETURN.

6 H=HB/2.

CALL SH (T, H, X, X (1,2), F (1,2), PRAV, N, R).

8 T1=T+H.

CALL SH (T1,H, X (1,2), X (1,3), F (1,3), PRAV, N, R).

T2=T+HB.

CALL SH (T, HB, X, X (1,4), F, PRAV, N, R).

EH=0.

DO 14 K=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.

GO TO 6 21 T=T1.

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

KLP=0.

CALL OUT (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.

GO TO 33 30 IF (EH-E/50)31,37,37 31 IF (2*H.LE.-HM.AND.ABS (TK-T2).GE.ABS (2*HB))GO TO 32 37 DO 38 K=1,N.

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

GO TO 8.

32 HB=2*HB.

33 T=T2.

DO 35 K=1,N.

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

KLP=1.

GO TO 6.

END.

8 Тестова задача.

Вирішимо диференціальний уравнение.

[pic] з початковими умовами [pic]. Легко бачити, що рішенням це завдання є функція [pic]. Обчислимо рішення на інтервалі [pic], що становитиме майже [pic] періоду цієї функції. Якщо за такому тривалому інтегруванні амплітуда косинусоиды істотно не зміниться, то алгоритм чисельно стійкий. Можна ще порівняти рішення, у кінцевої точки [pic].

Підпрограма правих частин при цьому рівняння буде такой.

SUBROUTINE PRAV (T, X, F, N).

DIMENSION X (N), F (N).

F (1)=X (2).

F (2)= -X (1).

RETURN.

END.

У подпрограмме виведення передбачимо заповнення результатами масиву D для побудови графіків на інтервалі до п’яти періодів, і навіть заповнення масиву З позитивними максимумами вычисляемой функції по всьому інтервалі інтегрування. Ці масиви передамо у Ганно-Леонтовичевому чільний модуль через загальну область.

SUBROUTINE OUT (T, X, F, N, IER).

DIMENSION X (N), F (N), D (3,1000), C (300).

COMMON K, 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.

У головному модулі введемо вихідні дані, звернімося подпрограмме методу, отпечатаем отримані через загальну область максимуми функції і звернімося подпрограмме побудови графика.

EXTERNAL PRAV, OUT.

DIMENSION X (2,4), F (8), R (2,3), D (3,1000), C (300).

COMMON K, 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.

CALL ARK (HM, TN, TK, X, R, F, N, E, PRAV, OUT, IER).

PRINT 1, (C (J), J=1,L) 1 FORMAT (I4/(5E15.7)).

CALL KRIS (D, 3, K, 2,0,0., 0.).

END.

9 Результати тестирования.

Графіки вирахуваних вирішенням диференціального рівняння функцій наведено малюнку 4. Очевидно, що вони близькі функцій [pic] і [pic].

[pic].

Малюнок 4.

Амплітуди коливань рівні одиниці, період [pic] .

Вихідний файл рішення наведено нижче. T=270 X (270)= 9.81 0482E-01 0 .100 0000E+01 .999 4009E+00 .997 6879E+00 .994 8635E+00 .993 0583E+00 .996 3406E+00 .998 5125E+00 .999 5713E+00 .999 5162E+00 .998 3473E+00 .996 0660E+00 .992 6749E+00 .994 5613E+00 .997 2748E+00 .998 8768E+00 .999 3657E+00 .998 7408E+00 .997 0031E+00 .994 1545E+00 .992 5186E+00 .995 7730E+00 .997 9174E+00 .998 9495E+00 .998 8685E+00 .997 6745E+00 .995 3687E+00 .991 9540E+00 .994 0073E+00 .996 6935E+00 .998 2686E+00 .998 7311E+00 .998 0807E+00 .996 3180E+00 .993 4454E+00 .991 9787E+00 .995 2052E+00 .997 3223E+00 .998 3279E+00 .998 2209E+00 .997 0015E+00 .994 6712E+00 .991 2329E+00 .993 4532E+00 .996 1117E+00 .101 5252E+00.

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

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

1 Опис алгоритма.

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

У циклі з індексом J визначається той кінцевий елемент, у сфері якої перебуває вхідні величина Z3, та був обчислюється вихідна величина за такою формулою Лагранжа з допомогою L-той рядки масиву W. [pic].

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

2 Блок — схема алгоритму моделі усилителя.

[pic].

Малюнок 5.

3 Підпрограма — модель усилителя.

SUBROUTINE US (L, Z, U).

З Підпрограма — модель усилителя.

DIMENSION W (3,11).

З характеристики підсилювача з таблиць 2 і трьох по столбцам.

DATA W /-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. /.

З Пошук інтервалу, заключающего 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.

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

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 Рішення тестової задачи.

Як тестової завдання обчислимо малим кроком, і побудуємо графіки характеристик усилителя.

DIMENSION D (3,1000).

READ*, XN, XK, DX.

K=0.

DO X=XN, XK, DX.

K=K+1.

З Обчислення значення вхідний перемінної U1.

CALL US (2,X, U1).

З Заповнення рядки аргументу U1.

D (1,K)=U1.

З Обчислення значення вихідний перемінної підсилювача U2.

CALL US (3,X, U2).

З Заповнення рядки перемінної U2.

D (2,K)=U2.

ENDDO.

CALL KRIS (D, 3, K, 1, 1,0., 0.).

END.

[pic].

Малюнок 6.

З малюнка видно, що характеристика підсилювача відтворюється моделлю правильно.

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

У подпрограмме збережені найменування змінних моделі. Результати обчислень заносяться, як і вимагають підпрограми кроку та управляючого модуля, у стовпець масиву F, що тут, для простоти оголошено одномірною. Для передачі у цей модуль змінюваного від експерименту до експерименту параметра генератора [pic] у загальну область включена змінна TAU .Інші перемінні загальної області потрібні для зв’язку головного модуля з підпрограмою виведення результатів шага.

SUBROUTINE FUN (T, Z, F, N) З Підпрограма обчислення правих частин системи рівнянь моделі автогенератора.

DIMENSION Z (N*4), F (N*4), D (4,15 000).

COMMON K, TZ, TAU, D З Виклик підпрограми — моделі підсилювача для обчислення вхідний величини U1.

CALL US (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.

12 Підпрограма вывода.

У подпрограмме збережені найменування змінних модели.

А, щоб матимуть можливість хоча б якісно, але швидко, оцінювати правильність роботи моделі потрібно здійснити візуалізацію рішення. Тож у модулі виведення кожен крок обчислимо вхідну і вихідну перемінні підсилювача і заповнимо цими даними черговий стовпець масиву D. У цей самий стовпець запишемо поточні значення часу Т. Масив D передамо через загальну область у Ганно-Леонтовичевому чільний модуль, а звідти подпрограмме побудови графіків KRIS. У автогенераторе кілька днів триває процес самозбудження. Нас цікавить процес встановлених коливань, тому запис даних в масив робитимемо лише починаючи з певного моменту часу TZ. Ця величина і лічильник точок також ввімкнемо у загальну область.

SUBROUTINE PRIN (T, Z, F, N, IER) З Підпрограма виведення результатів шага.

DIMENSION Z (N*4), F (N*4), D (4,15 000).

COMMON K, TZ, TAU, D.

IF (T.GE.TZ)THEN.

K=K+1 З Обчислення значення перемінної входу U1.

CALL US (2,Z (3), U1) З Обчислення значення перемінної виходу U2.

CALL US (3,Z (3), U2) З Заповнення массива.

D (1,K)=T З Вихід підсилювача буде зображуватися на графіках кривою номер 1.

D (2,K)=U2 З Вхід підсилювача буде зображуватися на графіках кривою номер 2.

D (3,K)=U1.

ENDIF.

RETURN.

END.

13 Головний модуль рішення системи уравнений.

У головному модулі у відповідність із вимогами підпрограми методу Рунге — Кутта ARK оголосимо масиви на вирішення системи третього порядку. Імена масивів збережемо так само, як імена формальних параметрів підпрограми ARK. Поставимо нульові початкові умови й однакові всім інтегральних змінних вагові коефіцієнти похибки. З вихідного файла вводитимемо: час початку записи даних в вихідний масив TZ, параметр [pic], час початку інтегрування ТN, час кінця інтегрування ТК, максимальний крок інтегрування ПМ задаваемую похибка ЕР.

DIMENSION Z (12), RAB (9), F (12), D (4,15 000) З Головний модуль рішення системи уравнений.

EXTERNAL FUN, PRIN.

COMMON K, TZ, TAU, D З Завдання початкових умов і вагових коефіцієнтів погрешности.

DO 1 K=1,3.

Z (K)=0.

1 Z (K+3)=0.33 333.

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

K=0 З Рішення системы.

CALL ARK (HM, TN, TK, Z, RAB, F,3,EP, FUN, PRIN, IER) З Висновок успіхів у формі графіків і таблиц.

CALL KRIS (D, 4, K, 2,1,0., 0.).

END.

ЧИСЕЛЬНА МОДЕЛЮВАННЯ АВТОГЕНЕРАТОРА.

1 Пробні решения.

Пробне рішення виконаємо з параметрами, зазначеними в таблиці 6.

Таблиця 6 |TZ |[pi|TN |TK |HM |EР | | |з] | | | | | |0 |1 |0 |370 |1 |0.0001|.

[pic].

Малюнок 7.

З малюнка видно, що порушення автогенератора триває приблизно 20 періодів коливань, період коливання приблизно дорівнює 16с., що становить [pic].

Друге рішення виконаємо те щоб запис почалася режимі встановлених коливань і тривала близько двох періодів. Тоді таблиці рішення з достатньої точністю встановити амплітуду і період коливань. Дані на другому рішення наведені у таблиці 7.

Таблиця 7 |TZ |[pi|TN |TK |HM |EР | | |з] | | | | | |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.

[pic].

Малюнок 8.

2 Рішення для спектрального аналізу вихідного напряжения.

Виділимо один період коливань і зробимо третє решение.

Таблиця 9 |TZ |[pi|TN |TK |HM |EР | | |з] | | | | | |379,5|1 |0 |395|1 |0.000| | | | | | |1 |.

[pic].

Малюнок 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.

3 Рішення задля встановлення залежностей параметрів від [pic].

Змінюючи величину [pic], робимо рішення, аналогічні другому, і результати, витягнуті у вихідних файлів, заносимо в таблицю 10.

Таблиця 10.

|TZ |[pic|TN |TK |HM |EР |Т |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| |16 000|50 |0 |20 000|40 |0,0001|780 |0,1128|2,886|25,85| |32 000|100 |0 |36 000|80 |0,0001|1560 |0,1129|2,886|25,62|.

Аналізуючи цих результатів, доходимо висновку, що період коливань пропорційний [pic].

[pic].

(17).

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

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

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

Для обчислення амплітуди An n-ой гармоніки [pic] вихідного напруги від неї номери n необхідно кілька разів обраховувати певний інтеграл [pic],.

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

[pic].

(18) де М=33 , — кількість точок в таблиці 9.

Тоді амплітуду n-ой гармоніки можна визначити, как.

[pic].

(19).

Обчислимо в циклі амплітуди дев’яти гармонік і занесемо в масив D для побудови графіка з допомогою підпрограми KRIS.

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

DIMENSION T (200), U2(200), F (200), A (9), D (2,9).

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

DO N=1,9.

DO K=1,M, L.

F (K)=U2(K)*SIN (N*0.405 366*T (K)).

ENDDO.

S=0.

DO K=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.

CALL KRIS (D, 2,9,1,0,0., 0.).

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

16 FORMAT (I4,E14.6).

END.

Зміна кроку L дозволяє оцінити похибка інтегрування. Змінні X і Y потрібні списку введення для зчитування даних безпосередньо з вихідного файла третього решения.

2 Блок — схема алгоритму обчислення амплітуд гармоник.

[pic].

Малюнок 10.

3 Результати гармонійного анализа.

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

Таблиця 11.

1 .28 4373E+01.

2 .22 2451E-02.

3 .10 3735E-01.

4 .49 8333E-03.

5 -.75 1302E-02.

6 .19 1248E-03.

7 .31 8412E-02.

8 -.10 7523E-04.

9 .14 5544E-03.

[pic].

Малюнок 11.

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

Таблиця 12.

1 .28 4373E+01.

2 .22 2451E-02.

3 .10 3735E-01.

4 .49 8333E-03.

5 -.75 1302E-02.

6 .19 1248E-03.

7 .31 8412E-02.

8 -.10 7523E-04.

9 .14 5544E-03.

Інтегрування проведено з точністю, оскільки обидва рішення совпадают.

Парні гармоніки практично рівні нулю, а найбільша з непарних, — третя не перевищує 0,36% з першої. За цих умов апроксимація цієї характеристики немає смысла.

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