Реферат: Исследование RC-генератора синусоидальных колебаний
--PAGE_BREAK--1. ПОСТАНОВКА ЗАДАЧИВыполнить исследование RC-генератора синусоидальных колебаний(Рисунок. 1)
<img width=«541» height=«170» src=«ref-1_1930778051-2621.coolpic» v:shapes="_x0000_i1027">
Рисунок 1
Генератор состоит из пассивной линейной части, включающей резисторы с сопротивлением Rи конденсаторы с емкостью С, и электронного усилителя с нелинейной характеристикой.
Передаточная функция линейной части
<img width=«412» height=«61» src=«ref-1_1930780672-1563.coolpic» v:shapes="_x0000_i1028">,
где <img width=«81» height=«21» src=«ref-1_1930782235-305.coolpic» v:shapes="_x0000_i1029">.
Нелинейная зависимость выходного напряжения <img width=«28» height=«28» src=«ref-1_1930782540-192.coolpic» v:shapes="_x0000_i1030">усилителя от его входного напряжения <img width=«27» height=«28» src=«ref-1_1930782732-192.coolpic» v:shapes="_x0000_i1031"> приведена в таблице 1
Таблица 1
Численными экспериментами на ЭВМ найти зависимости:
· периода Т установившихся автоколебаний от параметра <img width=«16» height=«17» src=«ref-1_1930777911-140.coolpic» v:shapes="_x0000_i1032"> ,
· амплитуды U2maxвыходного напряжения U2(t
)от <img width=«16» height=«17» src=«ref-1_1930777911-140.coolpic» v:shapes="_x0000_i1033"> ,
· амплитуды An
n-ойгармоники <img width=«193» height=«53» src=«ref-1_1930783204-766.coolpic» v:shapes="_x0000_i1034"> выходного напряжения от ее номера n <img width=«289» height=«63» src=«ref-1_1930783970-1228.coolpic» v:shapes="_x0000_i1035">,
· коэффициента усиления электронного усилителя в режиме установившихся автоколебаний от <img width=«16» height=«17» src=«ref-1_1930777911-140.coolpic» v:shapes="_x0000_i1036">.
Найденные экспериментально зависимости аппроксимировать степенными многочленами.
Из зависимости <img width=«79» height=«24» src=«ref-1_1930785338-320.coolpic» v:shapes="_x0000_i1037"> найти значение <img width=«16» height=«17» src=«ref-1_1930777911-140.coolpic» v:shapes="_x0000_i1038">, необходимое для получения периода автоколебаний <img width=«83» height=«20» src=«ref-1_1930785798-339.coolpic» v:shapes="_x0000_i1039">, и расчетом колебаний проверить правильность полученного значения параметра <img width=«16» height=«17» src=«ref-1_1930777911-140.coolpic» v:shapes="_x0000_i1040">.
Для вывода графиков и таблиц разрешается использовать библиотечную подпрограмму KRIS. Все остальные программные модули разработать самостоятельно.
2. МАТЕМАТИЧЕСКАЯ ФОРМУЛИРОВКА ЗАДАЧИ продолжение
--PAGE_BREAK--2.1 ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ ЛИНЕЙНОЙ ЧАСТИ
Запишем систему дифференциальных уравнений линейной части RC-генератора. Для этого преобразуем ее передаточную функцию
<img width=«405» height=«31» src=«ref-1_1930786277-1230.coolpic» v:shapes="_x0000_i1041"> ( 1 )
<img width=«488» height=«31» src=«ref-1_1930787507-1522.coolpic» v:shapes="_x0000_i1042"> ( 2 )
Введем первую вспомогательную переменную <img width=«51» height=«28» src=«ref-1_1930789029-313.coolpic» v:shapes="_x0000_i1043">, определяемую из уравнения
<img width=«172» height=«28» src=«ref-1_1930789342-604.coolpic» v:shapes="_x0000_i1044"> ( 3 )
Подставляя ( 3 ) в ( 2 ), получаем
<img width=«500» height=«31» src=«ref-1_1930789946-1594.coolpic» v:shapes="_x0000_i1045"> ( 4 )
Сокращая на <img width=«27» height=«21» src=«ref-1_1930791540-190.coolpic» v:shapes="_x0000_i1046"> и группируя в правой части члены, не содержащие <img width=«16» height=«20» src=«ref-1_1930791730-161.coolpic» v:shapes="_x0000_i1047">, получаем
<img width=«475» height=«31» src=«ref-1_1930791891-1520.coolpic» v:shapes="_x0000_i1048"> ( 5 )
Введем вторую вспомогательную переменную <img width=«52» height=«28» src=«ref-1_1930793411-310.coolpic» v:shapes="_x0000_i1049">, определяемую из уравнения
<img width=«248» height=«28» src=«ref-1_1930793721-830.coolpic» v:shapes="_x0000_i1050"> ( 6 )
Подставляя ( 6 ) в ( 5 ), получаем
<img width=«415» height=«31» src=«ref-1_1930794551-1317.coolpic» v:shapes="_x0000_i1051"> ( 7 )
Снова сокращая на <img width=«27» height=«21» src=«ref-1_1930791540-190.coolpic» v:shapes="_x0000_i1052"> и группируя в правой части члены, не содержащие <img width=«16» height=«20» src=«ref-1_1930791730-161.coolpic» v:shapes="_x0000_i1053">, получаем
<img width=«341» height=«28» src=«ref-1_1930796219-1072.coolpic» v:shapes="_x0000_i1054"> ( 8 )
Введем третью вспомогательную переменную <img width=«52» height=«28» src=«ref-1_1930797291-312.coolpic» v:shapes="_x0000_i1055">, определяемую из уравнения
<img width=«249» height=«28» src=«ref-1_1930797603-846.coolpic» v:shapes="_x0000_i1056"> ( 9 )
Подставляя ( 9 ) в ( 8 ) и сокращая на <img width=«27» height=«21» src=«ref-1_1930791540-190.coolpic» v:shapes="_x0000_i1057">, получаем
<img width=«201» height=«28» src=«ref-1_1930798639-723.coolpic» v:shapes="_x0000_i1058"> ( 10 )
Переходя в уравнениях ( 10 ), ( 9 ), ( 6 ), ( 3 ) от изображений переменных к их оригиналам, получаем систему уравнений
<img width=«113» height=«28» src=«ref-1_1930799362-355.coolpic» v:shapes="_x0000_i1059"> ( 11 )
<img width=«112» height=«53» src=«ref-1_1930799717-498.coolpic» v:shapes="_x0000_i1060"> ( 12 )
<img width=«176» height=«53» src=«ref-1_1930800215-741.coolpic» v:shapes="_x0000_i1061"> ( 13 )
<img width=«177» height=«53» src=«ref-1_1930800956-750.coolpic» v:shapes="_x0000_i1062"> ( 14 )
Здесь <img width=«28» height=«28» src=«ref-1_1930782540-192.coolpic» v:shapes="_x0000_i1063"> — функция, определяемая нелинейной характеристикой усилителя.
Так как генератор должен самовозбуждаться, то решение системы ( 11 ) — ( 14 ) можно выполнять от любых начальных условий, в том числе и от нулевых.
2.2 Уравнение усилителя
Уравнение ( 11 ) представляет собой нелинейное уравнение, которое необходимо решать при каждом вычислении правых частей системы.
Можно решать это уравнение методом итераций. Но есть более простой путь.
Найдем из характеристики усилителя разности <img width=«113» height=«28» src=«ref-1_1930801898-343.coolpic» v:shapes="_x0000_i1064">, а затем построим характеристику <img width=«105» height=«28» src=«ref-1_1930802241-415.coolpic» v:shapes="_x0000_i1065"> Значение <img width=«21» height=«28» src=«ref-1_1930802656-175.coolpic» v:shapes="_x0000_i1066"> известно сначала из начальных условий, а затем при каждом обращении к вычислению правых частей системы и из построенной нами характеристики всегда можно вычислить <img width=«27» height=«28» src=«ref-1_1930782732-192.coolpic» v:shapes="_x0000_i1067"> для подстановки в правые части остальных уравнений.
Вычисленная характеристика представлена в таблице 2.
Таблица 2
продолжение
--PAGE_BREAK--2.3 Конечно-элементная модель усилителя
Для построения квадратичного конечного элемента используем интерполяционную формулу Лагранжа
<img width=«539» height=«125» src=«ref-1_1930803023-3766.coolpic» v:shapes="_x0000_i1068"> ( 15 )
Для вычисления выходной величины автогенератора необходимо также по формуле Лагранжа по заданному значению <img width=«21» height=«28» src=«ref-1_1930802656-175.coolpic» v:shapes="_x0000_i1069"> находить <img width=«28» height=«28» src=«ref-1_1930782540-192.coolpic» v:shapes="_x0000_i1070">.
<img width=«552» height=«125» src=«ref-1_1930807156-3788.coolpic» v:shapes="_x0000_i1071"> ( 16 )
Данные в этом случае необходимо выбирать из таблицы 3, полученной из таблиц 1 и 2.
Таблица 3
3. ФУНКЦИОНАЛЬНОЕ ПРОЕКТИРОВАНИЕ ПРОГРАММНОГО КОМПЛЕКСА
Функционально программный комплекс должен состоять из двух независимых частей:
* программы — модели RC- генератора;
* набора программ обработки результатов моделирования автогенератора.
Модель RC- генератора должна, в свою очередь, включать:
* модуль, вызывающий подпрограмму метода Рунге — Кутта;
* модули метода Рунге — Кутта;
* модуль — модель усилителя;
* модуль правых частей ;
* модуль вывода результатов одного шага интегрирования.
Для программной реализации метода Рунге — Кутта удобно использовать два модуля:
* модуль, выполняющий один заданный шаг метода;
* модуль, управляющий величиной шага в зависимости от получаемой погрешности решения.
Взаимодействие этих модулей таково. Вызывающий модуль вводит значение параметра <img width=«16» height=«17» src=«ref-1_1930777911-140.coolpic» v:shapes="_x0000_i1072"> , начало и конец интервала интегрирования, максимальный шаг, начальные условия и заданную погрешность. Затем этот модуль обращается к модулю управления метода Рунге — Кутта. Последний задает величину шага подпрограмме одного шага и ведет процесс интегрирования системы уравнений, удерживая погрешность в заданных пределах. При выполнения шага, в соответствие с методом Рунге — Кутта, модуль шага четырежды обращается к модулю правых частей, а тот, в свою очередь, — к модели усилителя в виде функции <img width=«53» height=«28» src=«ref-1_1930811084-310.coolpic» v:shapes="_x0000_i1073">. После выполнения шага, удовлетворяющего условиям точности, модуль управления вызывает подпрограмму вывода результатов шага, а она, в свою очередь обращается к модели усилителя в виде функции <img width=«65» height=«28» src=«ref-1_1930811394-338.coolpic» v:shapes="_x0000_i1074">. Модуль управления заканчивает свою работу после достижения конца интервала интегрирования. Тогда вызывающий модуль обращается к подпрограмме вывода таблиц и графиков KRIS.
В набор подпрограмм обработки результатов моделирования необходимо включить две независимые программы:
* программу численного интегрирования по методу трапеций;
* программу аппроксимации экспериментальных зависимостей степенными многочленами методом наименьших квадратов.
продолжение
--PAGE_BREAK--4. МОДУЛИ РЕШЕНИЯ СИСТЕМЫ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ
4.1 Описание метода Рунге — Кутта четвертого порядка
Сначала рассмотрим применение метода для решения дифференциального уравнения, а затем для случая системы уравнений.
Пусть имеется уравнение
<img width=«107» height=«53» src=«ref-1_1930811732-569.coolpic» v:shapes="_x0000_i1075"> или <img width=«99» height=«25» src=«ref-1_1930812301-422.coolpic» v:shapes="_x0000_i1076">
где
· <img width=«15» height=«15» src=«ref-1_1930812723-164.coolpic» v:shapes="_x0000_i1077"> — неизвестная функция от независимой переменной <img width=«12» height=«19» src=«ref-1_1930812887-120.coolpic» v:shapes="_x0000_i1078">;
· <img width=«59» height=«24» src=«ref-1_1930813007-335.coolpic» v:shapes="_x0000_i1079"> — известная функция.
Все численные методы решения задачи Коши основаны на приближенной замене искомой функции степенными многочленами.
В методе Рунге-Кутта четвёртого порядка отыскиваетсяприращение, которое даёт приближающий многочлен на шаге интегрирования. Приращение искомой функции вычисляется ввиде произведения длины шага на значение производной от этойфункции. В качестве производной берется средневзвешенное от значений производных <img width=«155» height=«28» src=«ref-1_1930813342-462.coolpic» v:shapes="_x0000_i1080">
вычисленных в специально подобранных четырёх точках.
В качестве первой точки берут начальную точку шага<img width=«61» height=«28» src=«ref-1_1930813804-327.coolpic» v:shapes="_x0000_i1081">.
Производная в этой точке равна
<img width=«117» height=«28» src=«ref-1_1930814131-453.coolpic» v:shapes="_x0000_i1082"> ,
где <img width=«59» height=«24» src=«ref-1_1930813007-335.coolpic» v:shapes="_x0000_i1083"> — правая часть уравнения .
В качестве второй точки на плоскости решения <img width=«48» height=«24» src=«ref-1_1930814919-303.coolpic» v:shapes="_x0000_i1084"> выбираютточку с координатами<img width=«164» height=«53» src=«ref-1_1930815222-696.coolpic» v:shapes="_x0000_i1085">.
Производная во второйточке равна
<img width=«237» height=«53» src=«ref-1_1930815918-840.coolpic» v:shapes="_x0000_i1086">
Для третьей точки берут координаты <img width=«181» height=«53» src=«ref-1_1930816758-708.coolpic» v:shapes="_x0000_i1087"> ивычисляют производную
<img width=«239» height=«53» src=«ref-1_1930817466-836.coolpic» v:shapes="_x0000_i1088">
Наконец, для четвёртой точки берут координаты <img width=«185» height=«28» src=«ref-1_1930818302-544.coolpic» v:shapes="_x0000_i1089"> ивычисляют производную
<img width=«241» height=«28» src=«ref-1_1930818846-659.coolpic» v:shapes="_x0000_i1090">
По полученным четырём значениям производной находят средневзвешенное значение
<img width=«295» height=«53» src=«ref-1_1930819505-961.coolpic» v:shapes="_x0000_i1091">
Теперь, находят координаты конечной точкишага.
<img width=«149» height=«64» src=«ref-1_1930820466-683.coolpic» v:shapes="_x0000_i1092">
При решении системы уравнений вычисления ведут параллельно для каждого из уравнений.
4.2 Описание алгоритма одного шага
В алгоритме используются следующие идентификаторы
Таблица 4
Имя
Тип
Назначение
PRAV
Подпрограмма.
Подпрограмма, возвращающая значения производных.
N
Целый.
Порядок решаемой системы.
XN
Вещественный.
Исходный массив начальной точки шага.
XK
Вещественный.
Результирующий массив конечной точки шага.
F
Вещественный.
Массив возвращаемых подпрограммой РRAVпроизводных.
TN
Вещественный.
Начальное на шаге значение независимой переменной.
K
Целый.
Номер переменной.
J
Целый.
Номер частичного приращения.
T
Вещественный.
Независимая переменная.
H
Вещественный.
Задаваемая величина шага.
P
Вещественный
Массив размера (4,2), содержащий необходимые для вычисления и накопления приращений константы (0,.5,.5,1,6,3,3,6)
R
Вещественный
Рабочий массив размера (N,3)
Блок-схема алгоритма изображена на Рисунке 2. Номер переменной записан как верхний индекс.
В цикле с номерами блоков 2, 3, 4, 5 обнуляются второй и третий столбцы рабочего массива R
.
Внешний цикл с номерами блоков 6 — 18 вычисляет производные в четырех им формируемых точках и накапливает средневзвешенное значение приращений в третьем столбце рабочего массива R. Вдоль столбца расположены значения, соответствующие всем Nискомым переменным.
Блок 7 задает в цикле последовательно значения независимой переменной: TN, TN+0.5H, TN+0.5H, TN+H.Константы0, 0.5, 0.5 и 1 содержатся в первом столбце массива Р.
Цикл 8 — 11 записывает в первый столбец рабочего массива значения переменных для вычисления производных. Для этого к начальному значению переменной прибавляется сначала нулевое приращение, затем половина приращения, получаемого на шаге со значением производной в начальной точке, потом половина приращения, получаемого на шаге с значением производной во второй точке, и, наконец, полное приращение, получаемое на шаге со значением производной в третьей точке.
В блоке 12 выполняется обращение к подпрограмме вычисления производных. Подпрограмме передается значение независимой переменной и первый столбец рабочего массива, содержащий значения зависимых переменных в задаваемой точке. Подпрограмма возвращает массив Fзначений производных.
В цикле 13 — 16 во второй столбец рабочего массива по вычисленным значениям производных записываются приращения, а в третьем столбце накапливаются суммы четырех приращений с весовыми коэффициентами 1/6, 1/3, 1/3, 1/6 . Константы 6, 3, 3, 6 для этого хранятся во втором столбце массива Р.
В цикле 19 — 22 полученные приращения прибавляются к начальной точке и результат записывается в выходной массив.
В блоке 23 вычисляются производные в конечной точке шага.
продолжение
--PAGE_BREAK--4.3 Блок — схема алгоритма одного шага по методу Рунге — Кутта
<img width=«656» height=«888» src=«ref-1_1930821149-11331.coolpic» v:shapes="_x0000_i1093">
Рисунок 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
Вещественный.
Массив размера (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 выполняется подготовка к продолжению расчетов с прежним шагом. Из трех последних точек средняя делается начальной для выполнения контрольного шага удвоенной величины НВ, а последняя, — начальной для очередного малого шага Н.
продолжение
--PAGE_BREAK--4.6 Блок — схема алгоритма метода Рунге — Кутта с автоматическим выбором шага
<img width=«657» height=«902» src=«ref-1_1930832480-16890.coolpic» v:shapes="_x0000_i1094"> Рисунок 3
4.7 Подпрограмма метода Рунге — Кутта с автоматическим выбором шага
Подпрограмма ARK позволяет решать произвольную систему N-го порядка с автоматическим выбором шага интегрирования. Эта подпрограмма обращается:
· к подпрограмме одного шага- SH,
· к подпрограмме вычисления правых частей системы,
· к подпрограмме вывода.
Подпрограмма SH записана в универсальном виде и приведена выше. Головной вызывающий модуль, а также подпрограммы правых частей и вывода Пользователь должен составить самостоятельно.
В главном модуле оператором EXTERNAL должны быть объявлены имена подпрограмм правых частей и вывода. ОператоромDIMENSION должны быть объявлены массивы — фактические параметры подпрограммы ARK. Эти массивы, по желанию, могут объявляться как одномерные или как двухмерные. Размеры массивов (N,4),(N,3),(N,4), где N-порядок системы. Формальные имена этих массивов в подпрограмме A
RK, соответственно, 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,6
4 IER=1
GO TO 6
5 IER=2
6 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 14K=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,3
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 38K=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 35K=1,N
35 X(K,1)=X(K,3)
KLP=1
GO TO 6
END
продолжение
--PAGE_BREAK--4.8 Тестовая задача
Решим дифференциальное уравнение
<img width=«84» height=«56» src=«ref-1_1930849370-437.coolpic» v:shapes="_x0000_i1095">
с начальными условиями <img width=«179» height=«27» src=«ref-1_1930849807-660.coolpic» v:shapes="_x0000_i1096"> . Легко видеть, что решением этой задачи является функция <img width=«91» height=«24» src=«ref-1_1930850467-433.coolpic» v:shapes="_x0000_i1097"> . Вычислим решение на интервале <img width=«100» height=«20» src=«ref-1_1930850900-444.coolpic» v:shapes="_x0000_i1098"> , что составит почти <img width=«25» height=«20» src=«ref-1_1930851344-203.coolpic» v:shapes="_x0000_i1099"> периода этой функции. Если при таком длительном интегрировании амплитуда косинусоиды существенно не изменится, то алгоритм численно устойчив. Можно также сравнить решение в конечной точке <img width=«247» height=«24» src=«ref-1_1930851547-1017.coolpic» v:shapes="_x0000_i1100">
Подпрограмма правых частей для этого уравнения будет такой.
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.9 Результаты тестирования
Графики вычисленных путем решения дифференциального уравнения функций приведены на рисунке 4. Видно, что они близки к функциям <img width=«91» height=«24» src=«ref-1_1930850467-433.coolpic» v:shapes="_x0000_i1101"> и <img width=«111» height=«25» src=«ref-1_1930852997-430.coolpic» v:shapes="_x0000_i1102">.
<img width=«614» height=«336» src=«ref-1_1930853427-13199.coolpic» v:shapes="_x0000_i1103">
Рисунок 4
Амплитуды колебаний равны единице, период <img width=«80» height=«24» src=«ref-1_1930866626-404.coolpic» v:shapes="_x0000_i1104"> .
Выходной файл решения приведен ниже.
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% . При этом следует учесть, что в эту погрешность вошла и погрешность процедуры нахождения максимума с шагом, равным шагу интегрирования. Тенденции к затуханию или раскачиванию колебаний нет. Все это доказывает работоспособность алгоритма и программы.
продолжение
--PAGE_BREAK--4.10 Квадратичная конечно-элементная модель усилителя 4.10.1 Описание алгоритма
Функция этого модуля заключается в том, что по заданной входной величине (обозначим ее Z3) выдается или величина U1, определяемая из таблицы 2, или величина U2, определяемая из таблицы 3. Эти таблицы представим в виде одного двухмерного массива W, в первой строке которого запишем табличные значения входной переменной Z3, а во второй и третьей строках — им соответствующие табличные значения переменных U1и U2 . Значение еще одного входного параметра L ,— номера строки, будет определять, какую выходную переменную вычисляет модель (L=2 или L=3). Выходную переменную модуля обозначим U, а для модуля назначим имя US. Блок — схема алгоритма приведена на рисунке 5.
В цикле с индексом J определяется тот конечный элемент, в области которого находится входная величина Z3, а затем вычисляется выходная величина по формуле Лагранжа с использованием L-той строки массива W.
<img width=«644» height=«125» src=«ref-1_1930867030-4962.coolpic» v:shapes="_x0000_i1105">
Если значение входной переменной Z3выходит за пределы таблицы, определяющей характеристику усилителя, выводится сигнал об ошибке.
4.10.2 Блок — схема алгоритма модели усилителя
<img width=«596» height=«462» src=«ref-1_1930871992-4439.coolpic» v:shapes="_x0000_i1106">
Рисунок 5
4.10.3 Подпрограмма — модель усилителя
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
4.10.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
<img width=«614» height=«336» src=«ref-1_1930876431-5810.coolpic» v:shapes="_x0000_i1107">
Рисунок 6
Из рисунка видно, что характеристика усилителя воспроизводится моделью правильно.
продолжение
--PAGE_BREAK--4.11 Подпрограмма вычисления правых частей системы уравнений
В подпрограмме сохранены наименования переменных модели.Результаты вычислений заносятся, как это требуют подпрограммы шага и управляющего модуля, в первый столбец массива F,который здесь, для простоты объявлен одномерным. Для передачи в этот модуль изменяемого от эксперимента к эксперименту параметра генератора <img width=«16» height=«17» src=«ref-1_1930777911-140.coolpic» v:shapes="_x0000_i1108"> в общую область включена переменная 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
4.12 Подпрограмма вывода
В подпрограмме сохранены наименования переменных модели.
Для того, чтобы иметь возможность хотя бы качественно, но быстро, оценивать правильность работы модели необходимо осуществить визуализацию решения. Поэтому в модуле вывода на каждом шаге вычислим входную и выходную переменные усилителя и заполним этими данными очередной столбец массива 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
4.13 Главный модуль решения системы уравнений
В главном модуле в соответствие с требованиями подпрограммы метода Рунге — Кутта ARK объявим массивы для решения системы третьего порядка. Имена массивов сохраним такими же, как имена формальных параметров подпрограммы ARK. Зададим нулевые начальные условия и равные для всех интегральных переменных весовые коэффициенты погрешности. Из исходного файла будем вводить:
* время начала записи данных в выходной массив TZ ,
* параметр <img width=«16» height=«17» src=«ref-1_1930777911-140.coolpic» v:shapes="_x0000_i1109"> ,
* время начала интегрирования Т
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
продолжение
--PAGE_BREAK--5. ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ АВТОГЕНЕРАТОРА 5.1 Пробные решения
Пробное решение выполним с параметрами, указанными в таблице 6
Таблица 6
<img width=«614» height=«336» src=«ref-1_1930882661-15407.coolpic» v:shapes="_x0000_i1111">
Рисунок 7
Из рисунка видно, что возбуждение автогенератора длится примерно 20 периодов колебаний, период колебания примерно равен 16с., что составляет <img width=«36» height=«21» src=«ref-1_1930898068-233.coolpic» v:shapes="_x0000_i1112">.
Второе решение выполним так, чтобы запись началась в режиме установившихся колебаний и длилась около двух периодов. Тогда по таблице решения можно с достаточной точностью установить амплитуду и период колебаний. Данные для второго решения приведены в таблице 7.
Таблица 7
Графики решения приведены на Рисунке 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
<img width=«614» height=«336» src=«ref-1_1930898441-9816.coolpic» v:shapes="_x0000_i1114">
Рисунок 8
продолжение
--PAGE_BREAK--
еще рефераты
Еще работы по математике
Реферат по математике
Группы преобразований
3 Сентября 2013
Реферат по математике
Шпоры по математическому анализу
3 Сентября 2013
Реферат по математике
Свойства пространства с некоторыми компактифицированными измерениями
3 Сентября 2013
Реферат по математике
Исследование взаимосвязей показателей финансовых рынков
20 Июня 2015