Жесткая задача с экспоненциальной нелинейностью для тестирования программ математического и схемотехнического моделирования
А.М. Пилипенко Южный федеральный университет, Таганрог
Аннотация: Представлена оригинальная тестовая задача для оценки эффективности методов численного решения обыкновенных дифференциальных уравнений и методов численного анализа переходных процессов в электронных цепях и системах. Тестовая задача описана как в виде математической модели, так и в виде модели электрической цепи. Предложенные в работе тестовая задача и методика оценки точности численного решения могут использоваться в программах математического и схемотехнического моделирования для тестирования как известных, так и новых методов численного моделирования электронных цепей и систем с различной динамикой, в том числе описывающихся жесткими и/или нелинейными дифференциальными уравнениями. Ключевые слова: численные методы, электронные цепи, анализ переходных процессов, обыкновенные дифференциальные уравнения, жесткие системы.
Для разработки высокотехнологичной радиоэлектронной аппаратуры в настоящее время важнейшее значение приобретают современные системы автоматизированного проектирования (САПР), позволяющие обеспечить высокую точность моделирования различных электронных устройств [1]. В связи с постоянным совершенствованием методов моделирования электронных устройств актуальной задачей является разработка тестов, позволяющих оценить эффективность указанных методов [2].
Моделирование динамики различных электронных устройств, например, радиотехнических цепей, как правило, состоит в решении системы обыкновенных дифференциальных уравнений (ОДУ), описывающей рассматриваемое устройство во временной области. В качестве среды для анализа радиотехнических цепей во временной области можно использовать как программы математического моделирования (Mathcad, MATLAB), так и программы схемотехнического моделирования (SPICE, Microcap, NI Multisim). В указанные программы включено большое количество методов предназначенных для численного решения различных типов ОДУ [3, 4].
Задачи моделирования сложных радиотехнических цепей относятся к жестким задачам, которые достаточно трудны для численного решения [5]. Решение жесткой задачи (жесткой системы ОДУ) содержит как быстрые, так и медленные компоненты, скорости изменения которых отличаются на порядок и более. Жесткость системы можно оценить отношением максимальной и минимальной постоянных времени системы.
В настоящее время для тестирования программ математического моделирования широко применяются уравнение Далквиста и уравнение Ван дер Поля [3 - 5]. Указанные задачи так же можно использовать и для тестирования программ схемотехнического моделирования, так как уравнению Далквиста соответствует линейная ЛС-цепь, а уравнению Ван дер Поля - нелинейный осциллятор. Основной недостаток задачи Далквиста заключается, в том, что ее решение дает только приближенные сведения о погрешности численных методов на комплексной плоскости [6], кроме того уравнение Далквиста линейное и не является жестким, т.е. не позволяет учесть основные трудности, возникающие при анализе радиотехнических цепей. Уравнение Ван дер Поля является нелинейным и позволяет моделировать как жесткие, так и колебательные цепи, но данное уравнение, как и ряд других тестов, основанных на прикладных задачах, имеет существенный недостаток - не известно явное выражение для его точного решения [7].
Цель настоящей статьи - разработать тестовую задачу в виде ОДУ и схемной модели с произвольной жесткостью и нелинейностью для исследования эффективности методов численного анализа переходных процессов в программах математического и схемотехнического моделирования. Кроме нелинейности, жесткости и возможности схемного описания, важнейшим требованием, предъявляемым к разрабатываемой тестовой задаче, является наличие точного решения в явном виде.
Наиболее распространенная в математике задача для тестирования методов численного интегрирования ОДУ представляет собой линейное дифференциальное уравнение первого порядка, также называемое уравнением Далквиста [8]:
— = Хи, и(0) = и0
(1)
Решение уравнения (1) описывается следующим выражением
и(г) = и0 ехр(Хг). (2)
Нетрудно показать, что уравнение Далквиста описывает линейную ЯС-цепь, схема которой приведена на рис. 1, а. Постоянная времени ЯС-цепи тС и параметр X, входящий в уравнение (1), связаны следующим соотношением тС = ЯС = -1 / X, где Я и С - параметры элементов ЯС-цепи.
Рис. 1. - Известная (а) и предлагаемая (б) схемы тестовых цепей
Для получения нелинейной тестовой цепи (рис. 1, б) заменим в схеме ЯС-цепи линейное сопротивление Я нелинейным элементом Будем полагать, что в качестве нелинейного элемента используется диод, вольт-амперная характеристика (ВАХ) которого описывается моделью Шокли [9]:
(ио ) = 13 [еХР(и£ / ф) - ^ (3)
где и0 - напряжение на диоде; 1з и ф - основные параметры диода.
Кроме того, будем полагать, что напряжение источника имеет вид:
/ ч 10, г < 0, е(( )=<
[б,г>0.
где Б -постоянная величина; г - время.
(4)
С учетом выражений (3) и (4) запишем систему уравнений электрического равновесия нелинейной тестовой цепи при t > 0:
ип + и = Е;
I = 13 ^хр(ил / ф) -1];
(5)
I = СШи / Ж.
Из системы (5) нетрудно получить ОДУ относительно напряжения на
емкости для тестовой цепи с экспоненциальной нелинейностью:
= 18 {{КЕ - и)/Ф] - 1}.
ш
(6)
К сожалению, получить из уравнения (6) явное выражение для и(0 не удается, поэтому далее рассмотрим ОДУ нелинейной цепи относительно тока. Из формулы (3) выразим напряжение на диоде через ток
ив = ф 1п
— +1
V ^ у
(7)
Напряжение на емкости связано с током следующим образом:
и
1 1
= — 11 й .
С
(8)
Подставив (7) и (8) в первое уравнения системы (5) получаем
ф 1п
г ■ \
1-+1
^ у
1 Х
+ — I" ¡Ж = Е.
С
(9)
Продифференцировав обе части уравнения (9) по t, получаем ОДУ заданной цепи относительно тока
Ф Ш ¡ ^ ■ + — = 0.
¡ +18 Ш С
Разделяя переменные в уравнении (10), получаем
7 фС 7. ш =--1-ш .
¡(¡' +18)
Интегрируя обе части уравнения (11), определяем
(10)
(11)
—с»
t = -Ь^di = ^lni^ - «Сln^ = Tmax In
l i(i + ^) is
i h
( I \ 1+
V i J
+ io, (12)
где т =
^ max
фС
77
; ^ = -Tmax ln
' I ^ 1 + -
V i0 J
; i0 = i(0) - значение тока цепи при t = 0.
Из выражения (12) нетрудно определить ток:
i(t) =
1
exp[(t - t0)/ тmax] - 1
(13)
Подставляя выражение (13) в (7) и учитывая, что и = Б - ид, получаем решение ОДУ нелинейной цепи в явном виде:
1
u (t) = E -ф ln
exP[(t - t0)/ Tmax ] - 1
+ 1} = E + фln{1 -exp[-(t-10)/тmaX]}. (14)
Значение t0 можно определить из выражения (12), полагая t = 0 и учитывая, что u(0) = u0, i(0) = IS {exp[(E - u0)/ф] -1} (см. формулу (3)):
'0 = -Tmax ^ + Is {exp[(E-Su0)/ф] - 1}} = Tmax ^ + Ч^-"0 / Ф) " ^ / 4
Для нелинейных цепей, как и для нелинейных ОДУ, вместо понятия «постоянная времени» используют термин «переменная временная постоянная» ~(t), которая для ОДУ первого порядка du/dt = f(u) определяется следующим образом [10]
1
т (t) = -
(15)
df / du
Учитывая, что правая часть тестового ОДУ (6) f (и) = 15 {ехр[(Б - и)/ф] -1} получаем с помощью формулы (15) «переменную временную постоянную» нелинейной тестовой цепи
s фС т (t) = — exp
IS
fu - Eл v ф j
= т max exp
E + фln{1 - exp[-(t - 10)/тmax ]} E
ф
(16)
тmax {1 - exp[-(t - t0)/тmax]}.
0
Жесткость нелинейной тестовой цепи можно оценить отношением максимального и минимального значений ~ (г) [10]:
П = Т тах / Т тт , (17)
ГДе Ттах = = ФС / ; Ттт = ~(0) = Ттах[1 - еХР(^0 / Ттах)].
Для изменения жесткости нелинейной цепи можно варьировать начального напряжения на емкости и0. На рис. 2, а показано напряжение на выходе нелинейной тестовой цепи, рассчитанное по формуле (14) при различных значениях и0 (кривая 1 соответствует и0 = 0,9 В, кривая 2 -и0 = 0,96 В, кривая 3 - и0 = 0,99 В), штриховая линия на этом же рисунке -напряжение на выходе линейной ЯС-цепи описывающееся следующим
выражением и(г) = Е - (Е - и0)е
-г / тС
Рис. 2. - Временные диаграммы напряжений на выходе тестовых цепей (а) и оценка жесткости нелинейной тестовой цепи (б)
Зависимость жесткости нелинейной тестовой цепи от начального напряжения на емкости приведена на рис. 2, б, из которого следует, что при приближении и0 к нулю значение жесткости может достигать 1016. При проведении расчетов были выбраны следующие параметры тестовых цепей: С = 1 Ф; Е = 1 В; Я = 1 Ом (для линейной цепи); ^ = 0,027 А, ф = 0,027 В (для нелинейной цепи). При этом тС = ттах = т = 1 с.
Важнейшим критерием эффективности численного метода является его точность. Таким образом, для оценки эффективности методов численного
анализа цепей во временной области (численных методов решения ОДУ) рассчитывают текущую относительную погрешность численного решения, которая имеет следующий вид
= U(tk ) - uk\/Umax , (18)
где и* -значение отклика цепи, полученное с помощью численного метода; u(t*) - точное значение отклика; umax - максимальное значение отклика.
Для проверки работоспособности предлагаемого в данной работе теста было проведено исследование точности следующих методов численного анализа переходных процессов, использующихся в пактетах математического и схемотехнического моделирования: метод трапеций (TR); методы Гира (Gear) и BDF, основанные на формулах дифференцирования назад; метод Radau [4 - 6]. Следует отметить, что метод Гира используется по умолчанию для анализа переходных процессов в большинстве программ схемотехнического моделирования. Методы BDF отличается от метода Гира только процедурой автоматического выбора шага интегрирования и реализован в программах компьютерной математики.
В данной работе тестирование методов TR и Gear проводилось в пакете NI Multisim, методов BDF и Radau - в программе Mathcad. На рис. 3 приведены текущие погрешности анализа переходных процессов в линейной и нелинейной тестовых цепях, параметры которых были указаны выше. Результаты на рис. 3 были получены при использовании установленных по умолчанию опций численного анализа: предельно допустимая относительная
- з
ошибка моделирования RELTOL = 10 ; временной шаг выбирается программой автоматически.
Из рис. 3 видно, что текущая погрешность моделирования для метода трапеций существенно возрастает при решении жесткой нелинейной задачи по сравнению с аналогичной погрешностью для линейной задачи. Данный факт подтверждает корректную работу предлагаемого теста, поскольку
известно, что метод трапеций теряет свою эффективность при решении жестких задач [5, 10].
Рис. 3. - Текущие погрешности моделирования линейной (а) и нелинейной
(б) тестовых цепей
Погрешности методов Gear и BDF существенно отличатся друг от друга для задач различной жесткости, что очевидно связано с различной процедурой автоматического выбора шага в методах Gear и BDF. Таким образом, предлагаемая тестовая задача позволяет сравнивать как точность разностных схем численных методов, так и эффективность алгоритмов автоматического выбора шага в программах численного анализа. При решении жесткой задачи на интервале t > т текущая погрешность моделирования для метода Radau, как и для метода трапеций, существенно возрастает по сравнению с аналогичной погрешностью для линейной задачи. По-видимому, данное свойство, а также низкая точность метода Radau при решении колебательных задач [6] препятствуют его применению в пакетах схемотехнического моделирования.
Таким образом, представленная тестовая задача и методика тестирования позволяют оценить как достоверность численного анализа, так и его точность. Предлагаемая тестовая задача также позволяет оценить влияние характеристик исследуемой цепи или системы (нелинейность, жесткость, скорость затухания переходных процессов) на свойства
10 102 )03 tlx
б
а
численного метода и установить предельную точность численных методов при решении практических задач.
Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 16-07-00631a.
Литература
1. Бегляров В.В., Берёза А.Н. Гибридный эволюционный алгоритм решения систем линейных алгебраических уравнений, описывающих электрические цепи // Инженерный вестник Дона, 2013, № 1. URL: ivdon.ru/ru/magazine/archive/n1y2013/1540.
2. Маничев В.Б., Жук Д.М., Витюков Ф.А. Метод математического тестирования программ анализа переходных процессов в САПР электронных схем // Проблемы разработки перспективных микро- и наноэлектронных систем (МЭС). 2014. № 1. С. 83-88.
3. Душенин Д.Ю. Численное моделирование функционирования ансамбля нейронов коры головного мозга // Инженерный вестник Дона, 2011, № 4. URL: ivdon.ru/ru/magazine/archive/n4y2011/533.
4. Hairer E., Wanner G. Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems. 2nd rev. ed. Springer-Verlag, Heidelberg, 1996. 614 p.
5. Калиткин Н.Н., Численные методы решения жестких систем // Математическое моделирование. 1995. Т. 7, № 5. С. 8-11.
6. Бирюков В.Н., Пилипенко А.М. Анализ погрешности численного моделирования автогенераторов во временной области // Известия ЮФУ. Технические науки. 2014. № 11 (160). С. 119-127.
7. Гужев Д.С., Калиткин Н.Н., Уравнение Бюргерса - тест для численных методов // Математическое моделирование. 1995. Т. 7, № 4. С. 99-127.
8. Dahlquist G.G. A special stability problem for linear multistep methods // BIT Numerical Mathematics. 1963, Vol. 3, № 1. pp. 27-43.
9. Jaeger R.C., Blalock T.N. Microelectronic circuit design. Fourth Edition. McGraw-Hill, New York, 2011. 1365 p.
10. Modern Numerical Methods for Ordinary Differential Equations / edited by G. Hall, J. M. Watt, Clarendon Press, Oxford, 1976. 336 p.
References
1. Beglyarov V.V., Bereza A.N. Inzenernyj vestnik Dona (Rus), 2013, № 1. URL: ivdon.ru/ru/magazine/archive/n1y2013/1540.
2. Manichev V.B., Zhuk D.M., Vitukov F.A. Problemy razrabotki perspektivnykh mikro- i nanoelektronnykh sistem (MES). 2014. № 1. pp. 83-88.
3. Dushenin D.Yu. Inzenernyj vestnik Dona (Rus), 2011, № 4. URL: ivdon.ru/ru/magazine/archive/n4y2011/533.
4. Hairer E., Wanner G. Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems. 2nd rev. ed. Springer-Verlag, Heidelberg, 1996. 614 p.
5. Kalitkin N.N. Matematicheskoe Modelirovanie. 1995. Vol. 7, № 5. pp. 8-11.
6. Biryukov V.N., Pilipenko A.M. Izvestiya SFedU. Engineering sciences. 2014. № 11 (160). pp. 119-127.
7. Guzhev D.S., Kalitkin N.N. Matematicheskoe Modelirovanie. 1995. Vol. 7, № 4. pp. 99-127.
8. Dahlquist G.G. BIT Numerical Mathematics. 1963, Vol. 3, № 1. pp. 27-43.
9. Jaeger R.C., Blalock T.N. Microelectronic circuit design. Fourth Edition. McGraw-Hill, New York, 2011. 1365 p.
10. Modern Numerical Methods for Ordinary Differential Equations. Edited by G. Hall, J. M. Watt, Clarendon Press, Oxford, 1976. 336 p.