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

СОДЕРЖАНИЕ

TOC \o “1-3” 1. ПОСТАНОВКА ЗАДАЧИ GOTOBUTTON _Toc439340136
PAGEREF _Toc439340136 3

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4.10.2 Блок – схема алгоритма модели усилителя GOTOBUTTON
_Toc439340154 PAGEREF _Toc439340154 18

4.10.3 Подпрограмма – модель усилителя GOTOBUTTON _Toc439340155
PAGEREF _Toc439340155 18

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

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

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

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

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

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

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

GOTOBUTTON _Toc439340163 PAGEREF _Toc439340163 25

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

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

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

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

7. ЛИТЕРАТУРА GOTOBUTTON _Toc439340168 PAGEREF _Toc439340168 29

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

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

Рисунок SEQ Рисунок \* ARABIC 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

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

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

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

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

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

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

( SEQ НФ\n \*
MERGEFORMAT 1 )

( SEQ НФ\n \* MERGEFORMAT 2 )

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

( SEQ НФ\n \* MERGEFORMAT 3 )

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

( SEQ НФ\n \* MERGEFORMAT 4 )

, получаем

( SEQ НФ\n \* MERGEFORMAT 5 )

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

( SEQ НФ\n \* MERGEFORMAT 6 )

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

( SEQ НФ\n \* MERGEFORMAT
7 )

, получаем

( SEQ
НФ\n \* MERGEFORMAT 8 )

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

( SEQ НФ\n \* MERGEFORMAT 9 )

, получаем

( SEQ НФ\n \* MERGEFORMAT 10 )

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

( SEQ НФ\n \*
MERGEFORMAT 11 )

( SEQ НФ\n \*
MERGEFORMAT 12 )

( SEQ НФ\n \* MERGEFORMAT 13 )

( SEQ НФ\n \* MERGEFORMAT 14 )

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

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

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

