Диссертация (Исследование волноведущих систем методами математической физики), страница 10

Для решения системы уравнений A   B  yk 1  Bxk(3.26)матрица  A   B  в левой части равенства факторизуется.Такимобразом,приреализациивсехрассмотренныхметодовнеобходимо производить факторизацию матриц.3.3.2Факторизация матрицы жесткостиКак уже отмечалось, в результате применения метода конечныхэлементовисходнаядифференциальнаязадачасводитсяк конечно-разностной задаче, то есть к системе линейных алгебраических уравнений(СЛАУ) (3.26), для решения которой матрица должна быть предварительнофакторизована. Причем выбор метода факторизации ограничен спецификой72задачи, которая имеет ряд особенностей.

Входящие в СЛАУ матрицыоказываются разреженными матрицами (очень часто ленточной структуры)весьма высокого порядка, но главное заключается в том, что матрицы,возникающие при решении электродинамических задач, как правило,оказываются незнакоопределенными. Все это затрудняет применять дляобработки таких матриц стандартные методы [92].Рассмотрим более подробно особенности решения системы (3.26).Во-первых,вследовательно, исилуA ,незнакоопределенностистановитсяматрициAB,а,невозможным применение методаХолецкого [93].

До публикации диссертации Джемса Банча [94] - [96] в 1969году (Калифорнийский университет в Беркли) не существовало устойчивогоалгоритма для обработки симметричных незнакоопределенных матриц.Во-вторых, нежелательным является применение методов с выборомглавного элемента по строкам либо по столбцам, ибо при перестановках,очевидно, утрачивается такое важное свойство матрицыAкак ееленточность.В-третьих, не представляется возможным использовать метод Гаусса безвыбора главных элементов, так как на главной диагонали в процессефакторизации могут появиться малые элементы.В нашем случае в силу особенностей матрицы А для ее факторизацииприменяется алгоритм, основанный на стратегии, разработанной Дж.

Банчеми Л. Кауфман.Дадим описание метода, позволяющего факторизовать матрицу A сучетом всех ее особенностей.ПустьсимметричнаяневырожденнаяматрицаA   aij  , i, j  1, nпредставлена в блочном виде: S CT  ,CBA  73где S- невырожденная подматрица размера k  k ; k<n. Тогда справедливотождество: Ek0  A -CS 1 E nk E k0-S 1CT   S 0 ,Enk   0 D1 или иначе Ek 0 S 0  E kA   1 CS E   0 D   01 nk  S 1C TEnk,где учтена симметрия матрицы S, следующая из симметрии матрицы A. ЗдесьEk - единичная подматрица порядка k; S T   S 1  , причем D1  B  CS 1C T .TТаким образом, матрица A может быть представлена в виде произведениятрех матриц:3.27A  L1D1L1Tпричем L1  нижнетреугольная матрица, а D1 - блочно-диагональная матрица, Ek 0S 0 ; L  0 D L E1 1 n kD1   ; L1  CS 1.Далее, производя аналогичную процедуру с матрицей D 1 , т.е.

представляяее в виде D1  L2 D2 LT2 , записываем матрицу A как произведение матрицL1L2 D2 LT2 LT1 , где74L2 , D2  матрицы n  k    n  k  ; L2 размера S1 0 D2    блочно-диагональная матрицы, а матрицы0D3 EkL2   000Ek1L3;En k1 k 00S 0 0 D2   0 S1 0  ;0 0 D 3нижнетреугольная,L2 , D2  размера n  n :S1 - матрица размера k1  k1.Повторяя затем эту процедуру, производим разложение матрицыA иполучаем ее представление в факторизованном виде:A  L1L2...Li1Li Di LTi LTi1...L2T L1T  LDi LT ,где L  нижнетреугольная матрица.Разработанный Банчем и Кауфман метод предполагает в качестве S k наk-м шаге выбирать матрицу 11 (скалярный шаг исключения) или 2  2(блочный шаг) в зависимости от следующего условия.Пусть элемент d11   max di1 , где   достаточно малое вещественное числоi 1 0    1 .

Тогда вычисляемL и D для k  1. В противном случае расчетыпроводятся при k  2 .Рассмотрим теперь непосредственно вычисления, производимые прискалярном и блочном шаге.1. k  1. Вначале определяются значения элементов D :dij  bij  cik c jk / s11; i, j  1, n 1,А затем пересчитывается столбец C длины n  1 :li  ci / s11; i  1, n 1; li   L2  CS 1.k  2. Рассмотрим T  S 1. Пусть  det S; тогда:75T1  S22 - S12  ;   S11S22  S122 ( матрица симметрична).

