Содержание к диссертации
Введение
Глава 1: Обзор известных методов анализа жестких радиотехнических цепей 15
1.1 Жесткие задачи анализа радиотехнических цепей 15
1.2 Точность и устойчивость численных методов решения систем обыкновенных дифференциальных уравнений 25
1.3 Жесткость систем обыкновенных дифференциальных уравнений 28
1.4 Современные численные методы решения жестких систем 30
1.5 Методы идентификации параметров моделей компонентов радиотехнических цепей 35
1.6 Основные результаты и выводы 42
Глава 2: Проблемы численного интегрирования уравнений жестких радиотехнических цепей и оптимизации параметров моделей компонентов 43
2.1 Исследование неявных методов Рунгс-Кутта 43
2.2 Исследование методов Розенброка, трапеций и Гира 48
2.3 Критический анализ известных методов параметрической идентификации 54
2.4 Основные результаты и выводы 60
Глава 3: Разработка эффективных численных методов анализа жестких радиотехнических цепей с сосредоточенными параметрами 61
3.1 Метод экспоненциальной нолгоики 62
3.2 Метод локальной подгонки 64
3.3 Метод глобальной подгонки 71
3.4 Исследование метода глобальной подгонки при решении жесткой нелинейной задачи 77
3.5 Основные результаты и выводы 81
Глава 4: Разработка эффективных численных методов анализа цепей с распределенными параметрами 82
4.1 Решение неявных разностных схем Рунге-Кутта итерациями Зейделя 83
4.2 Решение неявных разностных схем Рупге-Кутта блочными итерациями Зейделя 97
4.3 Основные результаты и выводы 103
Глава 5; Идентификация параметров жестких моделей компонентов радиотехнических цепей 105
5.1 Идентификация параметров статических моделей полупроводникового диода 106
5.2 Алгоритм оценки жесткости задач параметрической идентификации 115
5.3 Исследование обусловленности моделей диода 120
5.4 Идентификация параметров моделей барьерной емкости -и-перехода 128
5.5 Идентификация параметров статических моделей полевого транзистора 132
5.6 Основные результаты и выводы 142
Заключение 144
Литература 146
Приложение 152
- Точность и устойчивость численных методов решения систем обыкновенных дифференциальных уравнений
- Исследование методов Розенброка, трапеций и Гира
- Метод глобальной подгонки
- Решение неявных разностных схем Рупге-Кутта блочными итерациями Зейделя
Введение к работе
Актуальность работы. Известные программы схемотехнического проектирования используются чаще всего для проверки работоспособности радиотехнических цепей. Однако возможности современной вычислительной техники позволяют выполнить численными методами более детальные исследования, сопоставимые с результатами физического макетирования. Для этого необходимо учитывать большое число параметров радиотехнической цепи. Как показано в книге Ю. В. Ракитского [1], учет большого числа параметров при построении математических моделей радиотехнических цепей часто требует для полного описания процессов на любом отрезке времени функций двух типов, убывающих быстро и медленно. Функции первого типа убывают очень быстро, так что большую часть времени наблюдаются функции второго типа, которые убывают очень медленно. Однако, в любой момент времени сохраняется возможность появления быстрозатухающего процесса, описываемого функциями первого типа. Такое поведение радиотехнической цепи называется жесткостью, а системы обыкновенных дифференциальных уравнений (СОДУ), моделирующие цепи такого типа, называются жесткими СОДУ.
Жесткость наиболее часто проявляется в следующих задачах анализа радиотехнических цепей, которые рассмотрены в работах Л. О. Чуа [2], В. В. Дикусара |3j, В. Б. Дмитриева-Здорова |4|:
При анализе цепей с сосредоточенными параметрами, таких как генераторы гармонических сигналов, широкополосные усилители с обратными связями, устройства с нелинейными элементами, имеющими высокую вариацию дифференциальных параметров.
При анализе цепей с распределенными параметрами, представленных в виде многозвенных структур с сосредоточенными параметрами.
При идентификации параметров моделей компонентов радиотехнических цепей.
Решение первой и второй задачи состоит в численном интегрировании
5 СОДУ цепи. В работах Н. Н. Калиткина [5, 6] показано, что современные численные методы интегрирования СОДУ при решении жестких задач имеют ряд недостатков, связанных с необходимостью использовать мелкий шаг интегрирования, что приводит к росту времени анализа и погрешности округления. В результате в программах автоматизированного анализа электронных цепей (иногда в одной w той же программе) используются различные методы решения СОДУ. Чаще всего - метод Гира, метод трапеций и неявный метод Эйлера (как наиболее надежный при решении жестких задач). При моделировании указанных выше цепей, которые далее будем называть жесткими, методы Гира и трапеций теряют свою эффективность из-за особенностей их устойчивости, а неявный метод Эйлера дает весьма приближенные результаты, поскольку имеет минимальный порядок точности.
Решение третьей задачи состоит в оптимизации (минимизации) многомерных функционалов, имеющих, так называемую, овражную структуру. Задача минимизации таких функционалов является жесткой. Овражные функционалы характеризуются наличием некоторой области притяжения, в которой норма вектора-градиента минимизируемого функционала существенно меньше, чем в остальной части пространства. Практически все методы оптимизации, используемые в программах анализа цепей, позволяют достаточно быстро попасть в область притяжения из начальной точки, но далее процесс спуска резко замедляется и останавливается в некоторой точке, расположенной не обязательно вблизи точки истинного минимума. Кроме того, в результате недифференцируемости минимизируемого функционала вблизи дна оврага, потенциальная минимальная погрешность метода наименьших квадратов не может быть реализована, причем алгоритмы, использующие численное дифференцирование, оказываются наименее точными. Это связано с тем, что на практике характеристики компонентов цепей измеряются с погрешностью, намного превышающей погрешность представления чисел в компьютере, и при этом на дне овражной структуры наблюдается цифровой шум.
Фактически эффективность алгоритмов решения задачи оптимизации
6 падает не с увеличением размерности задачи, а с увеличением ее жесткости которая, как правило, растет с увеличением размерности, что показано в работе Р. П. Федоренко [7].
Исходя из всего выше сказанного, очевидно, что методы анализа радиотехнических цепей в существующих пакетах автоматизированного проектирования либо имеют ряд существенных недостатков, либо вообще непригодны при анализе жестких моделей. Поэтому требуется выявить недостатки известных методов и разработать новые методы анализа жестких цепей, более эффективные, то есть более точные и устойчивые, по сравнению с известными методами.
Значительный вклад в развитие теории методов решения жестких задач анализа радиотехнических цепей внесли И. Г. Черноруцкий, В. Н. Ильин, Я. К. Трохименко, В.Б.Дмитриев-Здоров, а также Ю. В. Ракитский, Н. Н. Калиткин, Р. П. Федоренко. Работы по данной теме ведутся отечественными исследователями, и интенсивно продолжаются за рубежом, однако существует много нерешенных проблем связанных с использованием известных методов для автоматизированного анализа жестких моделей, что говорит об актуальности темы диссертационной работы. Новые методы решения жестких СОДУ и алгоритмы параметрической идентификации для задач детального анализа радиотехнических цепей необходимы для повышения эффективности разработки радиоэлектронных устройств при полной или частичной замене натурного эксперимента развитым численным моделированием.
Целью работы является создание и обоснование новых эффективных методов и моделей для решения с гарантированной точностью жестких задач теории радиотехнических цепей.
В соответствии с поставленной целью в диссертационной работе решались следующие задачи;
1. Критический анализ эффективности (точности и устойчивости) известных численных методов решения СОДУ и идентификации параметров жестких моделей радиотехнических цепей и их компонентов.
Разработка и обоснование более эффективного, чем известные методы, численного метода решения СОДУ, предназначенного для анализа жестких целей с сосредоточенными параметрами.
Формулировка и обоснование более эффективных, чем известные методы, численных методов решения СОДУ высокого порядка, предназначенных для моделирования цепей с распределенными параметрами.
Разработка численных методов оценки жесткости и обусловленности моделей компонентов радиотехнических цепей.
Обоснование алгоритма оптимизации и моделей компонентов радиотехнических цепей, позволяющих уменьшить влияние жесткости при идентификации параметров моделей.
Основные положения и результаты, выносимые на защиту.
Результаты исследования точности и устойчивости современных методов численного анализа жестких радиотехнических цепей.
Модифицированные численные методы анализа жестких цепей с сосредоточенными параметрами и результаты исследования их точности и устойчивости при решении тестовых задач соответствующих жестким моделям линейных и нелинейных радиотехнических цепей.
Модифицированные численные методы анализа цепей с распределенными параметрами и результаты исследования их точности и устойчивости при решении тестовых задач соответствующих различным моделям звеньев с сосредоточенными параметрами, аппроксимирующих распределенную структуру.
Методы оценки жесткости и обусловленности моделей компонентов радиотехнических цепей.
Алгоритм оптимизации и результаты исследования его точности на известных и скомпилированных в данной работе статических и динамических моделях диода и полевого транзистора.
Научная новизна.
1. Разработан модифицированный метод решения СОДУ, предназначенный для анализа жестких цепей с сосредоточенными параметрами, разност-
8 ная схема которого учитывает все исходные данные для задачи Коши и обеспечивает квадратичную зависимость глобальной погрешности от временного шага.
2. Предложен метод решения жестких СОДУ высокой размерности (мо делирующих распределенные цепи с малыми потерями), построенный на осно ве итераций Зейделя.
3. Разработан алгоритм, позволяющий количественно оценить жест кость моделей компонентов радиотехнических цепей. Впервые определены жесткость и обусловленность ряда моделей полупроводникового диода и полевого транзистора.
4. Предложен метод параметрической идентификации моделей компонентов радиотехнических цепей, позволяющий уменьшить влияние жесткости па поиск минимума многомерной целевой функции и тем самым повысить точность определения параметров моделей.
Практическая значимость.
Разработанный метод численного анализа жестких цепей с сосредоточенными параметрами, обеспечивает увеличение порядка точности по сравнению с известными методами при той же вычислительной сложности.
Предложенные методы решения жестких СОДУ высокой размерности для анализа цепей с распределенными параметрами обеспечивают линейный рост вычислительных затрат с ростом размерности задачи при сохранении необходимой точности и устойчивости решения.
Результаты исследования жесткости моделей компонентов радиотехнических цепей позволяют обосновать выбор методов идентификации параметров для обеспечения требуемой точности определения параметров.
Результаты исследования обусловленности моделей полупроводникового диода позволяют сформулировать требования к точности и области измерения характеристик прибора при заданной точности определения параметров или при известной точности и области измерений получить оценку точности определения параметров.
5. Предложенный метод параметрической идентификации, позволяет определять параметры жестких моделей компонентов радиотехнических цепей с гарантированной точностью не менее 4 - 5-ти значащих цифр.
Методы исследовании. В работе использованы методы теории дифференциальных уравнений, вычислительной математики, теории цепей, математической статистики, линейного и нелинейного программирования.
Внедрение результатов работы. Основные результаты диссертационной работы нашли применение в разработках, выполненных в интересах предприятия «Научно-исследовательская лаборатория автоматизации проектирования» (г. Таганрог), при выполнении госбюджетных работ № 11170, №11055/1 на кафедре ТОР ТРТУ, а также используются в учебном процессе кафедры ТОР ТРТУ в курсах «Математические основы моделирования цепей и сигналов» и «Основы компьютерного проектирования и моделирования РЭА».
Апробация работы. Результаты работы докладывались и обсуждались: на семинарах кафедры теоретических основ радиотехники ТРТУ; на ежегодных научно-технических конференциях профессорско-преподавательского состава, аспирантов и сотрудников ТРТУ, 2002 - 2005 гг. на 6-й и 7-й Всероссийских научных конференциях «Техническая кибернетика, радиоэлектроника и системы управления». Таганрог. ТРТУ, 2002, 2004 г.г. па международной научной конференции «Анализ и синтез как методы научного познания», Таганрог: ТРТУ, 2004, на семинаре НИИ многопроцессорных вычислительных систем. Таганрог, 2004 г.
Публикации. По результатам диссертации опубликовано 16 печатных работ, в том числе 4 в центральной печати.
Личный вклад автора. Все включенные в диссертацию результаты получены в основном лично соискателем, либо при его непосредственном участии. В работах, выполненных в соавторстве, участие автора состояло в постановке экспериментов, анализе и обработке экспериментальных данных, разра-
10 ботке методов решения поставленных задач па ЭВМ.
Структура и объем диссертации. Диссертация состоит из введения, пяти глав, заключения и приложения, изложенных па 156-ти страницах и проиллюстрированных 68-го рисунками и 19-ю таблицами, а также списка литературы из 64-х наименований
Содержание работы.
Во введении обоснована актуальность работы, сформулированы цель и основные задачи работы, показана ее научная новизна и практическая значимость, выделены основные положения и результаты, выносимые на защиту.
В первой главе показано, что при анализе и оптимизации радиотехнических цепей часто приходится решать, так называемые, жесткие задачи. Рассмотрены наиболее характерные примеры таких задач.
Рассмотрены методы интегрирования СОДУ и методы оптимизации, применяемые для решения жестких задач. Отмечены известные недостатки этих методов, препятствующие их использованию или делающие их малоэффективными при анализе жестких радиотехнических цепей.
Во второй главе исследуются проблемы известных численных методов решения СОДУ и методов идентификации параметров при анализе жестких радиотехнических цепей.
Предложена методика исследования локальных и глобальных погрешностей методов интегрирования СОДУ, заключающаяся в численном решении тестовых задач различной жесткости, моделирующих радиотехнические цепи. Необходимость проведенных исследований обусловлена тем, что численное решение рассмотренных тестовых задач позволяет выявить описанные в предыдущей главе недостатки известных методов. Показано, что методы анализа жестких цепей с сосредоточенными параметрами целесообразно разрабатывать на основе простейших неявных одпошаговых методов - метода трапеций, обладающего наименьшей асимптотической локальной погрешностью и метода Эйлера, обеспечивающего /.-устойчивость решения. Для анализа цепей с распределенными параметрами, необходимы итерационные методы, позволяющие обеспечить линейный рост вычислительных затрат с ростом размерности задачи, поскольку размерность таких задач оказывается достаточно большой, а для классических неявных методов вычислительные затраты растут пропорционально квадрату размерности.
Недостатки известных методов параметрической идентификации показаны при численном анализе моделей компонентов радиотехнических цепей, используемых в программах автоматизированного анализа. Показано, что для эффективного применения методов оптимизации, при определении параметров моделей компонентов радиотехнических цепей необходимо иметь информацию о жесткости идентифицируемых моделей, то есть необходимо разрабатывать специальные алгоритмы оценки жесткости моделей. При высокой жесткости модели для обеспечения наибольшей точности идентификации целесообразно использовать наиболее медленный, но наиболее надежный в данном случае метод случайного спуска. Кроме того, для оценки точности определения параметров, необходимо исследование обусловленности, то есть зависимости параметров от диапазона и точности измерений характеристик исследуемого компонента.
В третьей главе сформулирован и обоснован более эффективный, чем известные методьі, численный метод решения СОДУ, предназначенный для анализа жестких цепей с сосредоточенными параметрами.
Рассмотрен известный, метод экспоненциальной подгонки. Этот метод имеет минимальный порядок точности и не обладает /".-устойчивостью, что препятствует его использованию для решения практических задач. Сформулирован метод локальной подгонки - одна из модификаций метода экспоненциальной подгонки, которая позволяет обеспечить L-устойчивость этого метода. Модифицированный метод локальной подгонки позволяет уменьшить погрешность неявного метода Эйлера, при сохранении /.-устойчивости метода. Но метод локальной подгонки, как и неявный метод Эйлера, имеет первый порядок точности. В связи с этим, на основе метода экспоненциальной подгонки был разработан еще один алгоритм, не имеющий недостатков метода локальной
12 подгонки - метод глобальной подгонки. L-устойчивость этого метода обеспечивается введением в разностную схему информации о величине интервала наблюдения.
Доказательство эффективности метода глобальной подгонки основывалось на решении модельной задачи различной жесткости. При малых шагах зависимость глобальной погрешности от шага для предлагаемого метода и метода трапеций имеет одинаковый наклон, то есть порядок методов одинаков, кроме того, погрешность метода глобальной подгонки монотонно убывает и меньше погрешности неявного метода Эйлера. Для дополнительного доказательства эффективности метода глобальной подгонки была рассмотрена жесткая нелинейная задача анализа работы амплитудного детектора в режиме большого сигнала. При этом было показано, что метод трапеций с пассивным сглаживанием в данном случае, в отличие от линейных задач, не обладает никакими преимуществами перед классическим методом трапеций.
Полученные результаты позволяют сделать вывод о том, что при решении жестких задач анализа радиотехнических цепей метод глобальной подгонки но точности близок к методу трапеций, а но аппроксимации быстрых компонент решения - к неявному методу Эйлера. Таким образом, метод глобальной подгонки, обеспечивает большую эффективность по сравнению с известными методами при анализе жестких цепей с сосредоточенными параметрами.
В четвертой главе исследованы эффективные численные методы решения СОДУ высокого порядка, предназначенные для моделирования цепей с распределенными параметрами. Предлагаемые методы основаны па неявных разностных схемах Рунге-Кутта при использовании для их решения итераций Зейделя, что приводит к линейной зависимости времени решения от размерности системы. Рассмотрено непосредственное применение итераций Зейделя для решения неявных разностных схем Рунге-Кутта. Приведены разностные схемы для методов первого и второго порядка. Разностные схемы для методов более высокого порядка могут быть представлены аналогично, но при этом число итераций должно быть не меньше порядка точности метода. Были рассчитаны
13 области устойчивости и глобальные погрешности методов 1 - 3-го порядков при анализе й С-лини и.
Полученные результаты [[оказывают, что методы построенные на непосредственном применении итераций Зейделя Л(а)-устойчивы только в том случае, если матрица Якоби системы близка к диагопалыю-доминирующей и, наоборот, максимально неустойчивы при нулевых коэффициентах на главной диагонали. Это означает, что непосредственное применение итераций Зейделя целесообразно только для сравнительно узкого класса RL и ЯС-цепей.
Исследован новый алгоритм, трудоемкость которого также линейно зависит от размерности задачи. Этот алгоритм обеспечивает А -устойчивость даже при наличии нулевых элементов на главной диагонали матрицы Якоби. Последовательность действий, устанавливаемая новым алгоритмом, также приводит к симметричным итерациям Зейделя, но не для дифференциальных уравнений первого порядка, как ранее, а для СОДУ второго порядка, то есть к блочным итерациям. Для доказательства эффективности нового алгоритма были рассчитаны его области устойчивости, локальные и глобальные погрешности при анализе LC-линии с малыми потерями.
Полученные результаты показывают, что метод блочных итераций позволяет указать область практического применения неявным методам Рунге-Кутта. Метод блочных итераций може г использоваться для анализа и оптимизации устройств формирования ударных волн в нелинейных средах, где малое время решения играет решающую роль. Метод блочных итераций первого порядка позволяет определять устойчивость нелинейных распределенных структур без потерь.
В пятой главе исследованы проблемы параметрической идентификации моделей компонентов радиотехнических цепей. Предложен алгоритм идентификации параметров, основанный на уменьшении размерности оптимизируемой функции и применении метода случайного спуска. Результаты идентификации параметров статических моделей диода, показывают, что уменьшение размерности минимизируемой функции приводит к уменьшению жесткости за-
14 дачи. Жесткость задачи оптимизации, и, следовательно, жесткость оптимизируемых моделей оценивалась отношением максимального и минимального значений матрицы Гессе. Алгоритм оценки жесткости заключается в переходе от исходного базиса параметров к базису, содержащему направления наиболее быстрого и наиболее медленного изменения минимизируемой функции.
Показано, что понятие обусловленности может быть использовано для описания свойств задач параметрической идентификации. Были предложены критерии для численной оценки обусловленности задач параметрической идентификации. Полученные результаты, показывают, ччо вариация параметров предложенной в данной работе компактной модели полупроводникового диода значительно меньше, чем у iSYVCis-M одели, что говорит о лучшей обусловленности компактной модели.
Компактная модель имеет тот же порядок точности, что и известная SP/CE-модель и обладает по сравнению с ней рядом преимуществ, таких как явное выражение для вольт-амперной характеристики (ВАХ) диода, значительно меньшая жесткость и более высокая обусловленность, что позволяет увеличить скорость и точность идентификации и обеспечивает слабую зависимость определяемых параметров от диапазона измерений.
Исследованы свойства двух моделей барьерной емкости /т-л-перехода и ряда статических моделей полевого транзистора, при этом использовались описанные выше алгоритмы оптимизации и оценки жесткости. Показано, что в отдельных случаях без редуцирования размерности задачи численной оптимизации идентификация параметров становится невозможной.
В заключении сформулированы основные научные и практические результаты диссертационной работы.
В приложении приведена программа идентификации параметров полевого транзистора и тестовая ВАХ полевого транзистора КПЗ07А.
Точность и устойчивость численных методов решения систем обыкновенных дифференциальных уравнений
Анализ электрических цепей преимущественно заключается в решении дифференциальных уравнений, представляющих собой математическую модель цепи. В программах автоматизированного анализа цепей предлагаются десятки различных методов численного решения систем дифференциальных уравнений.
Рассмотрим задачу Ко щи для СОДУ, записанной в векторной форме:Требуется найти решение х(7) задачи Коши. Численные методы решения задачи (1.2) позволяют найти искомую функцию x(if) приближенно в отличие от аналитических методов решения. Но при этом численные методы дают возможность приближенно решать многие из задач, точные решения которых не могут быть найдены аналитическими методами.
Численное решение задачи (1.2) ищется на интервале [;0, 7], нижняя граница которого г() задана, а верхняя граница Т может быть либо задана, либо определена в процессе расчетов некоторым условием. Искомое приближенное решение ищется в отдельных точках t\, h,.-., tn,... интервала, называемых узлами сетки, в виде последовательности векторов Х, х?,..., х„,..., приближенно равных векторам x(f), х(Ь),..., x(f„),..., определяемых точным решением \(t). Расстояние между соседними узлами называется шагом интегрирования: h = „.,. і - tп. Шаг может быть задан заранее (интегрирование с постоянным шагом) или может меняться в ходе вычислений. Вектор х п _ і может вычисляться явно:где Ф(-) — некоторая функция, зависящая от конкретного алгоритма, или неявно:где искомая величина х„ і входит одновременно и в левую, и в правую часть. Соответственно численные методы (разностные схемы) делятся на явные и неявные.
Локальная погрешность численного метода па (« + 1 )-м временном шаге определяется выражениемгде х(ґ„_і) - точное решение при / = t„ - і, х„ -1 - приближенное решение, получаемое по формулам (1.3) или (1.4), при условии, что вместо приближенных значений х,„ х,,_ і,..., x„_f/ -1 используются значения, соответствующие точному решению, то есть \(tn), \(tfl. і),..., x(r„_t7 - і).Глобальная погрешность численного метода определяется выражением
При вычислении глобальной погрешности, получается, что локальные погрешности "переносятся" в точку t ,\- и складываются. В работе Д. Холла и Д. Уатта [13] величину Д(/г) также называют полной погрешностью метода и представляют в виде суммы погрешностей:где Q - погрешность округления, обусловленная числом разрядов ЭВМ и зависящая от числа операций, выполняемых для нахождения следующей точки приближенного решения; cw - методическая погрешность, обусловленная погрешностью выбранного алгоритма; ctl - переходная погрешность, обусловленная тем, что при расчете точки x.v, вместо х(/Л . i), x(t\-2),...., x(/.v ц і) берутся приближенные значения x,.v_i,x.v 2-.----. х д._t/ , і полученные па предыдущей итерации.
Одной из центральных проблем в численных методах для обыкновенных дифференциальных уравнений является получение надежных оценок глобальной погрешности. Естественно требовать, чтобы эту погрешность можно было сделать сколь угодно малой, выбирая достаточно малый шаг h. В этом заключается понятие сходимости.
Пусть А{Н) - некоторая функция переменной h с конечной областью определения D на полуоси h 0, причем И может принимать сколь угодно малые значения. Тогда, если существуют такие положительные числа /?0, С, р, что при всех И є D, удовлетворяющих условию 0 h h0, выполняется неравенството пишут Д(/г) - 0(hp) и говорят, что функция A{h) есть О большое от h1 . Число р называется порядком численного метода, если его глобальная погрешность есть О большое от h1 . Можно показать, что если локальная погрешностьимеет порядок {р + 1), то есть А(И) - 0(h J+ ), то глобальная погрешность имеетна единицу меньший порядок, то есть A(h) - 0(hp).
Величина шага интегрирования определяется необходимостью обеспечения устойчивости численного решения системы уравнений (1.2), поскольку при числовой неустойчивости в течение каждого последующего шага интегрирования ошибка интегрирования неуклонно растет. Анализ устойчивости численных методов проводят на скалярном уравнении, которое в работе А. В. Пантелеева [14] называют "тестовым":где X = с +j(a - комплексная константа.Из работы Л.О. Чуа [2] следует, что результаты исследования устойчивости численного метода при решении уравнения (1.8) применимы для линейных и нелинейных СОДУ более высокого порядка. Линейность уравнения (1.8) позволяет получить значимые критерии устойчивости в явной форме. Точное решение уравнения (1.8) имеет вид х(!) = ехр(л/), оно асимптотически устойчиво, если Re X 0, неустойчиво при Re X 0, и устойчиво, если Re X = 0.
Устойчивость численного метода определяется устойчивостью соответствующего численного решения уравнения (1.8). Численное решение уравнения (1.8) асимптотически устойчиво, если оно стремится к пулю при tH - nh — ос.
Устойчивость численного решения зависит в общем случае от h и X. Множество значений ИХ, удовлетворяющих условию асимптотической устойчивости численного решения, называется областью устойчивости метода в комплексной плоскости ИХ. Виды устойчивости численных методов решения СОДУ подробно рассмотрены в работах 10. В. Ракитского, Л. О. Чуа, Н. Ы. Калиткина, Д. Холла и Д. Уатта, А. В. Пантелеева [1, 2, 5, 13, 14].Рассмотрим наиболее известные определения жесткости СОДУ. Аппроксимируем систему (1.2) в некоторой окрестности точного решения x(t) линеаризованной системой:dx где J(t) - матрица Якоб и, составленная из вычисленных в точке (г, х(г)) частных производных 8f /Эх. Если изменение J(/) на некотором интервале г достаточно мало, то локальные фундаментальные решения (1.9) приближенно равны exp(k;t), где "kj = At(t) - локальные собственные значения матрицы Якоби J(() (которые предполагаются различными).Таким образом, решение х СОДУ (1.2) в окрестности точного решения х(/) представляется в видегде С, - постоянные, а ,, - собственные векторы матрицы J(t). Фундаментальное решение схр(?4/) характеризует локальный отклик системы на малые изменения или возмущения \(t).
Производные компонент вектора х(ґ) на интервале t e[t0,T] могут достигать величин порядка L max \x.-(t)\, где L - константа Липшица, определяемаянеравенствами:Рассмотрим два наиболее известных определения жесткой СОДУ.
Исследование методов Розенброка, трапеций и Гира
При решении нелинейных СОДУ неявные методы требуют решения нелинейных алгебраических уравнений. Наиболее распространенным методом решения нелинейных алгебраических уравнений является итерационный метод Ньютона, в западных источниках его называют методом Ньютона-Рафсона. Численное решение чисто неявными методами Рунге-Кутта 2-го порядка и выше требует довольно трудоемкого итерационного процесса. Вследствие указанных выше недостатков чисто неявных методов Рунге-Кутта итерации могут сходиться плохо или даже не сходиться. Поэтому многие вычислители предпочитают использовать безытерационные схемы Розенброка.
На рисунке 2.5 приведеш.! зависимости от шага локальных погрешностей методов Розенброка со следующими параметрами:роста которых, соответственно равны:
Графики локальной погрешности, приведенные на рисунке 2.5, рассчитывались так же, как графики локальной погрешности неявных методов Рунге-Кутта. Как видно из рисунка 2.5 методы Розенброка не всегда обладают L-устойчивостью, а если и обладают, то степень L-устойчивости, меньше порядка метода. Это означає ]-, что методы Розенброка могут иметь те же недостатки, что и метод трапеций, который будет рассмотрен ниже.х0 - [\ О] г, точное решение которой определяется по формуле (2.6) и состоит из двух экспоненциальных составляющих с постоянными времени т і = 1 /jX] = 1 и т2= 1 /1 21= 10 Пилообразная линия соответствует решению, полученному с шагом h — 0,025. Во втором случае шаг был переменным; на первых І 0-ти шагах h ] - 0,0025, на последующих - // 2 = 0,025.
Если в пограничном слое (при г тт) выбрать большой по сравнению с Х2 шаг, то быстрая компонента аналитического решения будет аппроксимироваться в численном решении с большой и весьма слабо уменьшающейся от шага к шагу ошибкой, обусловленной отсутствием /.--устойчивости метода, Именно эта ошибка приводит к пилообразному, «звенящему» характеру кривой 1 на рисунке 2.6. Если на первом интервале выбрать малый шаг к т2, то быстрая компонента решения будет воспроизведена в численном решении достаточно точно и хотя при последующих больших шагах эта ошибка аппроксимации быстрой составляющей не будет уменьшаться, но и влиять на общее решение вследствие малости самой составляющей не сможет. 5! Локальные погрешности методо» Гира, как и ранее для других методов, рассмотрим на примере численного решения тестового уравнения (1.8). Разностная схема неявного метода Гира первого порядка совпадает с разностной схемой неявного метода Эйлера, локальная погрешность которого определялась выше. Локальные погрешности методов Гира 2-го порядка и выше определяются не так, как это делалось ранее, поскольку понятие функции роста справедливо только для одношаговых методов. Рассмотрим подробно методику расчета локальной погрешности для метода Гира 2-го порядка. Разностная схема этого метода при решении уравнения (1.8) принимает видпогрешности заменяют на x(t\) = ехр(/Л), тогда локальная погрешность метода Гира 2-го порядка)
Аналогично можно определить локальную погрешность метода Гира 3-го порядка, которая имеет вид:На рисунке 2.7 приведены зависимости от шага локальных погрешностей методов Гира 1 - 3-го порядков при решении тестового уравнения (1.8).Рисунок 2.7 - Зависимости от шага локальных погрешностей методов Гира 1 -3-го порядков (сверху вниз)
Из рисунка 2.7 следует, что методы Гира /,-устойчивы. Отметим, что зависимости А(/г) j получены при точных стартовых значениях, что позволяетоценить локальную погрешность методов Гира в наилучшем случае. Приведенные на рисунке 2.1 и рисунке 2.7 результаты показывают, что метод Гира сильно уступает по устойчивости чисто неявным методам Рунге-Кутта, особенно при больших порядках точности. Этот недостаток метода Гира связан с тем, что численное решение уравнения (1.8) методами Гира 2-го порядка и выше описывается во времени не одной, а двумя экспоненциальными составляющими, одна из которых стремится к точному решению, а другая - нет, то есть является паразитной. Наличие паразитной составляющей можно наглядно проиллюстриро вать при решении уравнения (1.8) методом Гира с большим шагом интегрирования. На рисунке 2.8 приведены зависимости погрешности от независимой переменной / при решении уравнения (1.8) с постоянным шагом h = 1 для методов Гира 1 и 2-го порядка. Для сравнения на этом же рисунке приведены аналогичные зависимости для неявных метода Рунге-Кутта 2-го порядка и метода трапеций.
Из рисунка 2.8 видно, что численное решение методом Гира 2-го порядка имеет выраженную колебательную составляющую, а для неявных методов Эйлера, трапеций и Рунгс-Кутта этот недостаток отсутствует. Колебательная составляющая численного решения методами Гира 2-го порядка и выше, проявляется и при анализе глобальной погрешности. На рисунке 2.9 приведены зависимости от шага глобальной погрешности A(h) і - х, п - x{(tri)\ методов Гира1-3-го порядков для первой компоненты решения (2.5) при fl}= 1. Для сравнения на этом же рисунке приведены аналогичные зависимости для неявных методов Рунге-Кутта 2-го и 3-го порядков и метода трапеций. Параметры СОДУX 2 = - 102), и х0 = [і О] т, коэффициент жесткости А. \ тах/1 л т[п = 102, а решение определяется по формуле (2.6). Начальные условия (разгонные точки) для методов Гира 2-го и 3-го порядков вычислялись неявными методами Рунге-Кутта соответствующего порядка.Рисунок 2.9 - Зависимости от шага глобальной погрешности: 1 - неявного метода Эйлера; 2 и 3 - методов Гира 2-го и 3-го порядков соответственно; 4- метода трапеций; 5 и 6 - чисто неявных методов Рунге-Кутта 2-го и 3-го порядков соответственно.
Рассмотрим работу известных методов оптимизации на примере задачи идентификации параметров статических моделей полупроводникового диода. Будем определять параметры моделей диодов FR102 и MURD315, ВАХ которых приведены на рисунке 5.2. Сначала рассмотрим модель представленную в виде управляемого собственным током, источника напряжения
Параметры модели определим из условия минимума суммы квадратов относительных погрешностей где u(ij) - рассчитанные с помощью модели значения напряжения; щ - измеренные значения напряжения (/ = 1, 2, ..., N, N - число точек экспериментальной ВАХ).
Исходная оптимизируемая функция S S(IS, р), является двумерной, поскольку зависит от двух параметров. Используя для нее необходимое условие минимума dS/d(.p = Q, зависимость р(Д-),можно определить аналитически и перейти от двумерной функции S(Is, ср) к одномерной S(JS) = S\Js, ф(Д)] Жесткость задачи г\2 и Пз ПРИ минимизации двумерных и трехмерных функций, соответственно, можно оценить по методике приведенной в главе 5. В таблице 2.1. приведена жесткостьзадачи, а также начальные приближения (До и ро) для методов оптимизации, при определении параметров модели (2.12) для диодов FR102 и MURD315. Параметры модели (2.12) и соответствующая им величина S, определенные различными методами и при различной размерности минимизируемой функции S для диодов FR102 и MURD315 приведены в таблице 2.2 и таблице 2.3 соответственно.
Метод глобальной подгонки
В данном подразделе предлагается модифицированный -устойчивый метод экспоненциальной подгонки - метод глобальной подгонки, порядок точ ности которого для глобальной погрешности совпадает с порядком метода трапеций.
При решении сингулярпо-возмущепых задач при фиксированной допустимой погрешности, шаг интегрирования быстро изменяется в ходе решения в широких пределах, причем увеличение шага чередуется с его уменьшением. Рассмотренный выше метод экспоненциальной подгонки не обладает L-устойчивостью и поэтому скорость увеличения его шага определяется не только допустимой точностью, но и скоростью затухания быстрых компонент решения. Его отличие от метода трапеций в этом отношении - только количественное.1-устойчивость численного метода при решении жестких задач проявляется, прежде всего, в сохранении качественного характера решения независимо от шага. Такой же эффект может быть получен, если при больших шагах решения использовать в (3.2) значения п, близкие к 0,5 и близкие к нулю при малых шагах, то есть использовать вместо постоянного значения г функцию rf = F(h). Такая функция может быть получена в предположении, что максимальный шаг интегрирования /гШХ (hmaK Г, обычно hmnx 7750) соизмерим с наибольшей постоянной времени анализируемой цепи тпт, а порядок величины минимального шага hmn равен порядку минимальной постоянной ттіп при очевидном ограничении на отношение х1ШХ! ттш, характеризующем жесткость системы.
В качестве функции r\ - F(h) с целью экономии вычислительных ресурсов автором данной работы [29] предложено использовать рациональнуюпфункцию в виде степенного многочлена Эрмита г (h) = (- (h/hmax) , коэффи циенты а, которого определяются из условий: г] (0) = 0, г (/гтач) = 0,5,а, для п 6 приведены в таблице 3.1 (для всех п а„= 0). Рисунок 3.8 - Граница области устойчивости предлагаемого метода при различных Г).
Область устойчивости предлагаемого метода при h— 0 стремится к области устойчивости метода трапеций, а при h — ІіІШХ— к области устойчивости неявного метода Эйлера, следовательно, метод Л-устойчив.
Рассмотрим следующую жесткую модельную задачу, которую применяли в своей работе В. X. Пресс и С. А. Тенкольский [30] для исследования методов решения жестких СОДУ:с начальными условиями w(0) = 1, v(0) = 0 и с собственными числами I = VTma4 = _1 И 2 =-1Amm =-1 000, Т. Є. lg(Tiliax / Tmin) = 3.
При введении зависимости v\ (h) понятие «асимптотический порядок точности» для локальной погрешности теряет содержательный смысл, поэтому далее будем рассматривать глобальную погрешность метода Д - uk -u(tL), 0 /А Г,далее /А=А1ШХ=ттах.
На рисунке 3.9 приведены зависимости глобальной погрешности Л от шага численного решения (3.9) предлагаемого метода при различных значенияхп, то есть для различных r[ (h). Для сравнения па этом же рисунке приведены зависимости A(h) неявного метода Эйлера, метода трапеций и метода трапецийс пассивным сглаживанием (последний случай соответствует пунктирной линии, остальные линии на рисунке 3.9 соответствуют тем же методам, что и на рисунке 3.1).
Из рисунка 3.9 видно, что при больших шагах метод трапеций, как и следовало ожидать, обладает наибольшей глобальной погрешностью. При малых шагах зависимость Д(7г) для предлагаемого метода и метода трапеций имеет одинаковый наклон (т. с. порядок глобальной погрешности методов одинаков). В случае п lg(TmLlx / тШП) погрешность предлагаемого метода монотонно убывает и, по-видимому, всегда меньше погрешности неявного метода Эйлера.
Решение неявных разностных схем Рупге-Кутта блочными итерациями Зейделя
Методы построенные па основе итераций Зейделя и рассмотренные в предыдущем подразделе - (а)-устойчивы только в том случае, если матрица А близка к диагонально-доминирующей и, наоборот, максимально неустойчивы при нулевых коэффициентах на главной диагонали. Это означает, что непосредственное применение итераций Зейделя целесообразно только для сравнительно узкого класса RL и Я С-цепей.
В работе В.Н. Бирюкова и А. М. Пилипенко [41], автором данной диссертации показано, что и при наличии нулевых элементов на главной диагонали матрицы А возможно существование .4-устойчивого метода, трудоемкость которого пропорциональна размерности задачи. Этот наихудший случаи, соответствует линии передачи без потерь.
Пусть вектор хп и сеточной функции {хп} определен следующим алгоритмом:1. Па первом полушаге определяем значение первого элемента x]lHi5 искомого вектора путем решении неявным методом Эйлера системы из первых двух уравнений исходной СОДУ, принимая х ; 10; - х3п .2. Определяем х 0 5 путем решения системы из 2-го и 3-го уравнений СОДУ, используя найденное ранее значение А- +() 5 и принимая x .,ti $ = х .3. Определяем аналогично пп. 1, 2 последовательно -т +05, х11+05, ...,
Последовательность действий, устанавливаемая изложенным выше алгоритмом, также приводит к симметричным итерациям Зейделя, но не для дифференциальных уравнений первого порядка, как в (4.1), (4.2) а для систем уравнений второго порядка, то есть к блочным итерациям. Очевидно, что СОДУ 2-го порядка могут быть решены и неявными методами Рунге-Кутта высоких порядков. Чисто неявные методы Рунге-Кутта в случае жестких СОДУ порождают плохо обусловленные системы алгебраических уравнений [1], однако увеличение пространственной дискретности приводит к резкому увеличению жесткости только исходпой СОДУ, жесткость решаемых СОДУ 2-го порядка мала и от дискретности не зависит.
На рисунке 4.13, а показана область устойчивости, определенная для приведенного алгоритма. На рисунке 4.13, б показана область устойчивости при замене метода Эйлера неявным методом Рунге-Кутта 2-го порядка. Поскольку предлагаемый метод осуществим для СОДУ имеющих порядок не меньше третьего, его области устойчивости определялась, путем численного решения тестовой СОДУ 3-го порядка, соответствующей звену многозвенной цепи (рисунок 4.9, б). Область устойчивости определялась аналогично тому, как это делалось при определении области устойчивости метода (4.1), (4.2). На рисунке 4.13, а граница метода 1-го порядка точности совпадает с мнимой осью, но при ненулевых элементах на главной диагонали матрицы А устойчивость метода растет и граница области с ростом шага отклоняется от мнимой оси в положительном направлении Re(hX). Отметим, что совпадение сохраняется при увеличении размерности тестовой СОДУ. В любом случае метод оказывается Л-устойчивым, что позволяет использовать его как для анализа диссипа-тивных структур, так и для анализа устойчивости структур, потерями в которых можно пренебречь. Метод 2-го порядка точности Л(а)-устойчив. Для вычисления локальной (на одном шаге Л = !х„ і - x(t „ ,-\)\) и глобальной (на ряде шагов A =j хп -x(tn) ) погрешностей предлагаемого методарассматривалась тестовая СОДУ (2.5) при т = 3. Для того, чтобы решение тестовой СОДУ, соответствующей цепи на рисунке 4.9, б, как и в общепринятомслучае имело вид экспоненты, необходимо соответствующим образом выбратьпараметры цепи и начальные условия. В этом случае
Приведенные на рисунке 4.14, а результаты показывают, что вид локальной погрешности исходного метода при малых шагах и блочных итерациях сохраняется, но при больших шагах погрешность близка к погрешности метода трапеций, то есть Піп Д — 1.
Отметим, что при решении СОДУ, соответствующей диссипативной структуре, когда на главной диагонали пет ненулевых элементов, локальная ошибка при больших шагах оказывается меньше, чем у метода трапеций, так как в этом случае 1ипА 1. Именно из-за сравнительно высокой локальнойошибки наблюдается быстрый рост глобальной ошибки A = j хп -x(tn) [ прибольших шагах (см. рисунок 4.14,6). Для исследования глобальной погрешности предлагаемого метода решалась та же тестовая СОДУ, что и при анализе локальной погрешности. Глобальная погрешность вычислялась для достаточно большого момента времени /„ = 10.
Рассмотрим свойства предлагаемого метода при диагонально-доминирующей матрице А, которой обладает СОДУ, описывающая, например, ЛС-линию (рисунок 4.15). (метод (4.1), (4.2) в этом случае /і(а)-устойчив). Из сравнения зависимостей глобальных погрешностей от шага, определенных для различных методов и показанных на рисунках 4.16 и 4.17 следует что, использование блочных итераций не позволяет получить /--устойчивости и при диагонально-доминирующей матрице А.
Рисунок 4.16 - Зависимость от шага глобальной погрешности метода Эйлера(кривая 1), того же метода при использовании симметричных итераций (2) исимметричных блочных итераций Зейделя (3) при диагонально-доминирующейматрице Л.
Рисунок 4.17 - Зависимость от шага глобальной погрешности метода трапеций (кривая 1), L -устойчивого неявного метода Рунге-Кутта второго порядка (2) и //-устойчивого чисто неявного метода Рунге-Кутта второго порядка при использовании блочных итераций Зейделя (3) при диагонально-доминирующей матрице А.
Предлагаемый метод с блочными итерациями оказывается достаточно эффективным в задачах анализа многозвенных цепей моделирующих линии передачи. Это можно увидеть на следующем примере. Пусть отрезок однородной линии с распределенными параметрами моделируется многозвенной эквивалентной схемой, показанной на рисунке 4.18 (потери в диэлектрике приняты малыми по сравнению с потерями в проводниках). Схема выбрана из условия, что весь отрезок длинной линии разбивается па равные по длине участки, каждому участку соответствует ячейка на схеме, параметры соответствующих элементов каждой ячейки одинаковы: Rt = R2 = ... - R,\, С) = С2 = ... = С,\, L{ = L2 = ... = L\, где .V - число ячеек многозвенной схемы. Очевидно, что чем больше TV, тем точнее будет многозвенная эквивалентная схема описывать свойства отрезка. Требуется найти переходную характеристику git) разомкнутого на конце отрезка такой линии (напряжение на выходе линии при подаче единичного скачка напряжения на вход). Подобная задача возникает при анализе многих устройств, в которых используются отрезки длинных линий.ется условие LkICk»Rk. Было принято, что /?/( = 0,01 Lk IС к . Можно выполнить расчет переходной характеристики для общего случая, не задаваясь конкретными значениями Lk, Q, Rk. Шаг интегрирования при этом определялсякак h= t,tt) , где /о = jLkCk , , - выбираемая при вычислении константа. Расчетпереходной характеристики был выполнен двумя методами: классическим методом трапеций и предложенным методом с блочными итерациями 2-го порядка. Результаты приведены на рисунке 4.19.
