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

1 0, 0  hx , 0 Рис. 3.1. Построение сетки путем триангуляции сечения волновода:элементы первого и второго типа ( TI и TII ).Будемиспользоватьлинейнуюполученных треугольников, то естьинтерполяциюнакаждомизприменять линейные конечныеэлементы [52].61Рассмотрим вначале элемент типа 1 , обозначив его через T1 . ОбозначимузлыT1черезкакi, m, n,показанонарисунке.Тогдаимеем:и строим линейные базисные функции:X n  Yn  X m  Yi  0; X i  hx ; Ym  hyi , m , n : i  x / hx ; m  y / hy ; n  1  x / hx  y / hy .Каждую из компонент вектора магнитного поля H в элементе T1 теперьаппроксимируем следующим образом:       H  ;    H  ;    H  .H I   H   Hxmxi xiI H y  i H y   m H yiH I   H   Hzmzi zi где H qmmmnxnynzn3.7 nn значение H q в   ì узле,   i, m, n; q  x, y, z.Вернемся теперь к уравнению  3.4a  .

Рассмотрим его на элементе T1 ,H I  H xI , H yI , H zI . Взаменяя вектор H аппроксимирующем его векторомрезультате получим матричное уравнение:3.8AH I   2 BH IС 1 Т I  дифференцируемых на элементе TI ,Рассмотрим пространствовообщеговоря,комплексныхфункций,введемнанемскалярноепроизведение (звездочка обозначает комплексное сопряжение):u, v    u  x, y  v  x, y  dxdy,u, v  C   T1 1(3.9)TIи нормуu   u( x, y)u  ( x, y)dxdy.2(3.10)TI62Очевидно,i , m , n ,функцииа,следовательно,ифункция H q I ; q  x, y, z принадлежат пространству С  TI  .1Рассмотримвектор-столбцы Ψ,составленные из базисныхi Ψm , Ψn ,функций: Ψi  i , i , i  , Ψm  m , m , m  , Ψn  n , n , n  .TTTУмножим скалярно левую и правую части уравнения (3.8) последовательнона вектор-столбцы Ψi , Ψm , Ψn .

В результате получим систему уравнений: AH , Ψ     BH , Ψ  ,  = i, m, n.2IIВыражая в уравнении (3.11) H I ,(3.11)черезнабор базисных функций  i, m, n, перепишем уравнение (3.11) в виде:AI H I   2 BI H I ,(3.12)где                  HI  Hx ; H y ; Hz ; HxiH qi,   i, m, n; q  x, y, zi-Tm; Hymузловые; Hzm; Hx ; H y ; Hznзначенияnкомпонентnвекторамагнитного поля, а матрицы AI и BI имеют блочную структуру: Aii I Aim I Ain I  Bii I Bim I Bin IAI   Ami I Amm I Amn I  ; BI   Bmi I Bmm I Bmn I I I Ani Anm I Ann I  Bni Bnm I Bnn I(3.13)63Блоками матриц AI и BI являются матрицы:A I   A ,   ; B I   B ,   ;  ,   i, m, n.Матрицы, определенные формулами(3.13), связанные с элементом T1будем называть элементными.Проведем теперь аналогичные построения для T2  элемента типа 2.Именно, обозначив его узлы через r , s, t ( см. рис 3.1), строим вначалебазисные функции: r  1 x / hx , s  1 y / hy , t  1 x / hx  y / hy .Каждую из компонент магнитного поля H в TII аппроксимируем аналогично(3.7): H q  r ,s ,t H   ; q  x, y, z.qТогда уравнение (3.4.а) на элементе TIIимеет видAH II   2 BH II ,где H II вектор-столбец : H II  H x , H y , H z(3.14).TУмножая скалярно обе части уравнения (3.6) на вектор-столбцыΨr , Ψs , Ψt , : Ψr  r , r , r  , Ψs   s ,  s ,  s  , Ψt  t , t , t  ,TTTперепишем это уравнение в виде системы: AH II , Ψ    2  BH II , Ψ  ,   r, s, t.(3.15)Соотношение (3.15) может быть, как и для элемента TI , записано спомощью матриц AII , BII  элементных матриц для базисных функцийr ,  s , t элемента типа 2 - следующим образом:64AII H II   2 BII H II ,(3.16)где:                  H II  H x ; H y ; H z ; H x ; H y ; H z ; H x ; H y ; H zrrrsssttTtОтметим, что, поскольку базисные функции  и   ,   i, m, n;   r, s, t,имеют один и тот же вид для всех элементов типа 1 и типа 2 соответственно,то и элементные матрицы AI , BI и AII , BII будут одинаковыми для всехэлементов соответствующего типа, что существенно упрощает процесссборки матриц.Рассмотрим всю совокупность элементов, на которые разбито исходноесечение S..Пусть протяженность сеченияS вдоль оси x Lx  hx  N x ; вдоль осиy : Ly  hy  N y ; тогда S разбито на N e элементов, N e = 2 N x  N y .Занумеруем их так, как показано на рисунке 3.2 (в данном случаеN x  4; N y  2; Ne  16).3691215(4)(3)(7)(2)(8)(11)(6)(12)(15) (16)(10)(5)(9)(14)(13)(1)1471013Рис .