Элементы матрицы  -S12 S11 D вычисляются по формулам:dij  ai 2, j 2  li1c j1  li 2c j 2  ai 2, j 2  ci1l j1  ci 2l j 2 ;(3.28) i, j  1, n  2; li, j  L  CL1При этом в виду симметрии D достаточно вычислить лишь элементыdij , i  j.Отметим, что единовременное вычисление всей матрицы D требуетдополнительного массива размера  n  2   2, т.к. для определения элементовd ij по формулам (3.28) необходима и матрица C .Для решения этой проблемы производится вычисление элементов постолбцам. Именно, фиксируем j  1 и определяем элементы di1 согласно (3.28)дляi  1, n  2.Сначалавычисляемl11  c11t12  c12t21; l12  c11t12  c12t22 ,элементыl11(tij )  L1 , присваиваяиl12 :их значениярабочим переменным.

Когда все элементы первого столбца D будутнайдены, элементы c11 и с12 не требуются для дальнейших вычислений, и наих место можно записатьДалее определяем :l11,l12 .l21  c21t11  c22t21;l22  c21t12  c22t22 ,находимпо(3.28)приj2элементы di 2 , записываем элементы l21 , l22 на место элементов c21 и c22 и такдалее.Применяя теперь описанный алгоритм к матрицам D1 , D2 и так далее,производим разложение матрицы A до тех пор, пока получаемая матрица Diне будет иметь порядок 11 или 2  2.В заключение следует отметить, что на практике такая стратегияоправдывает себя, так как, даже если dkk  0 , то определитель матрицы S k ,соответствующей блочному k-му шагу исключения, оказывается достаточно76велик и чрезмерного роста элементов, приводящего к потере устойчивости,не происходит.Рассмотрим вопрос об организации хранения матриц A и B .

Обе ониявляются ленточными разреженными, что наводит на мысль о хранении их впамяти ЭВМв строчном формате одним из стандартных методов [59].Однако в нашем случае компактное хранение оказывается затруднительным,поскольку при факторизации матрицы A происходит заполнение внутриленты матрицы. Кроме того, возможно уширение ленты, приводящее кувеличению ширины ленты на единицу. Оба эти явления носят нескольконепредсказуемый характер, так что при взгляде на исходную матрицуопределить, какие именно нулевые элементы на некотором шаге обратятся вненулевые, практически невозможно.Именно поэтому в данной работе использовалась довольно простая,однако как нельзя лучше подходящая к применяемому методу факторизациисхема хранения: в виде прямоугольного массива хранится нижняя половиналенты (полулента), что допустимо в силу симметрии матриц.Рассмотрим некоторые подробности, касающиеся реализации этой схемы.Пусть нам предстоит применить ее к матрице A размера N  N ,полуширина ленты которой Sh  N1.

Учитывая, что факторизация, возможно,увеличит Sh на единицу, создадим прямоугольный массив C размераN   N1  1 , и поместим в него полуленту исходной матрицы A   aij  ; i, j  1, n,перенося элементы так, как показано на рисунке 3.4.

Крестиками на рисункеобозначены ненулевые элементы матрицы A и массива B . В данном примереN=8, N1=3. Несложно заметить, что при такой форме записи строкипереходят в строки, столбцы переходят в диагонали, параллельные главнойдиагонали, а диагонали переходят в столбцы.Таким образом, при записи матрицы A в массив B пересчет индексовпроизводится следующим образом: пусть элемент aij  A записывается, какэлементbgh  B;тогда i, j и g, hсвязанысоотношениями77g  i; h  j  N11 i. Обратный переход от  g , h  к  i, j  осуществляется поформулам i  g , j  g  h   N11.X X X 0 0 0 0 0 0 0 0 X X X X X 0 0 0 0 0 0 X X X X X X X 0 0 00 X X X 0 X X X X X 0 00 X X XAB0 0 X X X X X 0 0 X X X 0 0 0 X X X X X0 X X X 0 0 0 0 X X X X 0 X X X 0 0 0 0 0 X X X 0 X X XРис. 3.4.

Организация схемы хранения.3.4Численная апробация: металло-диэлектрический волновод скусочно-постоянным заполнениемНа основании предложенного алгоритма была построена программарасчётаволноводовпрямоугольногосечениясдиэлектрическимзаполнением.Программавключаетвсебятриосновныхмодуляишестьвспомогательных блоков.Блок 1 содержит параметры сетки и основные характеристики волновода,необходимые для проведения расчётов.Блок 2 и блок 3 предназначены для хранения элементных матриц.В блоки 4 и 5 записываются собранные из матриц, хранящихся в блоках2 и 3, матрицы A и B.В блок6 помещаются вычисленные собственные значения задачи, атакже постоянные распространения78Модуль 1 имеет сервисное значение.

Он позволяет редактировать данныеиз блока 1 и вносит изменения в конфигурацию исследуемых волноводов, атакже строит элементные матрицы, записывая их в блоки 2и 3.В модуле 2 осуществляется алгоритм сборки матриц.В модуле 3 происходит реализация алгоритма метода обратных итераций истратегии Банча-Кауфман факторизации матрицы A  A   BТестирование программы происходило на волноводе без заполнения,когдаизвестно(точ )2  k 2   2 ,точноеатакжезначениенапостояннойволноводахсраспространениякусочно-постояннымдиэлектрическим заполнением .Результаты тестирования для пустого волновода приведены на рисунке 3.5и в таблице 1.Рисунок 3.5.