СОДЕРЖАНИЕ
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
Имя | Тип | Назначение |
PRAV | Подпрограмма. | Подпрограмма, |
N | Целый. | Порядок |
XN | Вещественный. | Исходный |
XK | Вещественный. | Результирующий |
F | Вещественный. | Массив |
TN | Вещественный. | Начальное |
K | Целый. | Номер |
J | Целый. | Номер |
T | Вещественный. | Независимая |
H | Вещественный. | Задаваемая |
P | Вещественный | Массив |
R | Вещественный | Рабочий |
Блок-схема
алгоритма
изображена
на Рисунке 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 вычисляются
производные
в конечной
точке шага.
- схема алгоритма
одного шага
по методу Рунге
- Кутта
Рисунок
2
4.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
4.5Описаниеалгоритма
метода Рунге
- Кутта с автоматическим
выбором шага
Блок
-схема алгоритма
приведена на
Рисунке 3.
В
алгоритме
используются
следующие
идентификаторы
Таблица
5
Имя | Тип | Назначение |
PRAV | Подпрограмма. | Подпрограмма, |
OUT | Подпрограмма. | Подпрограмма, |
N | Целый. | Порядок |
X | Вещественный. | Массив |
R | Вещественный. | Рабочий |
F | Вещественный. | Массив |
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
метода Рунге
- Кутта с автоматическим
выбором шага
Подпрограмма
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
продолжает
вычисление
следующей точки
решения.
SUBROUTINE
ARK(HM,TN,TK,X,R,F,N,E,PRAV,OUT,IER)
C
Подпрограм
автоматического
выбора шага.
C
HM
-Задаваемый
максимальный
шаг.
C
TN,TK
-Начало и конец
отрезка интегрирования.
C
N
-Порядок системы.
C
E
-Задаваемое
значение абсолютной
погрешности.
EXTERNAL
PRAV,OUT
C
PRAV
и OUT
имена составляемых
Пользователем
подпрограмм
правых частей
и вывода.
C
IER
-Выходной код
ошибки.
DIMENSION
X(N,4),R(N,3),F(N,4)
C
Первый столбец
массива X
должен при
входе содержать
начальные
условия,
С
на выходе в
нем содержится
решение.
C
Второй столбец
массива X
должен при
входе содержать
весовые коэффициенты
погрешности.
C
Первый столбец
массива F
должен заполняться
вычисляемыми
C
в подпрограмме
PRAV
значениями
правых частей
системы уравнений.
C
Остальные
элементы массивов
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)
C
Вызов составленной
Пользователем
подпрограммы
правых частей
системы уравнений.
C
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)
C
Вызов составленной
Пользователем
подпрограммы
вывода результатов
шага
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
задача
Решим
дифференциальное
уравнение
с
начальными
условиями
. Легко видеть,
что решением
этой задачи
является функция
. Вычислим решение
на интервале
, что составит
почти
периода этой
функции. Если
при таком длительном
интегрировании
амплитуда
косинусоиды
существенно
не изменится,
то алгоритм
численно устойчив.
Можно также
сравнить решение
в конечной
точке
Подпрограмма
правых частей
для этого уравнения
будет такой.
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
тестирования
Графики
вычисленных
путем решения
дифференциального
уравнения
функций приведены
на рисунке 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.1 Описание
алгоритма
Функция
этого модуля
заключается
в том, что по
заданной входной
величине (обозначим
ее Z3
) выдается или
величина U1,
определяемая
из таблицы 2,
или величина
U2,
определяемая
из таблицы 3.
Эти таблицы
представим
в виде одного
двухмерного
массива W,
в первой строке
которого запишем
табличные
значения входной
переменной
Z3, а во
второй и третьей
строках - им
соответствующие
табличные
значения переменных
U1 и U2
. Значение еще
одного входного
параметра L
,- номера
строки, будет
определять,
какую выходную
переменную
вычисляет
модель (L=2
или L=3). Выходную
переменную
модуля обозначим
U , а для
модуля назначим
имя US.
Блок - схема
алгоритма
приведена на
рисунке 5.
В
цикле с индексом
J определяется
тот конечный
элемент, в области
которого находится
входная величина
Z3 , а
затем вычисляется
выходная величина
по формуле
Лагранжа с
использованием
L-той строки
массива W.
Если
значение входной
переменной
Z3 выходит
за пределы
таблицы, определяющей
характеристику
усилителя,
выводится
сигнал об ошибке.
- схема алгоритма
модели усилителя
Рисунок
5
- модель усилителя
SUBROUTINE
US(L,Z,U)
С
Подпрограмма
- модель усилителя.
DIMENSION
W(3,11)
C
характеристики
усилителя из
таблиц 2 и 3 по
столбцам
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. /
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
тестовой задачи
В
качестве тестовой
задачи вычислим
с малым шагом
и построим
графики характеристик
усилителя.
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
Рисунок
6
Из
рисунка видно,
что характеристика
усилителя
воспроизводится
моделью правильно.
вычисления
правых частей
системы уравнений
В
подпрограмме
сохранены
наименования
переменных
модели.
Результаты
вычислений
заносятся, как
это требуют
подпрограммы
шага и управляющего
модуля, в первый
столбец массива
F,
который здесь,
для простоты
объявлен одномерным.
Для передачи
в этот модуль
изменяемого
от эксперимента
к эксперименту
параметра
генератора
в общую область
включена переменная
TAU
.Остальные
переменные
общей области
нужны для связи
главного модуля
с подпрограммой
вывода результатов
шага.
SUBROUTINE
FUN(T,Z,F,N)
С
Подпрограмма
вычисления
правых частей
системы уравнений
модели автогенератора.
DIMENSION
Z(N*4),F(N*4),D(4,15000)
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
вывода
В
подпрограмме
сохранены
наименования
переменных
модели.
Для
того, чтобы
иметь возможность
хотя бы качественно,
но быстро, оценивать
правильность
работы модели
необходимо
осуществить
визуализацию
решения. Поэтому
в модуле вывода
на каждом шаге
вычислим входную
и выходную
переменные
усилителя и
заполним этими
данными очередной
столбец массива
D.
В этот же столбец
запишем текущие
значения времени
Т
. Массив D
передадим через
общую область
в главный модуль,
а оттуда подпрограмме
построения
графиков KRIS.
В
автогенераторе
некоторое время
длится процесс
самовозбуждения.
Нас интересует
процесс установившихся
колебаний,
поэтому запись
данных в массив
будем делать
только начиная
с некоторого
момента времени
TZ.
Эта
величина и
счетчик точек
также включим
в общую область.
SUBROUTINE
PRIN(T,Z,F,N,IER)
С
Подпрограмма
вывода результатов
шага.
DIMENSION
Z(N*4),F(N*4),D(4,15000)
COMMON
K,TZ,TAU,D
IF(T.GE.TZ)THEN
K=K+1
С
Вычисление
значения переменной
входа U1.
CALL
US(2,Z(3),U1)
C
Вычисление
значения переменной
выхода U2.
CALL
US(3,Z(3),U2)
С
Заполнение
массива.
D(1,K)=T
С
Выход усилителя
будет изображаться
на графиках
кривой номер
1.
D(2,K)=U2
С
Вход усилителя
будет изображаться
на графиках
кривой номер
2.
D(3,K)=U1
ENDIF
RETURN
END
модуль решения
системы уравнений
В
главном модуле
в соответствие
с требованиями
подпрограммы
метода Рунге
- Кутта ARK объявим
массивы для
решения системы
третьего порядка.
Имена массивов
сохраним такими
же, как имена
формальных
параметров
подпрограммы
ARK. Зададим
нулевые начальные
условия и равные
для всех интегральных
переменных
весовые коэффициенты
погрешности.
Из исходного
файла будем
вводить:
время
начала записи
данных в выходной
массив TZ ,
параметр
,
время
начала интегрирования
ТN,
время
конца интегрирования
ТК,
максимальный
шаг интегрирования
НМ
задаваемую
погрешность
ЕР.
DIMENSION
Z(12),RAB(9),F(12),D(4,15000)
С
Главный модуль
решения системы
уравнений
EXTERNAL
FUN,PRIN
COMMON
K,TZ,TAU,D
С
Задание начальных
условий и весовых
коэффициентов
погрешности.
DO
1 K=1,3
Z(K)=0.
1
Z(K+3)=0.33333
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
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
для спектрального
анализа выходного
напряжения
Выделим
один период
колебаний и
сделаем третье
решение.
Таблица
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
для установления
зависимостей
параметров
от
Изменяя
величину
,
делаем решения,
аналогичные
второму, и
результаты,
извлеченные
из выходных
файлов, заносим
в таблицу 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-ой
гармоники
можно вычислить,
как
(
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.405366*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
нужны в списке
ввода для считывания
данных прямо
из выходного
файла третьего
решения.
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%
от первой. В
таких условиях
аппроксимация
этой характеристики
не имеет смысла.
Б.П.
ДЕМИДОВИЧ,
И.А. МАРОН, Основы
вычислительной
математики,
«Наука», М., 1966.
Б.П.
ДЕМИДОВИЧ,
И.А. МАРОН, Э.З.
ШУВАЛОВА, Численные
методы анализа,
«Наука», М., 1967.
И.С.
БЕРЕЗИН, Н.П.
ЖИДКОВ, Методы
вычислений,
Физматгиз,
1961.
Н.Н.
КАЛИТКИН, Численные
методы, «Наука»,
М., 1978.
Н.С.
БАХВАЛОВ, Численные
методы, «Наука»,
М., 1975.
Д.
ХИММЕЛЬБЛАУ,
Прикладное
нелинейное
программирование,
«Мир», М., 1975.
А.А.
ФЕЛЬДБАУМ,
А.Д. ДУДЫКИН,
А.П. МАНОВЦЕВ,
Н.Н. МИРОЛЮБОВ,
Теоретические
основы связи
и управления,
Физматгиз, М.,
1963.
З.С.
БРИЧ, Д.В. КАПИЛЕВИЧ,
Н.А. КЛЕЦКОВА,
ФОРТРАН 77 для
ПЭВМ ЕС, «Финансы
и статистика»,
М., 1991.
П.В.
СОЛОВЬЕВ, FORTRAN
для
персонального
компьютера,
«ARIST»,
М., 1991.
Г.Н.
РЫБАЛЬЧЕНКО,
Численные
методы решения
задач строительства
на ЭВМ, Киев
УМК ВО, 1989.
Г.
Н. РЫБАЛЬЧЕНКО,
Методические
указания к
курсовой работе
по дисциплине
«Основы вычислительной
математики»,
Кривой Рог,
КТУ, 1997.
Исследование
RC-генератора
синусоидальных
колебаний.