3.2. Нумерация элементов и узлов сетки.65В дальнейшем номера элементов будем заключать в круглые скобки, стем, чтобы не путать их с номерами узлов. Аналогичным образоми узлы нашей сетки (число узлов Nu   N x  1   N y  1 ).пронумеруемНапомним, что каждому узлу ставится в соответствие три величиныH x , H y , H z - компоненты вектора магнитного поля.До сих пор мы рассматривали каждый элемент ( T 1 или T2 ) отдельно, внесвязи с остальной частью области сеченияS. Теперь переформулируемуравнения (3.11) и (3.15) таким образом, чтобы полученная задача включалав себя все элементы.Для этого рассмотрим вначале произвольный элемент T 1 с номером k . Емусоответствуют узлы с номерами N1, N1 1, N1  N y 1. Уравнение (3.11) будетиметь для него следующий вид:A1H k   2 B1H kHk  H xN1  N y 1 ; HyN1  N y 1 ; HzN1  N y 1 ; HxN1 1(3.17) ; HyN1 1 ; HzN1 1     ; HxN1; HyN1; HzTN1Составим вектор       H u  H x ; H y ; H z ; H x ;...; H z1112Nu 1     ; HxNu; HyNu; HzTNuи перепишем уравнение (3.15) в следующем виде:Ak Hu   2 Bk Hu ,где(3.18)блочные матрицы Ak и Bk состоят из блоков A1 , и B1 , ,  ,   n, m, iэлементных матриц, стоящих на пересечении строк и столбцов с номерами66N1 , N1  1, N1  N y  1.

Вид матриц0 0 0=  0 0 0  0 0 0нулеваяматрицаAkи Bk показан на рисунке 3.3. Здесьразмером3  3.Проводяаналогичныерассуждения для всех элементов типа T I , приходим к N x  N y уравнениямтипа (3.12) для k  1,3,..., 2 N x  N y 1.Для элементов типа T II проводимподобную же процедуру, исходя из уравнения (3.16), и получаем сноваN x  N y уравнений для элементов с номерами k  2, 4,..., 2 N x  N y .Суммируя теперь полученные уравнения по всем элементам, получаем: 2 Nx N y  2 Nx N y 2  Ak  H u     Bk  H u , k 1 k 1после чего, обозначая: 2 Nx N yAk   A; k 1 2Nx N yBk   B, k 1перепишем это уравнение в виде:AHu   2 BHuПолученная обобщенная задача на собственные значения(3.19)позволяетопределить постоянные распространения и соответствующие им модовыерешения уравнения (3.3а) в каждом из узлов сетки.Отметим несколько особенностей матриц A и B , фигурирующих вуравнении (3.19).Во-первых, из вида матриц A и B (3.5а, б) следует, что в случае кусочнопостоянных коэффициентов матрицы A и B являются симметричными,однако не являются знакоопределенными.67 . . . ..Ak   . . . ....

... .... .... .... A nn1... A mn1... .... .... .... A in1... .... .... ... ...  .... . ... . .... . ... . .... . ... . A nm1 ... A ni1 ...  A mm1 ... A mi1 ...  . ... .... . ,. ... .... . . ... .... . A im1 ... A ii1 ...  . ... .... . . ... .... .

 ...  ...   . . . ..Bk =  . . . . . .... .... .... ......... Bnn1Bnm1 ... B ni1 ...... Bmn1Bmm1 ... Bmi1 .................................... .... ...... .... ............ ..... ........ Bin1... .... .... .... Bim1 ..................Bii1 ... ....

