CC BY
61
УДК 621.396.6:004.942 А.В. Ненашев
Сходимость метода последовательных приближений при моделировании нелинейных радиотехнических устройств. Приближенный метод Ньютона
Рассмотрены вопросы сходимости метода последовательных приближений при моделировании нелинейных радиотехнических устройств. Предложена приближенная реализация метода Ньютона, проведено сравнение его результатов с методами простой итерации и гармонического баланса.
Ключевые слова: моделирование, нелинейные устройства, последовательные приближения, метод Ньютона.
Методы последовательных приближений (итерационные) широко применяются для численного решения различных уравнений и их систем. При схемотехническом моделировании нелинейных радиоэлектронных устройств (РЭУ) к ним относится, например, метод гармонического баланса (ГБ). В [1] автором предложен новый метод последовательных приближений для моделирования таких устройств, при котором математические операции производятся в пространстве функционалов, а нелинейные элементы представляются как параметрические. В [2] рассмотрены вопросы сходимости при его реализации в виде метода простой итерации. Однако на практике в большинстве случаев применяют итерационный метод Ньютона, обеспечивающий более быструю сходимость.
Постановка задачи: рассмотреть возможность моделирования на основе итерационного метода Ньютона, провести сравнение его сходимости с другими алгоритмами.
Разработанные в литературе методы последовательных приближений [3], как правило, предназначены для решения нелинейных (трансцендентных) уравнений, в состав которых входят функции действительного переменного, например времени. В частности, в системах схемотехнического проектирования они используются для определения статических режимов устройств по постоянному току, т.е. без учёта инерционных свойств. При моделировании динамических режимов работа устройств описывается интегродифференциальными уравнениями. Для их решения также можно применить методы последовательных приближений, если оперировать в пространстве функционалов. При этом не возникает серьёзных препятствий для реализации метода простой итерации, но он не всегда обеспечивает достаточно быструю сходимость.
Рассмотрим принципы построения итерационного метода Ньютона. Он является развитием метода простой итерации, основы которого изложены в [2]. Пусть имеется уравнение, записанное в общем виде: ф(х) = 0 , где ф(х) - функция, определённая на отрезке [а,Ь]. Известна достаточно малая область, в которой располагается решение этого уравнения: х = X . В ней выбирается начальное приближение х(0), достаточно близкое к искомому корню. Выбирают функцию у(х) так, чтобы
корень х = х* совпадал с корнем уравнения х = у(х). Затем строится последовательность х(0),х(1),--.,х(и),... через соотношение х(п) =у(х(и-1)), образующее итерационный процесс. Здесь и далее нижний индекс в скобках означает номер вычисляемого приближения. Сходимость обеспечивается выбором вида функции у(х) и начального приближения х(0). Методу простой итерации
соответствует выбор функции у( х) в соответствии с равенством у( х) = х -8-ф( х), где 8 - кон -станта.
Для метода Ньютона эта функция выбирается иначе:
ф( х)
у( х) = х-8-
ф'( х)
где 8< 1. Тогда итерационная формула будет иметь следующий вид:
х х 5Ф(х(«-1)) (1)
Х(п) = Х(п -1) 8 ,, ч- (1)
ф( Х(п-1))
Будем считать коэффициент 8 постоянной величиной, хотя существуют варианты алгоритма, при которых он принимает различные значения на разных шагах вычислений. Если ф(х) действительная функция, монотонная вблизи корня х = х* , производная ф'(х) не обращается в нуль и (1) непосредственно можно использовать для вычислений. В более сложных случаях, например для систем уравнений, это невозможно, поэтому записывают иначе, в векторной форме:
Ф'(Х(И_1)) • АХ(И) = -8 • Ф(Х(И_1)), (2)
где АХ(И) =(Х(И)) —(Х(п-1)) - вектор приращений, ф(Х(и—)) играет роль вектора невязок, а Ф'(Х(и-1)) - матрица Якоби, составленная из частных производных. В частности, при моделировании нелинейных РЭУ методом ГБ, (2) строится как система линейных алгебраических уравнений, составленных для амплитуд всех гармоник тока или напряжения во всех узлах схемы. Поэтому способы ускорения моделирования методом ГБ сводятся к сокращению объёма вычислений при решении системы линейных уравнений.
При реализации метода моделирования [1, 2], предложенного автором, ф(х) представляет собой функционал, построенный на основе интегродифференциального уравнения, описывающего устройство. Непосредственное дифференцирование функционала невозможно, допустимо лишь вычисление дифференциала Гато или Фреше [4]. На практике используют дифференцирование Гато, как более просто реализуемое, имея в виду, что в достаточно малой окрестности корня уравнения производные Гато и Фреше совпадают. В этом случае система уравнений (2) примет вид
аФ(Х) = -8-Ф(Х), (3)
где dФ(X) = Ф'(Х)- АХ - вектор дифференциалов Гато от функций Ф(Х). Очевидно, что решение такой системы уравнений относительно АХ обычными способами невозможно.
Рассмотрим пример схемы диодного преобразователя частоты, для которого в [2] применён метод простой итерации. Там же приведена его упрощенная схема. Работа преобразователя описывается уравнением
г
Ыд (г) = еэкв (г) - {Ыд (т) • Яд [Ыд (т)] • (г - т) • Й т ,
0
где Ыд (г) - напряжение на диоде; Яд = я д (Ыд) - проводимость диода, еэкв (г) - эквивалентная эдс, включающая напряжения сигнала и гетеродина; 2^ (г) - импульсная характеристика суммарного сопротивления, включенного последовательно с диодом. Ёмкостью диода пренебрегаем. Тогда функция ф(г) запишется в виде
г
ф[Ыд (г)]=Ыд(г)+|Ыд(т)Яд[Ыд(т)}2Е(г—'тМт-еэкв(г) .
0
Соответственно, система (3) сведётся к одному уравнению
Йф(х)=-8^ф(х), (4)
где Йф(х) = ф'(х)•Ах - дифференциал Гато функции ф(х). Согласно правилам дифференцирования по Гато [4], он будет иметь вид
Йф[ыд (г)] =Амд (г)+!-|яд[Ыд(т)] + Ыд (тЭ-Яд-Ы^^Амд(т>(г-т><1т .
Тогда уравнение (4) запишется следующим образом:
г
Амд (?) +1 |яд [мд (т)] + Мд (т)-^ Яд [мд (т)]|-Амд (т)-гТ.(?-тМт = 5 ■
<?экв (?) - мд (0 - ^д (Ф-^Е(?-хМт
0
где для краткости сделана замена ід(?) = Мд(?)-Яд[мд(?)] . Для дальнейших операций перейдём в область многомерного преобразования Лапласа (МПЛ), как рекомендовано в [5]:
(5)
Шд ( 51,^2,-) + ^ (s1)•АUд (52,53,---)7'Е(51+s2 +53+--0 + Uд (s1)•GД Ы х хАид (s3, ^г-^Е (51+52 +53+- - -)^{^эквС^ -Uд (51)—^дС^^Е^].
Здесь неизвестная величина АЦд (sl ,S2,—) записывается как многомерное изображение, величины Uд (5), Gд(s), GД (5) считаются известными по результатам предыдущего цикла вычислений, как и Eэкв (5), поэтому записываются как одномерные. Многомерное изображение искомой величины представляется как сумма слагаемых различных порядков:
Аид (5 ^53,-) = Аид) (51) + Аид2; (51,52) + Аид3; (5,52,53) + - .
(2)(
к3)(
(6)
Для нахождения неизвестного приращения последовательно приравниваются его слагаемые, начиная с первого порядка, и известные слагаемые таких же порядков в правой части. При этом произведения, в которые входят уже найденные порядки, переносятся вправо. Таким образом, из (5) получим, обозначив одновременно номера итераций:
Аи $,)( •*) =8[Еэкв (*)-/д (п-фУг Е(*)-и д(„-!)(*)],
Аи &(,) ={-8 ^(п-1)(*) { {2> {(*1+^2)+(1-8)^Аи ^Ц*,
Аи ді)(5) ={-5
Аи д4(П)( 5) ={-5
G„
д (п—1) (51) •Аи д2(П) (52,53) + и д(и—1)(51>Ц (и-1)(52) х
хАи д(И) (s3)}•ZS(s1+52 +53)
+а—5>Аи д2(П)( 51)
(7)
^(п—1)(51) Аи д3(П) (52,53,54) + и д(п—1)(51) • ^д (п—1)(52) х
хАид^И)(53,54)
7х(51+52+53+54) + (1—5) •АидИ)(51)} ,
Здесь символом {о}* обозначена операция перехода от многомерного преобразования Лапласа к
одномерному. Такую операцию удобнее производить во временной области [5], заменяя все переменные времени одной. Также во временной области производится перемножение переменных (нелинейных) величин. Линейные операции, например, умножение на комплексные коэффициенты или сопротивления, дифференцирование, интегрирование производятся в частотной области. Если требуется моделирование устройства в установившемся режиме, от МПЛ следует перейти к многомерному преобразованию Фурье заменой комплексных аргументов 5^ = . При вычислениях необхо-
димо постоянно переходить из временной области в частотную и обратно, что делают с помощью БПФ. Это характерно и для других методов моделирования, например ГБ.
По найденному приращению вычисляют следующее приближение к напряжению на диоде из соотношения
ид(и)(5) = ид(и-1)(5) + {Аид(и)(51,52,52,”•)} . (8)
В (7) во всех выражениях, кроме первого, присутствует слагаемое, равное приращению предыдущего порядка, умноженному на (1 — 5). Это связано с тем, что вычисление суммарного приращения (6) фактически реализуется методом последовательных приближений и подчиняется алгоритму простой итерации, описанному в [2]. Конкретное значение коэффициента 5 выбирается так, чтобы обеспечить сходимость при заданной амплитуде воздействия, и в то же время добиться требуемого результата за наименьшее число итераций. Аналитического критерия для выбора 5 , по-видимому, в общем случае предложить невозможно. На практике следует использовать метод подбора на основе пробных вычислений.
Требуемое число слагаемых различных порядков в (6), вычисляемых в соответствии с (7), определяется необходимой точностью решения. При этом процесс состоит обычно из довольно большого числа итераций. Поэтому для сокращения общего количества вычислений можно на начальных
*
циклах вычислять малое число слагаемых, увеличивая его в каждой следующей итерации. Таким образом, можно записать общую итерационную процедуру, с начального, линеаризованного приближения:
Uд(0}(5) = Еэкв
(s)/[1+Сд0 Z2 (s)] ,
где Сд о - начальная проводимость диода. Приращение ДЦд (s) на этом шаге будет состоять из одного слагаемого первого порядка:
Дид(1) (s) = Ди®} (s) =8 • [Кэкв (s)-^{0)(s) ^(sb^s)].
Первое приближение к напряжению на диоде определяется как сумма U^q(s) = ^(0)(s) + ДЦд(1)(£), затем вычисляется следующее приращение, уже из слагаемых двух порядков, вычисляемых в соответствии с (7): ДЦд(2)(£)={дид^^О+Дид^)^^)} . Следующее приращение состоит уже из слагаемых трёх порядков: Дид(3)(я)={Ди^(я^+Дид^)^^) +Дид^^)(si,s2,s3)} и т.д. Соответствующее приближение к напряжению на диоде вычисляется по формуле (8).
Из теории известно [3], что метод простой итерации обеспечивает сходимость со скоростью геометрической прогрессии, т.е. остаточная погрешность уменьшается на каждом шаге в одинаковое число раз. В то же время метод Ньютона обеспечивает квадратичную скорость сходимости. Однако коэффициенты при этих законах могут иметь различную величину, поэтому реальный выигрыш можно оценить только сравнительными расчётами для конкретных примеров. В соответствии с приведёнными соотношениями была разработана программа и проведены расчёты, результаты которых сравнивались с результатами, полученными методом простой итерации [2]. Одновременно производилось сравнение с методом ГБ, широко применяемым на практике, например в пакете программ Microwave Office. Его существенное отличие от рассматриваемого в данной работе в том, что составляется система уравнений, каждое из которых соответствует амплитуде одной гармоники тока или напряжения. Размерность системы уравнений получается значительной, для её решения применяется метод Ньютона, с различными алгоритмами сокращения вычислений. Для расчётов по этому методу также была составлена компьютерная программа, с использованием принципов, описанных в [6]. Результаты, полученные тремя методами, отражены в таблице. Там же приведены результаты, полученные для схемы усилительного каскада, построенного на полевом транзисторе с барьером Шоттки [1]. Для детального сравнения методов необходимо проводить оценку сложности программ, общего количества операций и требуемого объёма памяти.
Количество итераций, необходимое для получения заданной точности
Моделируемое устройство Диодный преобразователь частоты Каскад усилителя на ПТШ
Метод моделирования Погреш- ность Степень нелинейности
Высокая Средняя Умеренная Умеренная
Простой итерации 0,01 96 50 19 12
0,001 148 72 28 20
Г армонического баланса 0,01 58 32 12 9
0,001 88 41 18 13
Приближенный метод Ньютона 0,01 13 11 5 5
0,001 18 13 7 7
В настоящей работе предполагалось, что поскольку все методы итерационные, обобщённое суждение об их эффективности можно сделать на основе сведений о количестве итераций, которое требуется для достижения заданной погрешности расчётов.
Степень нелинейности устройства приближенно оценивалась по уровню (амплитуде) воздействия. Усилительный каскад на ПТШ исследовался при уровне входного воздействия 10 дБм, что примерно на 1 дБм превышает точку сжатия и соответствует умеренной нелинейности. Диодный преобразователь частоты исследовался при амплитудах гетеродинного напряжения 1,5 В, 1 В и
0,7 В, что было оценено как режимы высокой средней и умеренной нелинейности соответственно. По полученным результатам можно сделать следующие выводы:
— метод гармонического баланса во всех случаях даёт сокращение количества циклов, относительно метода простой итерации, менее чем в 2 раза. Вопрос требует более подробного изучения, но этот факт ставит под сомнение целесообразность широкого применения метода ГБ, поскольку вычисление матрицы Якоби и решение системы уравнений, особенно для громоздких схем, неизмеримо сложнее, по сравнению с методом простой итерации;
— приближенный метод Ньютона при моделировании в пространстве функционалов даёт существенное сокращение числа циклов, по сравнению с методом простой итерации, причём выигрыш возрастает с увеличением степени нелинейности. В случае диодного преобразователя частоты выигрыш достигает 8 раз;
— по сравнению с методом ГБ приближенный метод Ньютона также даёт существенное сокращение количества циклов (максимально более чем в 4 раза) что, несомненно, связано с тем, что уравнения имеют интегральный характер, описывают процессы в целом, а не поведение отдельных гармоник, как это имеет место в методе ГБ.
Заключение. Рассмотрено применение метода Ньютона, предложен приближённый алгоритм его реализации при моделировании нелинейных РЭУ в пространстве функционалов. Проведено сравнение с методом простой итерации и традиционным методом ГБ. Полученные результаты свидетельствуют о перспективности предложенного метода.
Литература
1. Ненашев А.В. Метод моделирования нелинейных радиотехнических устройств // Доклады Томского государственного университета систем управления и радиоэлектроники. - 2010. - № 1 (21), ч. 2. - С. 50-54.
2. Ненашев А.В. Сходимость метода последовательных приближений при моделировании нелинейных радиотехнических устройств. Метод простой итерации // Доклады Томского государственного университета систем управления и радиоэлектроники. - 2010. - № 2 (22), ч. 1. - С. 244-248.
3. Березин И.С. Методы вычислений / И.С. Березин, Н.П. Жидков. - М.: Физматгиз, 1962. -Ч. 2. - 620 с.
4. Колмогоров А.Н. Элементы теории функций и функционального анализа / А.Н. Колмогоров, С.В. Фомин. - М.: Физматгиз, 1976. - 543 с.
5. Ненашев А.В. Применение многомерного преобразования Лапласа для анализа параметрических и нелинейных процессов и устройств // Радиотехника. - 2006. - № 12. - С. 18-21.
6. Актуальные проблемы моделирования в системах автоматизации схемотехнического проектирования / Отв. ред. А.Л. Стемпковский. - М.: Наука, 2003. - 430 с.
Ненашев Александр Васильевич
Канд. техн. наук, доцент каф. радиоэлектроники Сургутского государственного университета ХМАО-Югры Тел.: 8-912-810-44-75 Эл. почта: navas1@mail.ru
Nenashev A.V
Convergence of successive approximation method for modeling nonlinear radio engineering devices. Newton’s approximate method
The article considers the questions of convergence of successive approximation method in modeling nonlinear radio engineering devices. There have been offered the approximate implementation of Newton’s method, the results have been compared to the results of fixed point iteration method and describing function method. Keywords: modeling, a nonlinear device, successive approximation method, Newton’s method.
