Южный федеральный университет
ОЦЕНКА ТОЧНОСТИ ЧИСЛЕННОГО АНАЛИЗА РЕЛАКСАЦИОННОГО ГЕНЕРАТОРА
А.М. Пилипенко
В.Н. Бирюков
Компьютерное моделирование радиотехнических цепей и устройств во временной области заключается в численном решении систем обыкновенных дифференциальных уравнений (СОДУ). В современных пакетах схемотехнического проектирования и автоматизации инженерных расчетов реализовано большое количество численных методов решения СОДУ. Основными критериями эффективности численных методов являются их вычислительная сложность и точность получаемых решений. Точность компьютерного моделирования количественно характеризует степень отклонения результатов моделирования от истинных результатов. Для теоретической оценки точности численных методов решения СОДУ используют тестовые, как правило, линейные СОДУ, для которых известно точное решение. Реальные радиотехнические устройства описываются нелинейными СОДУ, поэтому оценка точности их численного анализа представляет собой достаточно сложную задачу.
В работе [1] было проведено исследование точности современных численных методов при анализе автогенераторов гармонических колебаний. Данная статья является продолжением работы [1] применительно к другому классу автоколебательных цепей - генераторам релаксационных колебаний. Целью настоящей работы является оценка точности современных численных методов решения СОДУ при моделировании релаксационных генераторов.
Выбор задачи анализа генератора релаксационных колебаний в качестве тестовой обусловлен следующими причинами. Во-первых, релаксационный генератор является неотъемлемой частью многих радиоэлектронных устройств, таких как функциональные генераторы, радиоизмерительные приборы, счетчики, таймеры, устройства, инициирующие измерения или технологические процессы [2]. Во-вторых, задача численного анализа релаксационного генератора во временной области является одной из наиболее сложных для компьютерного моделирования, поскольку СОДУ, описывающая автогенератор, часто оказывается одновременно колебательной и жесткой [3]. Таким образом, оценка точности численного анализа релаксационного генератора необходима для выбора наиболее эффективного метода моделирования автоколебательных систем с различной степенью жесткости.
1. Модель релаксационного генератора
Релаксационные колебания возникают в нелинейных электрических цепях, не обладающих частотной избирательностью. Наиболее простая схемная модель релаксационного генератора может быть представлена также как и модель генератора гармонических колебаний - в виде параллельного соединения трех элементов - линейных емкости C и индуктивности L и нелинейного резистивного элемента, дифференциальное сопротивление (проводимость) которого может изменять знак [1].
Математическая модель релаксационного генератора (уравнение автогенератора) имеет вид нелинейного дифференциального уравнения второго порядка:
где u - напряжение на выходе генератора, - резонансная частота LC-контура, - дифференциальная проводимость нелинейного элемента, i(u) - вольт-амперная характеристика (ВАХ) нелинейного элемента.
Если дифференциальная проводимость G(u) на начальном участке ВАХ вблизи нуля принимает отрицательное значение, то в LC-контуре возникают автоколебания. При аппроксимации ВАХ нелинейного элемента полиномом третьей степени i(u) = a1 u + a3 u3 (a1 < 0, a3 > 0) уравнение (1) называется уравнением Ван-дер-Поля и для установившегося режима известны его аналитические решения. К сожалению, решение в явном виде получить не удается, и существуют только эталонные численные решения, специально вычисленные с высокой точностью для некоторых параметров уравнения и начальных данных [4].
Для получения аналитического решения уравнения автогенератора в установившемся режиме можно использовать методику, предложенную в работе [5]. Эта методика заключается в применении кусочно-линейной аппроксимации ВАХ нелинейного элемента. Такой подход позволяет определить стационарное решение уравнения автогенератора с точностью до последнего знака разрядной сетки компьютера.
Для релаксационного генератора характерны значительные потери, которые можно оценить затуханием , где ? характеристическая проводимость LC-контура, G1- дифференциальная проводимость нелинейного элемента при di / du > 0. В случае значительных потерь d > 2 и решение линейного однородного дифференциального уравнения второго порядка описывается суммой двух экспоненциальных функций. Таким образом, при использовании симметричной кусочно-линейной функции для аппроксимации ВАХ нелинейного элемента аналитическое выражение для напряжения на выходе релаксационного генератора в установившемся режиме имеет вид [1]
где A1, A2, A3, A4 - постоянные интегрирования; t0 - момент времени, в который проводимость нелинейного элемента меняет знак; T - период колебаний; и ? постоянные времени автогенератора при G(u) = G0 < 0 и G(u) = G1 > 0 соответственно; .
На рис. 1 показаны точные стационарные решения уравнения автогенератора, полученные с помощью выражения (2) при двух различных значениях затухания d = 4 (сплошная линия) и d = 100 (штриховая линия), значения остальных параметров автогенератора были заданы следующим образом: L = 1 Гн, С = 1 Ф, G0 = ? 4 См.
Рис. 1. Точные стационарные решения уравнения автогенератора.
Следует отметить, что при d = 4 модель автогенератора (1) может считаться слабо жесткой, так как в этом случае коэффициент жесткости, равный отношению максимальной и минимальной постоянных времени автогенератора з = фmax / фmin ? 10 > 1. При d = 100 модель (1) является жесткой, так как в этом случае з ? 104 >> 1.
Полученные аналитические решения позволяют с точностью до последнего знака разрядной сетки компьютера оценить основные параметры колебаний на выходе автогенератора. В таблице 1 приведены точные максимальные значения Am0 и точные значения частот fm0 релаксационных колебаний, представленных на рис. 1.
Таблица 1.
|
d |
з |
Am0, В |
fm0, Гц |
|
|
4 |
10 |
2.765625375929214 |
0.09132317033942401 |
|
|
100 |
10 4 |
1.084334602144661 |
0.05506258511132066 |
2. Численное решение уравнения автогенератора
В данной работе проведена оценка точности моделирования релаксационного генератора во временной области при использовании трех численных методов решения СОДУ: метода трапеций (TR), Гира (BDF) и RADAU5. Подробное описание указанных методов приведено в работах [6, 7]. Выбор для численного анализа методов TR и BDF объясняется тем, что они являются основными методами для анализа переходных процессов в программах схемотехнического моделирования. Метод RADAU5 является одним из наиболее перспективных методов типа Рунге-Кутты для решения жестких задач [7]. Следует отметить, что метод RADAU5 является трехстадийным, поэтому его применение сводится к решению системы алгебраических уравнений, порядок которой в три раза больше порядка СОДУ заданной цепи. Методы TR и BDF - одностадийные, таким образом, вычислительная сложность этих методов, определяющая затраты машинного времени и оперативной памяти на получение решения, минимальна. Многостадийные методы обладают повышенной L-устойчивостью и рекомендуются для решения задач высокой жесткости [8, 9].
Для оценки точности перечисленных выше численных методов применим методику, предложенную в работе [10] и основанную на анализе точности оценки основных параметров генерируемого колебания - частоты и амплитуды (максимального значения). Текущие относительные погрешности оценки частоты и амплитуды соответственно определяются с помощью следующих соотношений
и ,
где fm(t) и Am(t) - оценки частоты и амплитуды генерируемого колебания, полученные из решения уравнения (1) численными методами.
На рис. 2 и рис. 3 приведены временные диаграммы относительных погрешностей оценки частоты и амплитуды стационарного колебания релаксационного генератора при различном затухании и соответственно жесткости колебательной системы (см. таблицу 1). Решение системы уравнений (1) осуществлялось перечисленными выше численными методами при заданном максимальном шаге hmax = T / 1000 и предельно допустимой относительной погрешности TOL = 10 ? 5. Для оценки частоты колебаний использовалась линейная интерполяция численного решения в точках изменения знака u(t), для оценки амплитуды - квадратичная в точках локальных экстремумов. Интервал наблюдения выбран соответствующим стационарному режиму генерации.
Из рис. 2 и рис. 3 следует, что параметры процесса fm(t) и Am(t) изменяются от периода к периоду существенно. Для количественного описания полученных результатов в таблице 2 приведены средние значения относительных погрешностей оценки частоты и амплитуды колебаний (mещ и mеA соответственно) и среднеквадратические отклонения (СКО) относительных погрешностей оценки частоты и амплитуды колебаний (уещ и уеA соответственно).
Рис. 2. Текущие относительные погрешности оценки частоты релаксационного колебания методами TR (черная кривая), BDF (синяя кривая) и RADAU5 (красная кривая) приhmax = T / 1000, TOL = 10 ? 5 и различной жесткости модели автогенератора з = 10 (а) и з = 10 4 (б).
Рис. 3. Текущие относительные погрешности оценки амплитуды релаксационного колебания методами TR (черная кривая), BDF (синяя кривая) и RADAU5 (красная кривая) приhmax = T / 1000, TOL = 10 ? 5 и различной жесткости модели автогенератора з = 10 (а) и з = 10 4 (б).
Таблица 2.
|
з |
10 |
10 4 |
|||||
|
Метод |
TR |
BDF |
RADAU5 |
TR |
BDF |
RADAU5 |
|
|
mещ |
2.7•10?5 |
3.1•10?5 |
1.2•10?4 |
8.7•10?5 |
5.5•10?5 |
1.7•10?3 |
|
|
уещ |
4.8•10?6 |
4.5•10?6 |
5.1•10?4 |
1.7•10?6 |
7.2•10?6 |
5.8•10?3 |
|
|
mеA |
5.6•10?6 |
1.4•10?5 |
4.7•10?5 |
4.0•10?4 |
2.6•10?5 |
2.2•10?4 |
|
|
уеA |
3.7•10?6 |
4.2•10?6 |
2.1•10?4 |
1.3•10?4 |
1.9•10?5 |
2.3•10?3 |
Приведенные на рис. 2, рис. 3 и в таблице 2 результаты численного анализа показывают, что метод RADAU5, как правило, имеет наибольшие значения относительных погрешностей как при оценке частоты, так и при оценке амплитуды релаксационных колебаний независимо от жесткости задачи. Кроме того, текущие значения погрешностей для метода RADAU5 могут на один-два порядка превышать аналогичные погрешности методов BDF и TR. Особенно ярко этот эффект проявляется при оценке частоты в случае высокой жесткости задачи (см. рис. 2, б). В случае слабой жесткости системы (1) (з = 10) относительные погрешности для методов BDF и TR имеют одинаковый порядок величины. При высокой жесткости (з = 10 4) текущие погрешности еf (t) и еA(t) оказываются наименьшими для метода BDF, причем среднее значение погрешности оценки амплитуды для метода BDF меньше аналогичного значения для метода TR примерно в 15 раз. Следует также отметить, что погрешности метода BDF наиболее слабо зависят от жесткости задачи по сравнению с погрешностями других рассмотренных методов. Как видно из таблицы 2 при увеличении жесткости задачи в 1000 раз относительные погрешности для метода BDF возрастают не более чем в 4.5 раза, а для методов TR и RADAU5 могут увеличиваться более чем в 50 раз.
Из рис. 2 и рис. 3 следует, что для методов TR и BDF может наблюдаться эффект «синхронизации» ошибки с периодом генерируемых колебаний. При этом текущие погрешности имеют, как правило, периодический характер, причем период изменения погрешностей кратен периоду генерируемых колебаний. Кроме «ошибки синхронизации» второй случайной составляющей является ошибка, наблюдаемая в методах BDF и RADAU5, вызванная алгоритмом изменения шага в процессе решения.
Для проверки корректности реализованной методики оценки точности на рис. 4 и рис. 5 приведены текущие относительные погрешности оценки частоты и амплитуды релаксационного колебания методами TR, BDF и RADAU5 при уменьшении hmax в 2 раза и TOL в 4 раза по сравнению со случаем представленным на рис. 2 и рис. 3.
Из сравнительного анализа погрешностей, приведенных на рис. 2, 3 и на рис. 4, 5 видно, что при уменьшении TOL в 4 раза средние значения и СКО погрешностей методов TR и BDF также уменьшаются примерно в 4 раза, а эффект синхронизации ослабевает. Для метода RADAU5 уменьшение погрешностей происходит только в случае слабо жесткой задачи, а в случае жесткой задачи наблюдается значительный рост текущих погрешностей еf (t) и еA(t). Такой эффект для метода RADAU5 можно объяснить резким ростом погрешности округления при уменьшении шага решения. Погрешность округления метода RADAU5 возрастает из-за плохой обусловленности соответствующей ему системы алгебраических уравнений, порядок которой в три раза больше порядка исходной СОДУ [7].