..... ..... . ... Рис. 3.3. Вид блочных матриц Ak и Bk , содержащих элементные матрицы.Во-вторых, из формулы (3.5а) следует, что матрица A - вырожденная.Для того чтобы избежать связанных с численной реализацией трудностей,являющихсяследствиемэтогообстоятельства,переформулируемобобщенную проблему собственных значений (3.19) следующим образом.Прибавим к обеим частям равенства (3.19) выражение k 2  BHu :   max  : A  k  B H2u  k 2    2 BHu .Вводя обозначения A  A  k 2  B, B  B,   k 2    2 , получим:AHu   BHu .(3.20)Матрицы A и B являются ленточными матрицами. Ширина ленты зависитот разбиения области на элементы и способа нумерации узлов и вычисляетсяпо формуле Sh  Q  R 1 , где R- максимальная по элементам величина68наибольшей разности между номерами узлов в отдельном элементе, а Q число степеней свободы (неизвестных) в каждом узле.

Для приведенногоразбиения прямоугольного сечения S принятая нумерация узлов являетсяоптимальной с точки зрения минимизацииширины ленты, котораяоказывается равнойSh  3( N y  2)(3.21)3.3 Процедура Банча-Кауфман3.3.1 Обратные итерацииДлянахождениясобственныхвекторовисобственныхзначенийобобщенной задачи на собственные значения применяется метод обратныхитераций [41].Рассмотрим стандартную алгебраическую проблему собственных значений:Ax   x(3.22)Метод обратных итераций по существу представляет собой примененный кматрицеA1степенной метод, суть которого состоит в следующем.Выбирается произвольный ненулевой вектор-столбец e1 и строится алгоритм:y  e ,1 0 xk  yk / yk yk 1  Axk ., 2 n Здесь введена евклидова норма n-мерного вектора x   x  , x ,..., x  :1691/22 nix    x    i 0  .На  m  1 -м шаге получается, что Ax / Ax  x  x и, следовательно,величина h  Ax  Ame1 / Am1e1 является приближением к собственномузначению матрицы A , а вектор x - приближением к соответствующемусобственному вектору.Таким образом, описываемый процесс сходится к собственному вектору,соответствующему доминирующему собственному значению, и позволяетопределить само это собственное значение.Однако, степенной метод, к сожалению, зачастую не в состоянииобеспечить достаточную точность вычислений в силу значительногонакопления погрешности.

Кроме того, алгоритмы, основанные на степенномметоде, имеют, как правило, не слишком высокую скорость сходимости.Более эффективным, поэтому оказывается применение метода обратныхитераций.Предполагая, что А – симметричная невырожденная матрица, помножим(3.22) слева на A1 ; тогда получим:A1x 1xОтсюда следует, что собственные значения матрицы(3.23)A1являютсявеличинами, обратными к собственным значениям матрицы А, а собственныевекторы этих матриц совпадают.Применяя теперь к задаче (3.23) степенной метод, получим алгоритм70y  e ,1 0 xk  yk / yk ,1 yk 1  A xk ,где e1 – начальный вектор – выбирается произвольным образом.При реализации алгоритма явное умножение обратной матрицы A1 навектор xk привело бы к возникновению таких же трудностей, что и прииспользовании степенного метода.

Поэтому на практике, как правило,каждый вектор yk 1получают в результате решения системы линейныхуравнений Ayk 1  xk . Для решения этой системы матрица А предварительнофакторизуется.Для решения обобщенной проблемы собственных значенийAx   Bx(3.24)алгоритм метода обратных итераций записывается следующим образом: y0  e1, xk  yk / yk , Ayk 1  Bxk .Методобратныхитерацийсходится(3.25)ксобственномувектору,соответствующему ближайшему к нулю собственному значению задачи(3.24). А именно, если 1  2  i1 , i  2, то обратные итерации сходятся какгеометрическая прогрессия со знаменателем c0  1 / 2 .Длявычисленияостальныхсобственныхзначенийэффективнымоказывается введение сдвига  . Рассмотрим алгоритм (3.25), где вместоматрицы A фигурирует  A   B  .71Данный итерационный процесс сходится к собственному вектору,соответствующему собственному значению задачи (3.24), наиболее близкомук  .

Обозначая это собственное значение через s , получаем сходимость соскоростью геометрической прогрессии, знаменатель которой равен:c  s   / min i  .isПри удачном выборе параметра  скорость сходимости будет оченьвысокой. Это свойство выгодно отличает обратные итерации от прямыхитераций (к числу которых относится и степенной метод), сходимостькоторых не может быть существенно улучшена за счет выбора параметра  .Таким образом, реализуяметод обратных итераций, мыполучаемвозможность определить минимальное собственное значение задачи (3.24),наиболее близкое к значению параметра  ..