Уравнение ( 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

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

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

( SEQ НФ\n \* MERGEFORMAT 15 )

( SEQ НФ\n \* MERGEFORMAT 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

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

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

программы – модели RC – генератора;

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

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

модуль, вызывающий подпрограмму метода Рунге – Кутта;

модули метода Рунге – Кутта;

модуль – модель усилителя;

модуль правых частей ;

модуль вывода результатов одного шага интегрирования.

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

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

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

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

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

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

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

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

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

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

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

где

;

– известная функция.

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

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

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

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

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

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

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

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

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

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

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

Таблица 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 вычисляются производные в конечной точке шага.

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

Рисунок SEQ Рисунок \* ARABIC 2

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

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

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

Блок -схема алгоритма приведена на Рисунке 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
выполняется подготовка к продолжению расчетов с прежним шагом. Из трех
последних точек средняя делается начальной для выполнения контрольного
шага удвоенной величины НВ, а последняя , – начальной для очередного
малого шага Н.

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

Рисунок SEQ Рисунок \* ARABIC 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

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

Рисунок SEQ Рисунок \* ARABIC 4

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Рисунок SEQ Рисунок \* ARABIC 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

Рисунок SEQ Рисунок \* ARABIC 6

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

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

в общую область включена переменная 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

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

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

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

Таблица 6

TN TK HM EP

0 1 0 370 1 0.0001

Рисунок SEQ Рисунок \* ARABIC 7

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

Таблица 7

TN TK HM EP

370 1 0 400 1 0.0001

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

Таблица 8

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

370.0 -1.753 .5084E-01 .0000

370.5 -1.291 .3469E-01 .0000

371.0 -.7804 .1970E-01 .0000

371.5 -.2281 .6177E-02 .0000

372.0 .3466 -.8225E-02 .0000

372.5 .9243 -.2303E-01 .0000

373.0 1.476 -.4105E-01 .0000

373.5 1.974 -.5888E-01 .0000

374.0 2.395 -.7481E-01 .0000

374.5 2.699 -.9564E-01 .0000

375.0 2.860 -.1103 .0000

375.5 2.885 -.1127 .0000

376.0 2.792 -.1037 .0000

376.5 2.600 -.8794E-01 .0000

377.0 2.324 -.7205E-01 .0000

377.5 1.961 -.5838E-01 .0000

378.0 1.527 -.4280E-01 .0000

378.5 1.038 -.2625E-01 .0000

379.0 .5052 -.1226E-01 .0000

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

380.0 -.6338 .1614E-01 .0000

380.5 -1.202 .3169E-01 .0000

381.0 -1.729 .4996E-01 .0000

381.5 -2.190 .6695E-01 .0000

382.0 -2.559 .8495E-01 .0000

382.5 -2.793 .1038 .0000

383.0 -2.885 .1127 .0000

383.5 -2.849 .1092 .0000

384.0 -2.706 .9619E-01 .0000

384.5 -2.472 .7926E-01 .0000

385.0 -2.152 .6553E-01 .0000

385.5 -1.753 .5082E-01 .0000

386.0 -1.290 .3467E-01 .0000

386.5 -.7795 .1968E-01 .0000

387.0 -.2272 .6154E-02 .0000

387.5 .3476 -.8250E-02 .0000

388.0 .9253 -.2306E-01 .0000

388.5 1.477 -.4108E-01 .0000

389.0 1.975 -.5892E-01 .0000

389.5 2.396 -.7484E-01 .0000

390.0 2.699 -.9568E-01 .0000

390.5 2.861 -.1103 .0000

391.0 2.885 -.1127 .0000

391.5 2.791 -.1037 .0000

392.0 2.600 -.8792E-01 .0000

392.5 2.323 -.7203E-01 .0000

393.0 1.960 -.5836E-01 .0000

393.5 1.526 -.4277E-01 .0000

394.0 1.037 -.2622E-01 .0000

394.5 .5042 -.1223E-01 .0000

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

395.5 -.6350 .1617E-01 .0000

396.0 -1.203 .3172E-01 .0000

396.5 -1.730 .4999E-01 .0000

397.0 -2.191 .6699E-01 .0000

397.5 -2.560 .8500E-01 .0000

398.0 -2.793 .1039 .0000

398.5 -2.885 .1127 .0000

399.0 -2.849 .1091 .0000

399.5 -2.705 .9616E-01 .0000

400.0 -2.472 .7922E-01 .0000

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

Рисунок SEQ Рисунок \* ARABIC 8

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

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

Таблица 9

TN TK HM EP

379,5 1 0 395 1 0.0001

Рисунок SEQ Рисунок \* ARABIC 9

Таблица 9

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

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

380.0 -.6338 .1614E-01 .0000

380.5 -1.202 .3169E-01 .0000

381.0 -1.729 .4996E-01 .0000

381.5 -2.190 .6695E-01 .0000

382.0 -2.559 .8495E-01 .0000

382.5 -2.793 .1038 .0000

383.0 -2.885 .1127 .0000

383.5 -2.849 .1092 .0000

384.0 -2.706 .9619E-01 .0000

384.5 -2.472 .7926E-01 .0000

385.0 -2.152 .6553E-01 .0000

385.5 -1.753 .5082E-01 .0000

386.0 -1.290 .3467E-01 .0000

386.5 -.7795 .1968E-01 .0000

387.0 -.2272 .6154E-02 .0000

387.5 .3476 -.8250E-02 .0000

388.0 .9253 -.2306E-01 .0000

388.5 1.477 -.4108E-01 .0000

389.0 1.975 -.5892E-01 .0000

389.5 2.396 -.7484E-01 .0000

390.0 2.699 -.9568E-01 .0000

390.5 2.861 -.1103 .0000

391.0 2.885 -.1127 .0000

391.5 2.791 -.1037 .0000

392.0 2.600 -.8792E-01 .0000

392.5 2.323 -.7203E-01 .0000

393.0 1.960 -.5836E-01 .0000

393.5 1.526 -.4277E-01 .0000

394.0 1.037 -.2622E-01 .0000

394.5 .5042 -.1223E-01 .0000

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

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

Таблица 10

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

( SEQ НФ\n \* MERGEFORMAT
17 )

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

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

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

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

(
SEQ НФ\n \* MERGEFORMAT 18 )

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

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

( SEQ НФ\n \*
MERGEFORMAT 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

†††††??????††††††?????†††††???????†††††????†††††??????????????††††††????
????????††???????????††††††??

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

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

Рисунок SEQ Рисунок \* ARABIC 10

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

Зависимость амплитуды гармоники от ее номера приведены в таблицах 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

Рисунок SEQ Рисунок \* ARABIC 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.

PAGE

PAGE 6

IEIENOA?NOAI IA?ACIAAIEss OE?AEIU

E?eai?iaeneee oaoie/aneee oieaa?neoao

Eaoaae?a iiaeaee?iaaiey e i?ia?aiiiiai iaania/aiey

EO?NIAAss ?AAIOA

ii aeenoeeieeia

« Iniiau au/eneeoaeueiie iaoaiaoeee »

Auiieiee no. a?oiiu
ICIO-96

I?EOIAe*AIEI A.

?oeiaiaeeoaeue
no. i?aiiaeaaaoaeue

?UAAEUe*AIEI A.I.

E?eaie ?ia

1998

IA*AEI

R2K=0

R3K=0

AeA

IAO

R1K=XNK+P1J. R2K

AeA

IAO

AeA

IAO

AeA

IAO

K=1

K=K+1

|2HB|

IA*AEI

US(L,Z,U)

W-IANNEA (3,11)

J=2

AeA

IAO

Z3>W(1,J-1) Z3

Нашли опечатку? Выделите и нажмите CTRL+Enter