Методы вычислений матриц переноса упругих деформаций Текст научной статьи по специальности «Математика»

2013 Механика № 3

УДК 539.3+534.2+512.64

Ю.Н. Беляев

Сыктывкарский государственный университет, Сыктывкар, Россия

МЕТОДЫ ВЫЧИСЛЕНИЙ МАТРИЦ ПЕРЕНОСА УПРУГИХ ДЕФОРМАЦИЙ

Дан обзор матричных методов описания распространения волн в слоистых средах. Развивается метод представления матрицы переноса (характеристической матрицы) в виде матричного решения системы обыкновенных дифференциальных уравнений первого порядка. Эта система уравнений называется определяющей. Метод получения определяющей системы уравнений показан на примере термоупругих волн. Рассмотрены традиционные методы нахождения матричной экспоненты: разложение в ряд Тейлора, рациональные аппроксимации Чебышева и Паде, метод масштабирования, методы численного интегрирования дифференциальных уравнений, методы преобразования матриц (метод собственных векторов, QR-алгоритм, жорданова каноническая форма, преобразование Шура, приведение матрицы к блочно-диагональной форме), формулы Лагранжа-Сильвестра, Бэкера, Ньютона, преобразование Лапласа. Представлен метод симметрических многочленов. Симметрические многочлены n-го порядка, введенные автором, использованы для выражения целых функций матриц, в том числе матричной экспоненты. Этот метод не требует вычисления или оценки собственных значений матрицы. Алгоритм вычисления целых степеней матриц, основанный на использовании симметрических многочленов, является наименее затратным по числу элементарных умножений и, следовательно, наиболее точным в сравнении с другими известными методами. Представлены формулы, аналитически выражающие матрицы переноса упругих деформаций второго и четвертого порядка через элементарные симметрические многочлены определяющей матрицы. Дана аналитическая оценка величин модулей симметрических многочленов. Метод симметрических многочленов позволяет контролировать ошибки округления и усечения при вычислении матрицы переноса. Выполнена оценка масштабирующего коэффициента, обеспечивающего надежное вычисление матричной экспоненты с допустимой погрешностью. Вычисление матрицы переноса упругих волн в слоистых средах методом симметрических многочленов имеет преимущества в сравнении с другими подходами по сочетанию таких параметров, как общность, надежность, стабильность, точность, простота, легкость использования и эффективность численного алгоритма.

Ключевые слова: упругие волны, слоистые среды, матрица, экспонента, многочлены, ошибка усечения, масштабирование.

Yu.N. Belyayev

Syktyvkar State University, Syktyvkar, Russian Federation

METHODS FOR COMPUTING A TRANSFER MATRICES OF ELASTIC DEFORMATIONS

Review of matrix methods, which describe the propagation of waves in layered media, is given. The method of representation transfer matrix (or characteristic matrix) as a solution of system of first order differential equations is developed. This system of equations is called the defining. The method of

obtaining the defining system of equations is shown by the example of thermoelastic waves. Traditional methods of finding the matrix exponential are considered: Taylor series expansion, Chebyshev and Pade rational approximations, scaling and squaring, methods based on numerical integration of ordinary differential equation, matrix decomposition methods (eigenvectors method, QR algorithm, Jordan canonical form, Schur decomposition, block diagonal method), Lagrange-Sylvester formula, Baker formula, Newton formula, Laplace transforms. The method of symmetric polynomials is presented. Symmetric polynomial of n-th order, introduced by the author, used to express the entire functions of matrices, including matrix exponential. This method does not require the computation (or estimation) of matrix eigenvalues. Algorithm to compute integer powers of matrices, based on the use of symmetric polynomials, is the least expensive in the number of elementary multiplications and, therefore, the most accurate in comparison with other known methods. The formulas, which analytically express transfer matrix of elastic deformation of 2-4-order through elementary symmetric polynomials, are presented. Analytical estimation of the modulus of symmetric polynomials is given. Symmetric polynomial method allows to control roundoff and truncation errors in the calculation of the transfer matrix. Evaluation of the scaling factor, which provides a reliable calculation of the matrix exponential with admissible error is made. The calculation of the transfer matrix of elastic waves in layered media by symmetric polynomials have advantages compared with other approaches on a combination of parameters such as generality, reliability, stability, accuracy, simplicity, ease of use and efficiency of the numerical algorithm.

Keywords: elastic waves, layered media, matrix, exponential, polynomials, roundoff error, scaling.

Введение

Среда, свойства которой постоянны на каждой плоскости, перпендикулярной к фиксированному направлению, называется плоскослоистой средой. Если выбрать это специальное направление за ось z декартовой системы координат, то плотность среды, коэффициенты Ламе, диэлектрическая проницаемость и другие параметры среды, характеризующие упругие и электромагнитные свойства вещества, будут функциями одной этой координаты z. Естественными средами, имеющими хорошо выраженную горизонтальную стратификацию, являются океан и атмосфера. Этим объясняется, в частности, сверхдальнее распространение звука под водой и радиоволн в атмосфере [1-5]. Во многих случаях земная кора также успешно может быть аппроксимирована слоистой средой, и это используется в теории распространения сейсмических волн [6, 7].

Широким спектром физических свойств обладают слоистые среды искусственного происхождения. Класс таких объектов постоянно расширяется, и области их применения уже в настоящее время исключительно велики.

Значительную роль в разнообразных технических и физических приложениях играют слоисто-периодические структуры. Исследования распространения волн различной физической природы в средах, свойства которых изменяются послойно периодическим образом, весьма обширны по тематике и охватывают большой круг важных вопросов.

В механике деформируемого твердого тела интерес к этой проблеме [8, 9] связан в первую очередь с широким внедрением в современную технику слоистых композиционных материалов, в том числе нанокри-сталлических [10-12].

Постоянно расширяющийся класс слоистых сред и многообразие их возможных применений поддерживают неослабевающий интерес исследователей к проблемам взаимодействия волн с такими структурами. Развитию теории распространения волн в слоистых средах способствуют два равных по силе обстоятельства. Во-первых, дифракция излучений различных диапазонов длин волн лежит в основе методов неразрушающего исследования строения и состава слоев, контроля качества изготавливаемых объектов. Во-вторых, исследования закономерностей распространения волн в слоистых объектах обосновывают создание на их основе новых приборов, устройств и технических конструкций и расширение областей применения существующих.

Стандартный подход математической физики к расчету распространения волн в слоистых средах (см., например, [2, 7, 13-15]) состоит в выводе соответствующего волнового уравнения, в решении этого уравнения в каждом из слоев рассматриваемой среды с последующим «сшиванием» полученных решений на границах слоев. Модель одномерного изменения параметров, хорошо аппроксимирующая многие слоистые среды, а также линейность используемых волновых уравнений обусловили широкое использование матричных методов в теории дифракции волн в слоистых средах. Одна из важных особенностей матричного подхода состоит в том, что он позволяет существенно упростить процедуру «сшивания» решений и, как следствие, значительно облегчить процесс моделирования спектров отражения от многослойных систем.

1. Метод характеристической матрицы и его аналоги

Суть метода характеристической матрицы (МХМ) состоит в следующем. Компоненты волнового поля, значения которых обуславливаются на границах слоя, составляют вектор-функцию ¥(2). Элементами у ^ 2) вектор-функции ¥(2) в зависимости от решаемой задачи

могут быть компоненты тензора напряжений, векторы смещения, составляющие электромагнитного и теплового полей. Значения этой вектор-функции ¥ (2) на противоположных границах слоя связаны друг

с другом матрицей М, элементы которой определяют величины коэффициентов отражения и пропускания волны слоистой средой, то есть характеризуют рассеивающие свойства последней. Поэтому данную матрицу М называют характеристической матрицей (ХМ) слоя.

Пусть пластина занимает область пространства, ограниченного плоскостями 2 = 0 и 2 = d. Тогда

¥ (0) = М ¥ ^). (1)

Если пластина состоит из N слоев, первый из которых лежит в области 0 < 2 < 21, второй - 21 < 2 < 22, ^й слой - 2N-1 < 2 < 2 N = d,

а слои характеризуются матрицами М1, М2,..., MN, то ХМ всей пластины равна произведению ХМ отдельных слоев, образующих пластину,

М = М1 (21)М2 (22 - 21)- ■ ■MN (22 - 21). (2)

Свойство (2) предопределяет наиболее значительные преимущества МХМ. Так как парциальные характеристические матрицы Му,

у = 1,...,N, зависят каждая только от структуры своего слоя, то изменение толщины и/или физических параметров у -го слоя оказывает влияние только на у-ю ХМ Му, а произведения Пё=Мё и П^++Мё остаются неизменными. Поэтому, если при расчете коэффициентов отражения и пропускания волны многослойной средой необходимо

учесть изменение структуры у-го слоя, частные произведения П] =11Мё

и П^у+1Мё вследствие ассоциативности произведения (2) могут

быть вычислены отдельно, и нет необходимости полного перерасчета произведения (2).

ХМ второго порядка была впервые предложена в 1950 г. Абеле

[16] для описания рассеяния света диэлектрическим слоем (без гиро-тропии). Данный подход, изложение которого имеется также в книге

[17] и обзорах [18-20], стал основным методом при расчете оптических свойств многослойных диэлектрических пленок. ХМ в этом случае еще называют интерференционной матрицей [21, с. 440], а сам метод - методом последовательного наращивания слоев [22].

Второй порядок ХМ достаточен также для описания распространения горизонтальной волны сдвига (так называемая БЫ волна [7]) и упругих волн в средах, в которых модуль сдвига ц = 0 [23]. Развитие МХМ применительно к двухволновой рентгеновской дифракции в слоистых кристаллах было выполнено как в геометрии Брэгга [24], так и Лауэ [25].

В исследованиях электронных состояний полупроводниковых сверхрешеток (см., например, [26-28]) очень плодотворным оказался метод матрицы переноса Т. Эта матрица является обратной по отношению к характеристической матрице М: Т = М- . Еще одно название, используемое для МХМ как в электронике, так и в радиотехнике, - это метод матрицы передачи [29, 30].

Наконец, отметим, что метод, заключающийся в получении характеристической матрицы второго порядка, связывающей выходные и входные параметры единичного элемента дискретной периодической структуры, успешно используется при разнообразных технических расчетах. Ссылки на соответствующие работы имеются в обзоре [31].

Матрицы четвертого и более высоких порядков. Характеристическая матрица четвертого порядка была введена для расчета распространения волн упругих напряжений в слоисто-неоднородных средах с плоскопараллельными границами в работах Томсона [32] и Хаскела [33]. Описание данного метода и развитие его в подход, использующий матрицы пятого и шестого порядка, дано в монографии [18].

Развитие МХМ до матриц четвертого порядка в оптике слоистых анизотропных сред было сделано Тейтлером и Хенвисом [34] и применено Берреманом к исследованию жидких кристаллов [35, 36].

Матрицы четвертого порядка были использованы, в частности, для расчета перерассеяния электромагнитных волн слоистыми анизотропными средами в оптическом [37-40] и СВЧ [41] диапазонах, для рассмотрения рентгеновской и мессбауэровской дифракции в условиях полного отражения [42, 43].

При дифракции излучения рентгеновского диапазона частот в кристаллах возможно возникновение трех и большего числа сильных волн, а при рассеянии электронов на тонких кристаллических пленках число дифракционных пучков, появляющихся одновременно, может достигать нескольких сотен. Соответственно, и-волновой дифракции соответствует матрица переноса и-го порядка, которая в теории электронографии носит название матрицы рассеяния [44].

Определяющее уравнение матрицы переноса. Теоретические исследования распространения волн, например упругих деформаций, основывается на решении соответствующих волновых уравнений. Обычно из исходной системы дифференциальных уравнений в частных производных первого порядка относительно компонент поля (компоненты тензора напряжений, вектора смещений, напряженностей электрического и магнитного полей, температуры и др.) выводится уравнение второго или более высокого порядка относительно одной из компонент у у поля или соответствующих потенциалов [2, 10, 15, 17].

Из решений полученных уравнений для каждого слоя структуры составляются элементы матрицы М = |тогу |, удовлетворяющей условию (1).

Аналитические решения в таком подходе были получены для характеристической матрицы второго порядка для случаев однородных по составу слоев [2, 16, 17], для слоев с линейным [21] и экспоненциальным [45] законами изменения параметров среды по толщине слоя. Кроме перечисленных существует еще небольшое число случаев (см., например, [2, 3, 15, 46, 47]), в которых соответствующие волновые уравнения имеют аналитические решения и, соответственно, может быть получен явный вид ХМ. Для нахождения ХМ второго порядка слоистой среды с непрерывно изменяющимися параметрами по толщине в работе [23] использован метод, основанный на решении соответствующих интегральных уравнений Вольтерра методом последовательных приближений. Аппроксимация неоднородного слоя набором слоев с линейным изменением упругих свойств была предложена в статье [48]. Основным теоретическим подходом к расчету ХМ непрерывно слоистой среды является аппроксимация последней наборами однородных слоев [2, 15, 17, 49].

Альтернативный подход к расчету характеристической матрицы слоя состоит в решении системы дифференциальных уравнений первого порядка

d ¥

^г = К¥, (3)

d2

которая может быть получена из исходных уравнений, например с помощью Фурье-преобразования компонент поля по времени и одной (или двух) координат. Порядок и системы (3) равен числу независимых компонент у 1 вектор-функции ¥ (2) и, соответственно, порядку матрицы Ж.

Решение ¥ ^) уравнения (3) выражается через начальные значения ¥ (0) с помощью матрицы Т:

¥ (d) = Т ¥ (0). (4)

В теории дифференциальных уравнений матрица Т, выражающая значение неизвестной вектор-функции ¥(d) через начальные данные ¥(0), называется фундаментальной. В теории дифракции матрица Т, «транслирующая» вектор ¥ (2) на расстояние d, известна как матрица переноса. Значения элементов этой матрицы определяются свойствами слоя рассматриваемой среды, то есть матрицей коэффициентов Ж системы дифференциальных уравнений (3). Поэтому матрицу Ж, как и систему (3), будем называть определяющей.

Из сравнения соотношений (1) и (4) очевидно, что матрица, обратная Т, является характеристической:

М = Т _1. (5)

Наиболее простой связью, определяющей матрицы Ж с матрицей переноса Т и, следовательно, с характеристической матрицей М, является их связь для однородного слоя. В этом случае [51, с. 115]

Т = ехр(Жг). (6)

Один из способов решения системы уравнений (3) для неоднородного по толщине слоя состоит в аппроксимации последнего достаточно большим числом однородных слоев и использовании для матриц переноса Ту каждого из этих слоев соотношения (6), так что

¥ ( ^ ^ = еХр [Ж; ( ^ +1 - ^ )]¥ ( ^ ),

где (2у+1 - 2у ) = dу - толщины слоев.

Существуют десятки методов численного расчета матричной экспоненты. Их обзор представлен в следующем разделе.

2. Методы вычислений матричной экспоненты

Эффективность некоторых алгоритмов вычислений экспоненты матрицы сильно зависит от вида матрицы. Существенную роль в некоторых подходах играют собственные значения и собственные векторы матрицы, даже если они не вовлечены в конкретный алгоритм. Напом-

ним их определения, а также некоторые другие понятия, обозначения и теоремы, связанные с рассматриваемым вопросом.

Некоторые сведения из теории матриц. Для обозначения элементов матрицы (таблицы чисел), скажем А, используются соответствующие строчные буквы ауу с двумя индексами. Сама матрица изображается в виде таблицы, заключенной в прямые двойные скобки и обозначается символами ||агу|| = А . Первый индекс элемента матрицы

указывает номер строки, а второй - номер столбца, на пересечении которых в таблице находится данное число. Элементы аи квадратной

матрицы называются диагональными. Сумма диагональных элементов -это след матрицы. Порядок квадратной матрицы и - это число ее строк (столбцов). С алгеброй матриц можно познакомиться по любому учебнику, например [52]. Определитель матрицы А обозначается ёе!А. Матрица А - невырожденная, если ёе! А Ф 0. Квадратная матрица I, у которой диагональные элементы равны единице, а остальные нули, называется единичной. Матрица А-1, удовлетворяющая условию АА_1 = А_1 А = I, называется обратной по отношению к матрице А. Если поменять в матрице А = ||агу|| строки и столбцы местами, то в резуль-

Т I I л *

тате получается транспонированная матрица А = а. . Матрица А

называется сопряженной по отношению к А, если а* = а. (здесь чертой обозначен переход к комплексно сопряженной величине).

Норма ||А|| квадратной матрицы А в данной работе определяется равенством

II .11

= итах а..

Пусть дана квадратная матрица А = | а^ и-го порядка. Ненулевой вектор, представляемый матрицей-столбцом ||уу|| у = 1,...,и , называется собственным вектором матрицы А, если имеют место равенства

аг1у1 =^г, 8,1 =1,к,и, (7)

где под повторяющимся индексом подразумевается суммирование. Числа X называются собственными значениями матрицы А. Данная

система однородных линеиных уравнении имеет ненулевое решение, если определитель, составленный из коэффициентов при неизвестных Уу, равняется нулю:

1 1 а12 • а1п

ёе! а 21 а 22 — Х п сч .. . а — 0

ап1 ап2 ■ • апп —Х

Это соотношение называется характеристическим уравнением матрицы А. Раскрыв определитель в последнем равенстве и расположив полученное выражение по убывающим степеням Х, получаем характеристическое уравнение в виде

(—1)п Х п + (—1)п 1 а1Х п 1 +... — а п—1Х + а п — 0.

(8)

Коэффициенты а у в этом уравнении равны суммам главных миноров у-го порядка определителя ёе! А, т.е.

а1 — а11 + а22 + к + апп, а 2

а.. а.

ау ау

= !(аиауу — а«ал^ ., ап — ёе!А.

у >

С другой стороны, как известно [52], коэффициенты а у являются элементарными симметрическими многочленами относительно корнеи Хг, I — 1,...,п, характеристического уравнения:

а1 — Х1 + Х 2 + ••• + Х п , а 2 —1Х1Х у , •••, а п —Х1Х 2 кХ п .

1 ^ У

Собственные значения X г являются решениями характеристического уравнения (8), поэтому их еще называют характеристическими числами матрицы А.

Теорема 1 (Гамильтона-Кэли). Каждая квадратная матрица удовлетворяет своему характеристическому уравнению [51, 52]:

Ап — Р1 Ап—1 — . — рп—1А — рп1 — 0, (9)

где I - единичная матрица. Здесь использованы обозначения

Р, — (—1) у 1а у, у — 1,.,п.

(10)

Теорема 2. Если функция f (X) разлагается в степенной ряд

в круге, |Х - Х0

<r

f(X)=Ха(х-Хо>', (ii)

j=0

то это разложение сохраняет силу, если скалярный аргумент X заменить любой матрицей A, характеристические числа которой лежат внутри круга сходимости. То есть если собственные значения X j матрицы A удовлетворяют условию |X j -X0| < r, j = 1,2,к,n, где r - радиус круга сходимости ряда (11), то функция f (A) определена и удовлетворяет следующему соотношению:

да

f (A) = £а, (A -X о I)j .

j=0

Следствие. Для любой квадратной матрицы A и скаляра z экспонента exp(Az) существует и может быть представлена сходящимся степенным рядом [50, 52-54]:

exp( Az) = I + Az +1( Az)2 + i( Az)3 + к (12)

Методы, основанные на аппроксимации ряда конечной суммой. Эта категория численных методов основана на прямом приложении к функциям матриц методов аппроксимации, разработанных первоначально для аналогичных функций скалярного аргумента. В этих методах ни порядок матрицы, ни ее собственные значения непосредственно не влияют на вычисления.

Усечение ряда Тейлора

да a j n a j

exPA=Z-r (13)

j=0 j=0

Этот алгоритм при некоторых значениях аргумента не является удовлетворительным даже в скалярном случае. Ошибки округления при выполнении арифметических операций по формуле (13) могут приводить к значениям элементов матрицы exp A , отличающихся от истинных значений не только по величине, но и по знаку. Оценки числа сла-

гаемых N при которых ошибка, обусловленная усечением ряда, не превышает допустимых пределов, сделаны в работах [55-57].

Рациональная аппроксимация по Чебышеву. Общая задача приближения функции /(х) в интервале (а,Ь) заключается в том, чтобы найти функцию g(х), наилучшим образом приближающую / (х) в этом интервале. Приближение по Чебышеву определено условием сделать как можно меньшим наибольшее отклонение между /(х) и g(х) (см., например, [58, с. 714]). Пусть сдд(х) - отношение

двух полиномов порядка д. Рассмотрим шах0<х<ж сдд - ехр(-х)|. Коэффициенты частного сдд для различных значений д, минимизирующие максимум, были определены в работе [59]. Эти результаты могут быть

|cqq (A) - exp(^

если

непосредственно распространены на оценку

матрица A эрмитова (т.е. выполняется равенство A' = A) и ее собственные числа Л, < 0. Аналогичные оценки для других видов матриц

неизвестны.

Аппроксимация Паде. Изложение метода, разработанного Паде [60] для приближения функций скалярного аргумента с помощью рациональных функций, имеется в обзоре [61]. Матричная экспонента expA аппроксимируется выражением [62-65]

Rpq =

Dpq ( A)Npq (A),

где

N (А)_у (р+д-/М А/ П (А)_у (р+д-П4 (-)-рд( ) /_0(р+д)!Др-/)! , рд( ) /_0(р+д)/(я-/)!( ) .

Невырожденность матрицы Брд (А) обеспечивается при достаточно

больших значениях р и д или если собственные значения матрицы А отрицательные. Как и в предыдущем методе, ошибки округления делают аппроксимацию Паде ненадежной. При больших значениях д матрица Бдд (А) аппроксимирует ряд, в который разлагается

ехр(-А /2), тогда как Ыдд (А) - ехр(-А /2). Поэтому аннулирующие

ошибки могут помешать точному определению этих матриц. Анало-

гичные замечания относятся к общим (р, д) аппроксимациям. Кроме того, матрица Орд (А) может оказаться плохо обусловленной к операции инверсии. Опасность этого особенно велика при большом разбросе собственных значений матрицы А.

Масштабирование. Трудности, связанные с ошибками округления, и затраты машинного времени при вычислении матричной экспоненты по предыдущим трем методам возрастают по мере увеличения нормы матрицы А или с возрастанием разброса собственных значений матрицы А. Обе эти трудности можно контролировать, если использовать фундаментальное свойство экспоненциальной функции:

Идея состоит в том, чтобы выбрать такое целое к, для которого матрица К _ ехр(А / к) может быть надежно и эффективно найдена,

а затем вычислить К. Один из часто используемых критериев для выбора к состоит в том, чтобы сделать его наименьшей степенью двойки, такой, что И / к < 1. При этом условии матрица К _ ехр(А / к) может быть удовлетворительно вычислена по методу Тейлора или Паде. Среди тех, кто предложил некоторые аналитические оценки данного метода или его усовершенствования, были Уорд [66] и Скрэтон [67].

Методы, основанные на численном интегрировании обыкновенных дифференциальных уравнений. Величина exp(Az)у0 является решением системы обыкновенных дифференциальных уравнений (3), удовлетворяющим начальным условиям ¥(0) _ ¥0, при условии,

что роль матрицы Ж в (3) играет матрица А с постоянными коэффициентами. Поэтому методы, основанные на численном интегрировании дифференциальных уравнений (см., например, [68-70]), являются естественным подходом к задаче вычисления матричной экспоненты. Все эти методы просты в использовании при вычислении ехр(Аг) и требуют очень небольшого дополнительного программирования. Основным недостатком данной группы методов является относительно большое время компьютерных вычислений.

Примеры вычисления матричной экспоненты с использованием программ, предназначенных для решения обыкновенных дифференци-

альных уравнений общего вида, рассмотрены в обзоре [71]. По оценке авторов указанной работы, трудоемкость такой процедуры примерно в 10 раз больше той, что требуется по методу масштабирования.

Более специализированные подходы оказываются более эффективными. Два классических подхода к решению дифференциальных уравнений - это методы Тейлора и Рунге-Кутты четвертого порядка с фиксированным шагом. Вся область изменения переменной г разбивается на т равных участков, так что 0 = г0 и г = гт. Применительно

к уравнению-----= Л¥ метод Тейлора дает

ёг

¥ =

т 1+1

( и2 И Л

I + ИЛ +—Л2 +—Л3 +—Л4

ч 2! 3! 4! ,

¥,= г4(4 .,

а метод Рунге-Кутты

¥ =¥ 1 1 1 1

¥ і+1 = ¥. +— с1 +— с2 +— с3 +— с4.

1+1 1 6 1 3 2 3 3 6 4

Здесь использованы следующие обозначения: ¥ . = ¥(2 і), И = 2І+1 - 2 і и

с1 = ИЛ¥ І , с2 = ИЛ [V І + 2 с1 ], с3 = ИЛ [V І +

Несложные преобразования показывают, что в этом случае оба метода приводят к формально одинаковому результату, но с разными ошибками округления. При фиксированном шаге И матрицу Т4(ИЛ) необходимо вычислять лишь один раз, и затем ¥ 1+1 находится в результате одного умножения матрицы на вектор столбец ¥ і . Стандартный

метод Рунге-Кутты требует 4 таких умножения на каждом шаге.

Если вычислить ¥(2) для одного частного значения г, скажем

2 = 1, то при И = 1/ т находим

¥(1) = ¥(тИ) = ¥т =[Т4(ИЛ)]т ¥0.

Это соотношение указывает на тесную связь рассматриваемых алгоритмов с методом масштабирования [72, 73]. «Масштабированной» матрицей в данном случае является ИЛ, а ее экспонента аппроксимируется матрицей Т4(ИЛ). Методы имеют одинаковые свойства

с точки зрения ошибок округления, поэтому метод Рунге-Кутты с фиксированной величиной шага не имеет преимуществ перед методом масштабирования.

Обсуждение возможности использования алгоритмов с переменным шагом интегрирования и многошаговых методов численного интегрирования к вычислению матричной экспоненты дано в обзоре [71].

Методы преобразования матриц основаны на преобразовании подобия вида

A = SBS-

При каждом преобразовании подобия сохраняется результат сложения матриц, умножения матриц и умножения матрицы на скаляр. Поэтому из разложения (12) экспоненты в ряд следует

exp Az = S (exp Bz)S_1.

Идея состоит в том, чтобы найти такую матрицу B, которая приводит к легко вычисляемой экспоненте exp(Bz ) .

Метод собственных векторов. Наивный подход заключается в выборе такой матрицы S, столбцами которой являются собственные векторы матрицы A, так что n уравнений (7) могут быть записаны в виде

AV = VD,

где матрица D - диагональная, элементы которой dii = 'Ki. Вычисление

экспоненты матрицы D выполняется элементарно, если имеется удовлетворительный метод вычисления экспонент от собственных значений. exp(Dz) - это диагональная матрица, элементы которой равны exp \. Поскольку матрица V неособенная,

exp( Az ) = V exp( Dz)V-1.

Этот метод хорошо работает для любой симметричной матрицы (aij = aji). Теоретическая трудность возникает, когда A не имеет полного набора линейно независимых собственных векторов и, таким образом является дефектной. В этом случае матрица V-1 не существует и алгоритм ломается.

Для расчетов собственных значений и собственных векторов матриц применяются специальные методы, среди которых наибольшую известность получили методы Крылова [74], Хессенберга [75, с. 273], Самуэльсона [76], Данилевского [77], эскалаторный метод [78]. Подробное изложение этого и других методов решения задачи о собственных значениях имеется в книге [75].

-алгоритм вычисления собственных векторов. Улучшение предыдущего метода с точки зрения эффективности и надежности может быть достигнуто, если собственные векторы вычисляются по ОЯ-алгоритму [79].

Предположим на время, что все собственные значения матрицы Л действительны. Идея состоит в том, чтобы найти ортогональную матрицу Q (такую, что QT Q = I) и треугольную Л, (у нее а^ = 0, если г > ] ) такие, что

QTЛQ = Л.

Это является преобразованием подобия, так как Q- = ^ . Поэтому искомые собственные значения матрицы Л находятся на главной диагонали матрицы Л . После этого нужно найти собственные векторы матрицы Л. Следующим шагом является вычисление матрицы К и диагональной Б, которые удовлетворяют условию

Л К = КБ.

После этого собственные векторы матрицы Л находятся простым перемножением матриц Q и К: V = QR . Ключевым здесь является то, что К - это верхняя треугольная матрица.

Одна из возможных реализаций указанного метода состоит в использовании специальных компьютерных программ, например ОКТНЕБ и НОЯ2 пакета Е1БРАСК [80].

Если матрица Л имеет комплексные собственные значения, то матрица К будет квазидиагональной (на главной диагонали будут стоять блоки размерностью 2 х 2, соответствующие каждой паре комплексно сопряженных собственных значений). В этом случае при вычислениях необходимо использовать комплексную арифметику.

Жорданова каноническая форма [50, с. 144], [75, с. 71]. Для любой квадратной матрицы Л существует такая невырожденная матрица С, которая преобразует Л к блочно-диагональному виду

J1 Ф J2 к.® Jl =

- квадратная матрица, порядок которой равен кратности собственного значения Xк матрицы Л. При таком преобразовании экспонента ехр( Л2) может быть вычислена по формуле

0 0 • Л ■ 0 0 = С-1 ЛС, где Jk = хк 1 0 0 хк 1 ••• 0 0 м

0 ". 0 1

0 0 ■■ Jг 0 ". 0 К 0 хк

еЛ =

С (е^ Ф.

JlZ

) С ".

в которой матрицы е к определяются согласно следующему примеру:

1 2 22/ 23/

1 2 /2! /3!

если Jk =

х к 1 0 0

0 хк 1 0

0 0 хк 1

0 0 0 хк

то е к = е

01 2 2

00 1 00 0

2!

Главная трудность такого подхода состоит в том, что при близких величинах собственных значений матрицы кратность собственных значений может быть найдена неверно, что приведет к изменению всей структуры матриц 3 и С [82, 83].

Преобразование Шура Л = QAОт состоит в нахождении треугольной (квазитреугольной) матрицы Л с помощью ортогональной (унитарной) матрицы Q. Сложности вычисления ехр( Лг) =

= Q [ ехР( Л г)] От связаны с проблемами вычисления экспоненты треугольной матрицы [84].

Тр еугольная блочно-диагональная матрица.

В этом методе используется не ортогональная, но хорошо обусловленная к операции инверсии матрица В результате преобразования получается одновременно матрица треугольная и блочно-диагональная

матрица В. С точки зрения нахождения преобразующей матрицы данный метод точнее метода «Жорданова каноническая форма», но уступает подходу «Преобразование Шура». Зато вычисляется ехр(Вг) легче, чем ехр(Лг) в предыдущем методе и, конечно, сложнее, чем в3'2 в методе канонической жордановой формы.

Полиномиальные методы. Из теоремы Гамильтона-Кэли (9) следует, что любая целая степень 1 > 0 матрицы Л может быть выражена через п степеней, например Л0 = I, Л, Л2, ..., Л”-1:

Л =”г п„Л'. (14)

1=0

Это представление и разложение (12) позволяют, в свою очередь, выразить матричную экспоненту ехр(Лг) в виде полинома

Методы вычислений ехр(Лг), которые представлены в данном разделе, основаны на использовании представления (15). Если собственные значения Х1, Х2, ..., Хп матрицы Л известны, для вычисления ехр( Лг) могут быть использованы следующие три метода.

Формула Лагранжа-Сильвестра

ехр( Лг)=2 ехр(х ^ )пЛ _У.

1 =0 1 * 1 К ]■ К1

Эта формула применима, если матрица Л не имеет кратных собственных значений. В иных случаях формула Лагранжа-Сильвестра имеет не столь простую форму [50]. Более того, если матрица имеет не равные, но близкие по величине собственные значения, расчет даже по специальным формулам не дает удовлетворительных результатов. Это замечание справедливо и по отношению к формулам двух следующих методов.

Формула Бэкера [58, с. 196] выражает функции матриц через собственные значения матрицы с помощью определителя Вандер-монда Ау = ёе1;| |Я, 1-1|. Применительно к матричной экспоненте данный метод дает

V І=0

где

А1 (ехР(Х))

1+1 1

1 Х1 ••• Х1 ехр(Х1і) Х

1 Х2 ••. X-1 ехр(Х21) Х2+1

1+1 2

X"

X 2

X 2

1 X п ■■■ Х!_' ехр(Х пг) Х2

это определитель, получаемый из Ау заменой Л1, X2, ., X” соответственно на ехр(Х1 г), ехр(Х2г), ..., ехр(Хпг). Применение данного метода к матрицам с кратными собственными значениями рассмотрено в работе [84].

Формула Ньютона

ехр( Лг) = ехр(Х1 г) + £[х„...,х 1 ]П ( л -V )■

І=2

1=1

где разностные отношения [Х^...,X^ ] зависят от г и определяются рекуррентными формулами

еХ 1 - еХ2

Х1 -Х 2

[X,,..., х„ ] = [Х|-' ХХ]-[ХХ 'Х І+1 ], і > 2.

Х1 -Х і+1

Обсуждение данного метода для случая близких собственных значений дано в книге [85].

П р е о б р а з о в а н и е Л а п л а с а. Изображение матричной экспоненты ехр(Аі), получаемое в результате преобразования Лапласа

и обозначаемое символами Л[ехр(Аі)], - это матрица (( - А) 1, элементы которой являются рациональными функциями переменной с;, а именно

Л [ехр( Аі)] = ( - А)-1 =2

|-1 гп- і-1

----------АІ..

і=о с(д)

Здесь матрицы Аі определяются рекуррентными формулами

А. = АА - р1 , Ао = 1,

п_1

с(д) = ёй (д1 _ Л) = дп _^ р^д1,

1=0

а коэффициенты рг могут быть вычислены по методу Д.К. Фаддеева (см., например, [50, с. 90] или [75, с. 295]):

Обратное преобразование Лапласа соотношения (16) дает искомую функцию

в степенной ряд по г [86]. Иные подходы к использованию преобразования Лапласа в задаче вычисления матричной экспоненты рассмотрены в статьях [87-90].

Общим для последних четырех методов является то, что все они требуют для своей реализации проведения вычислений с точностью 0(п4) . Это делает их применение, за исключением случаев матриц, для которых известны собственные значения, запредельно дорогими. Объем памяти, требуемый для хранения матриц, на порядок больше, чем требуется в методах, явно не зависящих от собственных значений матриц.

Оценка различных методов нахождения ехр(Лг) может быть проведена по таким параметрам, как общность, надежность, стабильность, точность, простота, легкость использования и эффективность численного алгоритма, объем памяти компьютера, требуемый для реализации алгоритма, и «аналитичность» метода. При этом надо иметь в виду, что на практике каждый из методов может быть реализован в виде различных компьютерных программ, отличающихся деталями и эффективностью использования.

Чем шире класс матриц, к которым применим конкретный алгоритм, тем выше общность этого метода вычислений. Именно этот параметр сыграл свою роль в отборе методов для данного обзора.

может быть разложено

При определении таких терминов, как надежность, стабильность и точность, важно проводить различие между ошибками, обусловленными самой задачей, и ошибками, свойственными конкретному алгоритму решения этой проблемы. Например, обращение почти вырожденной матрицы по своей природе является задачей очень чувствительной к малейшим изменениям параметров. Такие проблемы, как считают, относятся к разряду плохо обусловленных задач. Нет алгоритма работы с конечной арифметической точностью, который бы позволил провести вычисления на компьютере обратной матрицы, не отягощенные большими ошибками. Ошибки округления могут быть увеличены за счет больших множителей, что сделает компьютерный результат полностью ошибочным.

Алгоритм считается надежным, если он позволяет избежать чрезмерных ошибок вычислений или предупредить их возникновение.

Алгоритм устойчив, если не задает большую чувствительность к малым изменениям параметров, чем заложено самой задачей. Устойчивый алгоритм выдает ответ, который является точным для задачи с близкими данными. Метод может быть неустойчивым и при этом надежным, если нестабильность может быть обнаружена.

Точность алгоритма относится прежде всего к ошибке, которая возникает в результате усечения рядов. Эта ошибка является одним из компонентов, но не единственным, влияющим на точность полученного ответа. Часто при наличии большого компьютерного времени повышают точность вычислений, что делает метод устойчивым. Например, точностью итерационного метода решения системы уравнений можно управлять, изменяя число итераций.

Сравнение различных методов по погрешностям, вызванным округлением результатов арифметических операций, можно проводить по числу используемых элементарных операций умножений (чисел).

Эффективность измеряется количеством машинного времени, необходимым для решения конкретной задачи. Применительно к вычислению экспоненты следует различать расчеты exp A и exp(Az) для нескольких значений z. При оценке времени, необходимого на матричные вычисления, традиционно оценивается число используемых элементарных умножений (либо суммарное число элементарных сложений и умножений) и вводится некоторый коэффициент для учета остальных операций.

Наиболее конкурентоспособным по таким параметрам, как общность, надежность, стабильность и точность, среди представленных выше подходов является метод масштабирования.

Под «аналитичностью» метода понимается возможность получения аналитических оценок при проведении вычислений или нахождения аналитических решений с помощью данного метода. Принципиальная возможность анализа ехр(Лг), как функции собственных значений X. матрицы Л и параметра г заложена в формулах Лагранжа-

Сильвестра, Бэкера и Ньютона.

Еще один метод вычисления матричной экспоненты, который, с одной стороны, не уступает методу масштабирования по общности, надежности и стабильности и имеет более высокую точность, а с другой - превосходит полиномиальные методы по аналитическим возможностям, рассмотрим более подробно.

3. Метод симметрических многочленов (МСМ)

Этот метод относится к числу полиномиальных, но принципиально отличается от вышеперечисленных тем, что для вычисления целых функций матрицы не требует нахождения решений характеристического уравнения, а основан лишь на использовании коэффициентов этого уравнения.

Центральным понятием МСМ являются симметрические многочлены п-го порядка Ву матрицы Л. Они определяются рекуррентными

формулами [91, с. 63]

В =

0, если у = 0,1,..., п _ 2,

1, если 1 = п _ 1, (17)

_ Р1В1_1 + Р2В1 _2 + к + РпВ1 _п, если 1 ^ п

в которых коэффициенты р 1 выражаются через элементарные симметрические многочлены матрицы Л соотношениями (10).

Там, где необходимо уточнение, что речь идет о симметрических многочленах конкретной матрицы, скажем Л, вместо обозначений В 1 ,

О . , Ру соответственно используются Ву (Л), О . (Л), Ру (Л).

Представление функций матриц. Теорема 3 [92, с. 274]. Любая целочисленная степень у невырожденной матрицы Л п-го порядка может быть представлена с помощью ее симметрических многочленов Ву через п первых степеней Л0 = I, Л, Лп_1 формулой

Л = 1 (+п_1 _ (1В1+п_2 _ к _ Рп_1 В1 ) +

+Л (В;+п_2 _ Р1В]+п_3 _ . _ Рп_2В] ) + . + )п 3 [В]+2 _ Р1В] +1 _ Р2В] ) + (1 8)

+Лп_2 в _ Р1В1)+=£а ±Рп^В__,.

1=0 я=0

Для значений 1 = 0,1,., п _ 1, формула (18) легко проверяется, а для других у - доказывается по индукции. Если ёе! Л = 0, соотношение (18) выполняется для 1 > 0.

Соотношение (18) играет важную роль не только с точки зрения представления аналитических функций матриц, но и с точки зрения численных расчетов характеристических матриц слоисто периодических сред. Дело в том, что начиная с у = п +1 вычисление Л1 по формуле (18) становится существенно менее затратным в сравнении с обычным перемножением матриц как по числу элементарных операций умножения, так и по числу сложений.

Очевидным следствием Теорем 2 и 3 является представление целых функций матриц с помощью симметрических многочленов [93]. В частности [94],

ехр( Лг) = 2 ( Лг У-1

1=0 1 •

I го 1

1 +1 £ Р_+«(Лг)2 “Г 1« ( Лг)

я=0 1=п 1•

(19)

Если выразить симметрические многочлены п-го порядка Ву явным образом через собственные значения X матрицы Л, то суммирование по у

в формуле (19) может быть выполнено аналитически без особого труда. Результат суммирования в этом случае выражается через собственные значения матрицы. Его конвертация к функциональной зависимости от элементарных симметрических многочленов о у гарантированно выполняется для матриц второго, третьего и четвертого порядков.

Экспонента матрицы второго порядка. Симметрические многочлены второго порядка В у определяются соотношениями

Bo = 0,

Bi = 1;

(2O)

В1 О1В1 _1 О 2 В1 _ 2, где о1 и о 2 - соответственно след и определитель матрицы. Решением уравнений (20) являются функции

B л 2 - X л ©і ±. а2

B, = —------, где л12 = —±is, s = Л а2--

1 -А/ 1,2 2 V 2 4

г - мнимая единица. Подстановка этих функций в формулу (19) дает [95]

(21)

exp(Az) = exp I

Пример 1. Горизонтальная волна сдвига. Рассмотрим распространение плоской волны упругих деформаций в однородном слое твердого тела, если колебания происходят перпендикулярно плоскости падения хі волны (вдоль оси у). Определяющее уравнение и матрица переноса для данного случая выглядят следующим образом:

d

dz

uy = о I M uy ^ T = cos(kz cos 9) 1 sin(kz cos 9) Mk cos 9

Pyz -м^0 Pyz -Mk cos 9 sin(kz cos 9) cos(kz cos 9)

Здесь компонентами вектор-функции ¥ являются смещение uy и компонента тензора напряжений pyz; kz = к cos 9, к = ш^/р/р, к -

волновое число; 9 - угол падения; ю - циклическая частота; р - плотность; ц - модуль сдвига.

Столь же просто находится матрица переноса продольных колебаний в жидкости [23].

Экспонента матрицы четвертого порядка. Представим здесь выражение экспоненты матрицы четвертого порядка A в важном для практических приложений случае, когда ее элементарные симметрические многочлены ог. удовлетворяют соотношениям

ат = а3 = 0, о2 ^ 0, о4 ^ 0.

(22)

При этих условиях все многочлены четвертого порядка Бу с четными индексами равны нулю, а многочлены с нечетными индексами определяются равенствами

Б1 = 0, Б3 = 1 Б2у+1 = —°2Б2у-1 — ° 4Б2у-3’ У — 2 •

Решение этих уравнений можно представить в виде [96]

В2 ] +1 (°4)

бій [ і агссоБ Ь ]

VI - Ь2

Л

, где Ь = -

2>/°4

Использование этих функций в формуле (19) дает следующий результат:

ехр(Аг) = /50 + Л31 + А2Б2 + А3Б3,

где

(23)

£0 = сЬ(г$( ))соб(г$(+)) -

Ь

бЬ(г$( ))біп(г$(+))

£1 = 1 + 2Ь 8Ь(г$( ))со8(г$(+)) +1—'2+ЬьсЬ(г$( ))8Іп(г$(+)),

2$

23(-) 1

2_

^3 =

= КГ2

вЬ(г$( ))бій(г$(+)), " 1

(24)

1

^(+) сЬ(г$( ))8Іп(гд( ))--^-у8Ь(г$( ^сОБ^В^))

и 3( ±) =

±02 )/2.

Пр имер 2 • Волна Р - ЗУ типа. Если колебания происходят в плоскости падения, то определяющее уравнение при том же выборе системы координат, что и в предыдущем примере, имеет вид

ё_

ёг

Рхгг Рхг

Ргг = Ж Р гг

и и

X X

иг иг

0 -ікхХ 2д + X -И2р + к2 4^(Д + Х) х 2д + Х 0

-ікх 0 0 -ш2р

1 0 0 -ікх

0 1 -ікх 0

2д + X 2д + X

Р

Рг

иг

2

2

Здесь X и ^ - коэффициенты Ламе; кх = к sin 9. Нетрудно убедиться,

что определяющая матрица W в этом примере удовлетворяет условиям

(22). Поэтому матрица переноса T = exp( Wz) легко вычисляется в соответствии с формулами (23)-(24) [25, с. 78].

Масштабирование. Для экспоненты матрицы общего вида, порядок которой n > 4, получение аналитической зависимости от элементарных симметрических многочленов представляется невозможным. При сравнении различных методов численного расчета матричной экспоненты следует различать погрешность, вызванную усечением ряда (12), и погрешности, обусловленные методом вычисления каждого из слагаемых ряда и связанные с округлением результатов расчетов.

Относительная погрешность при аппроксимации ряда (12) первыми J слагаемыми, а именно

® Aj J-1 Aj

expA = У j j,

j=0 j ! j=0 j !

характеризуется величиной

У J4 Aj /(j!)

8j (A) = 1 -“^------- . (25)

J V ! ri® tin-l\

A /(j!)

Метод симметрических многочленов позволяет оценить погрешность усечения ряда аналитически. Покажем это на примере матрицы переноса.

Представим матрицу переноса T = exp(Wd) однородного слоя толщиной d по формуле (19) с применением масштабирующего коэффициента к и вспомогательной матрицы К = Wd / к :

T = Кк, (26)

где

К = Ч”)=|(")>.«,, <=7-

S=1 'У Pn-"< (W)% 7Bl-1-g (Т) ™

Величина S, определяемая равенством (28), имеет простой смысл. Легко заметить, что

8 = , «в■ /4)■

1 -еи (ЖУ / к)

Последнее равенство выполняется тем точнее, чем меньше относительная погрешность гп (Ж$ / к )■

Теорема 4 ■ Относительная погрешность е, обусловленная заменой точных формул (27)-(28) на приближенную

К = У (^)' 1

S? I к J ''

, (Wd)SN 1_ (Wd'

1+' 'У pn-'+‘ (Т )У j Bj-'-g {-,

(29)

■N+1

S N

не превосходит величины s <----------- —n------ при условии, что вели-

(N + 1)П ,=,(n + О

llWdll

чина параметра S = (2n -1)------ меньше единицы.

пк

Пример 3. Пусть коэффициент масштабирования к матрицы порядка n = 4 равен наименьшему целому, удовлетворяющему условию к > 10(2n - 1)d max |wgl| . В этом случае 0,1, и вычисление по

формуле (29) со значением N = 2 дает матрицу К с относительной погрешностью 8 менее 2 • 10-6.

Замечание. После того, как вспомогательная матрица К найдена, наиболее точная процедура вычисления матрицы переноса T = Кк, если к > п +1, состоит в использовании теоремы (18).

Пример на составление определяющего уравнения и нахождение матрицы переноса. Найдем матрицу, характеризующую распространение термоупругих волн в изотропном слое. Влияние тепловых возмущений на механические движения проявляются при

(P - SV )-поляризации вектора упругих перемещений. Выберем ту же

ориентацию декартовой системы координат, что и в предыдущих примерах. Ненулевыми компонентами вектора упругих перемещений будут их и uz. Уравнения движения в этом случае имеют вид

рд Ч = дРхх , dPxz

dt2 йх dz ’ (30)

pd 2uz = dPxz + dPzz

dt2 dx dz

гДе Рх, Рхг, Р^2 - компоненты тензора напряжении; р - плотность слоя; і - время.

Обозначим через X и р коэффициенты упругости Ламе, аТ - коэффициент линеИного расширения, А = 2ц + Х, АТ = (3Х + 2д)аТ , 9 = Т - Т0 - отклонение температуры тела Т от его температуры Т0 при

нулевых деформациях и напряжениях. Тогда, используя определение линеИного тензора деформации и закон Дюамеля-Неймана, получаем следующие три уравнения:

Р. = А ^ + Ф - Ат в-

дх дг

Р- = "&+!) (3')

Р. = Х^+А ^ - А в

дх дг

Уравнения (30) и (31) будем решать совместно с уравнениями

= , д0 Ух X Т ~ ,

дх (32)

= ,50 4 '

02

выражающими закон Фурье, и уравнением теплопроводности

|^ + 4То Г|£ + І.£Ц § = 0. (33)

дх дг удх ді дг ді ] ді

Здесь цх и - компоненты вектора плотности теплового потока; X, -коэффициент теплопроводности; сТ - теплоемкость.

Если деформации в слое порождаются падающей на него плоской волной вида ехр [/(кхх + кгі - юі)], то функции

Рхг, Р2, Рх, их, и2,9, , цх имеют одинаковую функциональную за-

висимость от координаты х и времени і :

11Р» Р2 их и2 0 ^ Рх Ях\\ = |к ^2 ^3 ^4 ^5 ^6 ^7 ^8 || еХР [/(кхх - Юі)] .

Из восьми неизвестных уг = уг (г) независимыми являются

шесть. Функции уг = уг (г), г = 1, к, 6, находятся в результате решения

уравнения вида (3), определяющая матрица которого выглядит следующим образом:

Ж =

0 ,кхА к24д(д + А) 2 -г — х 4 -рш2 А А 0 Д А |тА 0

-гкх 0 0 -рш2 0 0

1 Д 0 0 -гкх 0 0

0 1 А А Т 0 Ат А 0

0 0 0 0 0 -1 Ат

0 гш т0 Ат А -2шк,дТ„ А 0 . (т0 а? гш т ( А + с? У- У 0

(34)

А функции и выражаются через у2, у3 и равенствами

^7 = А V2 + гК

А -

2

V + АТ (А - 5, У8 =~1кх Хт V 5 .

Нетрудно получить значения элементарных симметрических многочленов матрицы Wd

^(Ш) = 0,

о3Ж) = 0,

а 2(Ш) = а 4(Wd) =

ш р . 2 ш“р | . ш

----------к,. + '

+ I

2д А Ат

( Т Ат2 ^ + ст - 2к2

1 А Т У

а5Ж) = 0,

d2

Д

У

ш2р + . ш (Т0А

+ ст

ш р . ш —+г—

( А Ат

т а?

+ Ст

' \ ЛЛ

уу

- 2к2

d4 + К d4 -

к?d4 + /рш ШСтd4,

А Ат

2

Л

- к

( Д У

к4 - к2 - г—Гт^ + с.

)1р

А

Ат

^к2 + /рш? шст

У

А Ат

d6.

Из определения (17) и равенств (35) следует, что все симметрические многочлены n -го порядка с четными индексами равны нулю, а Bj (Wd) с нечетными индексами выражаются равенствами

B1 (Wd) = B3 (Wd) = 0; B5 (Wd) = 1; B} (Wd) = -£ o2l (Wd)B}._2l (Wd), j > 7.

i=1

Дальнейшее вычисление характеристической матрицы по формуле (19) или методом масштабирования (26)-(28) не представляет труда.

Рассмотрим частный случай. Обозначим А0 - длину волны и будем считать, что выполняется условие

П = ^ < 1. (36)

X

Элементарный симметрический многочлен о j (Wd) матрицы Wd n-го порядка равен сумме главных миноров j-го порядка определителя матрицы Wd. Таких миноров Cj. Минор j-го порядка содержит j! слагаемых, каждое из которых является произведением j элементов матрицы Wd. Обозначим наибольшее по модулю значение элемента ajt матрицы A через max|а^|. Очевидно, что модуль |оj(Wd)| элементарного симметрического многочлена о j (Wd) не превосходит

Cjj!(max aJ{) . Поэтому в соответствии с определением (10) коэффициентов Pj (Wd) получаем следующие соотношения:

\Pj (Wd ^ ^ PjM (Wd) = (n П !j)! d (maX \agl I) j , j = ^ ^ ‘ ‘ ‘ , n.

Используя эту оценку и определение (17) симметрических многочленов, можно показать, что

\в, (Wd) < P1m (A) ^ ~nn_~ P1m (A) Т .

Очевидно, что при выполнении условия (36) модули симметрических многочленов |,Bj (Wd)| с индексами j = 5,7, к являются величинами (j - 5) -го порядка малости по безразмерному параметру Q, определенному соотношением (36). Поэтому

I!

I го 1

£ л-.«Ж )£ т В , (wd)

Я=0 1=п J !

<о3

и с точностью до членов третьего порядка малости по величинам О матрица переноса термоупругих волн может быть аппроксимирована

формулой Т ~ £ (Wd)1 / (1!).

1=0

Заключение

Метод численного расчета целых степеней 1 > п матрицы А п-го порядка общего вида, основанный на применении симметрических многочленов п-го порядка Ву, использует наименьшее число элементарных операций умножений и поэтому является наиболее точным по ошибкам округления в сравнении с другими известными подходами. Этот показатель при численном расчете характеристической матрицы многослойной периодической структуры является решающим при выборе метода вычислений.

МСМ применительно к вычислению экспоненты матрицы общего вида по надежности и точности вычислений, простоте программирования и использования имеет существенные преимущества перед другими подходами.

Эти преимущества обусловлены тем, что:

а) не требуется вычисление (оценка) собственных значений матрицы;

б) надежно контролируются ошибки усечения рядов;

в) аналитически определяется коэффициент масштабирования к;

г) вычисление Кк при к > п проводится с использованием симметрических многочленов Ву (К).

Для матрицы А второго-четвертого порядков МСМ позволяет получить формулы ехр(Аг) и другие аналитические функции /(Аг), определяемые через экспоненту, аналитически выражающие /(Аг) с помощью элементарных симметрических многочленов о у (А).

Для слоисто-периодических сред, описываемых характеристическими матрицами второго порядка, МХМ дает точные аналитические решения, выражающие, в частности, коэффициенты отражения и про-

пускания через элементы характеристической матрицы одного слоя [23, 25, 96-98] с помощью ортогональных полиномов. В отличие от теоремы Флоке [99] (аналогом которой в физике твердого тела является теорема Блоха [100, 101]), доказанной для бесконечных периодических сред, МХМ применим к среде произвольной толщины, состоящей из одного, двух и т.д., вплоть до бесконечно большого числа одинаковых слоев.

Большой интерес представляют эффекты, связанные с распространением упругих волн в слоистых структурах, состоящих из слоев с различными упругими, тепловыми, электрическими и магнитными свойствами. Обширная библиография, посвященная этим вопросам, имеется, например, в работах [102-107]. Применение метода матрицы переноса (характеристической матрицы) в сочетании с методом симметрических многочленов представляется весьма перспективным к расчетам распространения электроупругих, магнитоупругих и термоэлектромагнитоупругих волн в слоистых средах.

Библиографический список

1. Горелик Г.С. Колебания и волны. - М.; Л.: Гостехтеоретиздат, 1950. - 551 с.

2. Бреховских Л.М. Волны в слоистых средах. - М.: Изд-во АН СССР, 1957. - 503 с.

3. Гинзбург В.Л. Распространение электромагнитных волн в плазме. - М.: Наука, 1967. - 684 с.

4. Бреховских Л.М., Годин О.А. Акустика слоистых сред. - М.: Наука, 1989. - 416 с.

5. Фейнберг Е.Л. Распространение радиоволн вдоль земной поверхности. - М.: Наука, 1999. - 496 с.

6. Bullen K.E. An introduction to the theory of seismology. - Cambridge: Cambr. Univ. Press, 1963. - 381 p.

7. Ewing W.M., Jardetzky W.S., Press F. Elastic waves in layered media. - New York: McGraw-Hill Book Company. 1957. - 380 p.

8. Mead D.J. A general theory of harmonic wave propagation in linear periodic systems with multiple coupling // J. Sound Vibr. - 1973. - Vol. 27. -P. 235-260.

9. Сибиряков Б.П., Максимов Л.А. Татарников М. А. Анизотропия и дисперсия упругих волн в слоистых периодических структурах. -Новосибирск: Наука, 1980. - 72 с.

10. Шульга Н.А. Основы механики слоистых сред периодической структуры. - Киев: Наукова думка, 1981. - 200 с.

11. Meyers M.A., Mishra A., Benson D.J. Mechanical properties of nanocrystalline materials // Progress in Materials Science. - 2006. - Vol. 51. -P. 427-556.

12. Лаптева Т.В., Тарасенко О. С., Тарасенко С.В. Эффекты магнитоупругого взаимодействия при распространении сдвиговой волны в одномерном магнитном акустически гиротропном фононном кристалле // ФТТ. - 2007. - Т. 49. - Вып. 7. - С. 1210-1216.

13. Стретт Дж.В. (Лорд Рэлей). Теория звука. T. II. - М.: Гостех-теоретиздат, 1955. - 476 с.

14. Виноградова М.Б., Руденко О.В., Сухоруков А.П. Теория волн. - М.: Наука, 1979. - 383 с.

15. Молотков Л.А. Матричный метод в теории распространения волн в слоистых упругих и жидких средах. - Л.: Наука, 1984. - 201 с.

16. Abeles F. Recherehes sur la propagation des ondes electromag-netiques sinusoidals dans les milieux stratifiees. Application aux conches minces // Ann. Phys. - 1950. - Vol. 5. - P. 596-640, 706-782.

17. Борн М., Вольф Э. Основы оптики. - М.: Наука, 1970. - 856 с.

18. Berning P.H. Theory and calculations of optical thin films // Physics of thin films. Vol. 1. / Ed. G. Hass. - San Diego: Academic, 1963. -Pp. 69 -132.

19. Telen A. Design of multilayer interference filters // Physics of thin films. Vol. 5. / Ed. G. Hass, R.E. Thun. - San Diego: Academic, 1969. -Pp. 47 - 86.

20. Jacobson R. Inhomogeneous and coevaporated homogeneous films for optical applications // Physics of thin films. Vol. 8 / Ed. G. Hass, M.Y. Francombe, R.W. Hoffman. - New York: Academic, 1975. - P. 51-98.

21. Knittl Zd. Optics of thin films (An optical multilayer theory). -New York: J. Wiley & Sons, 1975. - 548 p

22. Бондарь Е.А., Шадрина Л.П. Метод определения оптических постоянных поглощающей пленки в составе слоистой системы // ОС. -

2004. - Т. 96, № 1. - С. 128-132.

23. Беляев Ю.Н. Рассеяние волн непрерывно слоистыми упругими средами // Вестн. СыктГУ. Сер. 1. Математика, механика, информатика. - 2011. - Вып. 13. - С. 75-88.

24. Беляев Ю.Н. Теория дифракции рентгеновских лучей в слоистых кристаллах: автореф. дис. ... канд. физ.-мат. наук. - М.: Изд-во МГУ, 1982. - 17 с.

25. Беляев Ю. Матричный подход теории волн к слоистым средам. - Saarbrucken: LAP Lambert Academic Publishing. 2012. - 156 c.

26. Tsu R., Esaki L. Tunneling in a finite superlattice // Appl. Phys. Lett. - 1973. - Vol. 22. - P. 562-564.

27. Trzeciakowski W. Boundary conditions and interface states in heterostructures // Phys. Rev. B. - 1988. - Vol. 38. - No. 6. - P. 4322-4325.

28. Владимиров М.Р., Кавокин А.В. Краевые электронные состояния в полупроводниковых сверхрешетках // ФТТ. - 1995. - Т. 37, № 7. - С. 2163-2188.

29. Басс Ф.Г., Булгаков А.А., Тетервов А.П. Высокочастотные свойства полупроводников со сверхрешетками. - М.: Наука, 1989. - 287 с.

30. Баскаков С.И. Радиотехнические цепи с распределенными параметрами. - М.: Высшая школа, 1980. - 152 с.

31. Элаши Ш. Волны в активных и пассивных периодических структурах // Труды Института инженеров по электротехнике и радиоэлектронике. - 1976. - Т. 64, № 12. - С. 22-59.

32. Thomson W.T. Transmission of elastic wave through a stratified solid material // J. Appl. Phys. - 1950. - Vol. 21. - No. 1. - P. 89-93.

33. Haskell N.A. The dispersion of surface waves on multilayered // Bul. Seism. Soc. Amer. - 1953. - Vol. 43. - No. 1. - P. 17-34.

34. Teitler S. Henvis B. Refraction in stratified anisotropic media // JOSA. - 1970. - Vol. 60. - P. 830-834.

35. Berreman D.W. Optics in stratified and anisotropic media: 4 x 4 matrix formulation // JOSA. - 1972. - Vol. 62. - P. 502-510.

36. Berreman D.W. Optics in smoothly varying anisotropic planar structures: application to liquid-crystal twist cells // JOSA. - 1973. -Vol. 63. - P. 1374-1380.

37. Lin-Chung P.J. Teitler S. 4 x 4 Matrix formalisms for optics in stratified anisotropic media // JOSA A. - 1984. - Vol. 1. - P. 703-705.

38. Abdulhalim I., Benguigui L., Weil R. Selective reflection by heli-coidal liquid crystals. Results of an exact calculation using 4 x 4 characteristic matrix method // J. Phys. (Paris). - 1985. - Vol. 46. - P. 815-825.

39. Abdulhalim I. Analytic propagation matrix method for anisotropic magneto-optic layered media // J. Opt. A: Pure Appl. Opt. - 2000. - P. 557-564.

40. Wohler H., Haas G., Fritsch M., Mlynski D.A. Faster 4 x 4 matrix method for uniaxial inhomogeneous media // JOSA. A. - 1988. - Vol. 5. -P.1554-1557.

41. Аржанников А.В., Кузнецов С. А. Методы расчета спектральных свойств многослойных структур на основе скрещенных решеток-поляризаторов // Журнал технической физики. - 2001. - Т. 71. -Вып. 12. - С. 1-5.

42. Andreeva M.A., Rosete K., Khapachev Yu.P. Matrix analog of Takagi equations for grazing-incidence diffraction // Phys. stat. sol. (a). -1985. - Vol. 88. - P. 455-462.

43. Андреева М.А., Росете К. Теория отражения от мессбауэров-ского зеркала. Учет послойных изменений параметров сверхтонких взаимодействий вблизи поверхности // Вестн. Моск. ун-та. Сер. 3. Физика. Астрономия. - 1986. - Т. 27, № 3. - С. 57-62.

44. Каули Дж. Физика дифракции. - М.: Мир. 1979. - 432 с.

45. Monaco S.F. The use of Chebyshev polynomials in the computation of inhomogeneous thin dielectric films with exponentially varying index // Thin Solid Films. - 1980. - Vol. 73. - P. L1-L4.

46. Tolstoy I. Effects of density stratification on sound waves // J. Geophys. Res. - 1965. -Vol. 70. - P. 6009-6015.

47. Robins A.J. Reflection of plane acoustic waves from a layer of varying density // JASA. - 1990. - Vol. 87. - P. 1546-1552.

48. Huang G.Y., Wang Y.S., Yu S.W. Stress concentration at apenny-shaped cracks in a nonhomogeneous medium under torsion // Acta Mech. -2005. - Vol. 180. - P. 107-115.

49. Itou S. Transient dynamic stress intensity factors around two rectangular cracks in nonhomogeneous interfacial layer between two elastic halfspaces under impact load // Acta Mech. - 2007. - Vol. 192. - P. 89-110.

50. Гантмахер Ф.Р. Теория матриц. - 4-е изд. - М.: Наука, 1988. -

552 с.

51. Курош А. Г. Курс высшей алгебры. - 9-е изд. - М.: Наука, 1971. - 432 с.

52. Ланкастер П. Теория матриц. - М.: Мир, 1969. - 272 с.

53. Маркус М., Минк Х. Обзор по теории матриц и матричных неравенств. - М.: Наука, 1972. - 232 с.

54. Беллман Р. Введение в теорию матриц. - М.: Наука, 1976. -

352 с.

55. Liou M.L. A novel method of evaluating transient response // Proc. IEEE. - 1966. - Vol. 54. - P. 20-23.

56. Everling W. On the evaluation of eAt by power series // Proc. IEEE. - 1967. - Vol. 55. - P. 413.

57. Bickart T.A. Matrix exponential: Approximation by truncated power series // Proc. IEEE. - 1968. - Vol. 56. - P. 372-373.

58. Анго А. Математика для электро- и радиоинженеров. - М.: Наука, 1967. - 780 с.

59. Cody W.J., Meinardus G., Varga R.S. Chebyshev rational approximation to exp(-x) in [0;ro] and applications to heat conduction problems // J. Approx. Theory. - 1969. - Vol. 2. - P. 50-65.

60. Pade H. Sur la representation approcheee d’une fonction par des fractions rationnelles. These de Doctorat presentee a l’Universite de la Sor-bonne, 1892 // Annales Scientifiques de l'Ecole Normale. (3ieme Serie). -1892. - Vol. 9. - P. 1-93.

61. Gragg W.B. The Pade Table and Its Relation to Certain Algorithms of Numerical Analysis // SIAM Review. - 1972. - Vol. 14. - No. 1. -P. 1-62.

62. Fair W., Luke Y. Pade approximations to the operator exponential // Numer. Math. - 1970. - Vol. 14. - P. 379-382.

63. Wragg A., Davies C. Computation of the exponential of a matrix I: Theoretical considerations // J. Inst. Math. Appl. - 1973. - Vol. 11. -P. 369-375.

64. Wragg A., Davies C. Computation of the exponential of a matrix II: Practical considerations // J. Inst. Math. Appl. - 1975. - Vol. 15. - P. 273-278.

65. Arioli M., Codenotti B., Fassino C., The Pade method for computing the matrix exponential // Lin. Alg. Appl. - 1996. - Vol. 240. - P. 111-130.

66. Ward R.C. Numerical computation of the matrix exponential with accuracy estimate // SIAM J. Numer. Anal. - 1977. - Vol. 14. - P. 600-610.

67. Scraton R.E. Comment on rational approximants to the matrix exponential // Electron. Lett. - 1971. - Vol. 7. - P. 260-261.

68. Камке Э. Справочник по обыкновенным дифференциальным уравнениям. - М.: Наука, 1965. - 704 с.

69. Современные методы решения обыкновенных дифференциальных уравнений / под ред. Дж. Холла, Дж. Уатта. - М.: Мир. 1979. - 312 с.

70. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. - М.: Бином. Лаборатория знаний, 2003. - 632 с.

71. Moler C., Van Loan C. Nineteen dubious ways to compute the exponential of matrix, Twenty-five years later // SIAM Review. - 2003. -Vol. 45. - No. 1. - P. 3-49.

72. Mastascusa E.J. A relation between Liou's method and fourth order Runge-Kutta method for evaluation of transient response // Proc. IEEE. -1969. - Vol. 57. - P. 803-804.

73. Whitney D.E. More similarities between Runge-Kutta and matrix exponential methods for evaluating transient response // Proc. IEEE. - 1969. -Vol. 57. - P. 2053-2054.

74. Крылов А.Н. О численном решении уравнения, которым в технических вопросах определяются частоты малых колебаний материальных систем // Изв. АН ОМЕН. - 1931. - № 4. - С. 491-539.

75. Фаддеев Д.К., Фаддеева В.Н. Вычислительные методы линейной алгебры. - М.: Наука, 1963. - 656 с.

76. Samuelson P.A. A method of determining explicitly the coefficients of the characteristic equation // Ann. Math. Statistics. - 1942. -Vol. 13. - P. 424-429.

77. Данилевский А. М. О численном решении векового уравнения // Математический сборник. - 1937. - Т. 2. - С. 169-171.

78. Morris J., Heed J.W. Lagrangian frequency equation. Fn “escalator” method for numerical solution // Aircraft Eng. - 1942. - Vol. 14. -P. 312-314, 316.

79. Уилкинсон Дж. Х. Алгебраическая проблема собственных значений. - М.: Наука. 1970. - 564 с.

80. Smith B.T., Boyle J.M., Dongarra J.J., Garbow B.S., Ikebe Y., Klema V.C., Moler C.B. Matrix Eigensystem Routines: EISPACK Guide, 2nd ed. Lecture Notes in Comput. Sci. 6. - New York: Springer-Verlag, 1976. - 551 p.

81. Golub G.H., Wilkinson J.H. Ill-conditioned eigensystems and the computation of the Jordan canonical form // SIAM Review. - 1976. -Vol. 18. - P. 578-619.

82. Kagstrom B., Ruhe A. An algorithm for numerical computation of the Jordan normal form of a complex matrix // ACM Transactions on Mathematical Software. - 1980. - Vol. 6. - P. 437-443.

83. Parlett B.N. A recurrence among the elements of functions of triangular matrices // Linear Algebra Appl. - 1976. - Vol. 14. - P. 117-121.

84. Vidysager M. A novel method of evaluating exp(At) in closed form // IEEE Trans. Automatic Control. - 1970. - Vol. AC-15. - P. 600-601.

85. MacDuffee C.C. The Theory of Matrices. - New York: Chelsea, 1956. - 128 p.

86. Liou M.L. Evaluation of the transition matrix // Proc. IEEE. -1967. - Vol. 55. - P. 228-229.

87. Bierman G.J. Finite series solutions for the transition matrix and covariance of a time invariant system // IEEE Trans. Automatic Control. -1971. - Vol. AC-16. - P. 173-175.

88. Chen C.F., Parker R.R. Generalization of Heaviside’s expansion technique to transition matrix evaluation // IEEE Trans. Educ. - 1966. -Vol. E-9. - P. 209-212.

89. Rao K.R., Ahmed N. Heaviside expansion of transition matrices // Proc. IEEE. - 1968. - Vol. 56. - P. 884-886.

90. Zakian V. Solution of homogeneous ordinary linear differential equations by numerical inversion of Laplace transforms // Electron. Lett. -1971. - Vol. 7. - P. 546-548.

91. Беляев Ю.Н. Алгебра тензоров. - Сыктывкар: Изд-во Сыктывкар. гос. ун-та, 2009. - 180 с.

92. Беляев Ю.Н. Векторный и тензорный анализ. - Сыктывкар: Изд-во Сыктывкар. гос. ун-та, 2010. - 298 с.

93. Belyayev Yu.N. Representation of matrix functions by means of symmetric polynomials // Book of abstracts of the International Conference on Algebra, Аugust 20-26 2012; Institute of mathematics NAS of Ukraine. -Kiev, 2012. - P. 30.

94. Беляев Ю.Н. Применение симметрических многочленов в решении задачи Коши // Алгебра и линейная оптимизация: тез. между-нар. конф., Екатеринбург, 14-19 мая 2012. - Екатеринбург: Изд-во УМЦ-УПИ, 2012. - С. 20-22.

95. Belyayev Yu.N. Calculations of transfer matrix by means of symmetric polynomials // Days on Diffraction 2012. Proceedings of the International Conference May 28 - June 1 2012. - Sankt Peterburg, 2012. - P. 36-41.

96. Беляев Ю.Н. Матричный метод расчета перерассеяния волн в периодической структуре // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки. - 2011. - № 2(23). - С. 142-148.

97. Belyaev Yu.N., Kolpakov A.V. On the theory of X-ray diffraction in a perfect crystal with distorted surface layer // Phys. stat. sol.(a). - 1983. -Vol. 76. - P. 641-646.

98. Беляев Ю.Н. Характеристическая матрица слоисто-периодической структуры // Вестн. СыктГУ. Сер. 1. Математика, механика, информатика. - 2010. - Вып. 11. - С. 86-91.

99. Floquet G. Sur les equations differentielles lineaires a coefficients periodiques // Ann. scientifiques de l'Ecole Normale superieure, ser. 2. -1883. - Vol. 12. - No. 1. - P. 47-88.

100. Bloch F. Uber die Quantenmechanik der Elektronen in Kristall-gittern // Z. Physik. - 1928. - Vol. 52. - P. 555-600.

101. Ашкрофт Н., Мермин Н. Физика твердого тела. Т. 1. - М.: Мир, 1979. - 400 с.

102. Chen Z., Yu S., Meng L., Lin Y. Effective properties of layered magneto-electro-elastic composites // Compos. Struct. - 2002. -Vol. 57. -P.177-182.

103. Багдасарян Г.Е., Даноян З.Н. Электромагнитоупругие волны. -Ереван: Изд-во Ереван. гос. ун-та, 2006. - 492 с.

104. Bhangale R.K., Ganesan N. Static analysis of simply supported functionally graded and layered magneto-electro-elastic plates // Int. J. Solids Struct. - 2006. - Vol. 43. - P. 3230-3253.

105. Guo S. H. A fully dynamic theory of piezoelectromagnetic waves // Acta Mech. - 2010. - Vol. 215. - P. 335-344.

106. Guo S. H. The thermo-electromagnetic waves in piezoelectric solids // Acta Mech. - 2011. - Vol. 219. - P. 231-240.

107. Xiaoshan Cao, Junping Shi Feng Jin Lamb wave propagation in the functionally graded piezoelectric-piezomagnetic material plate // Acta Mech. - 2012. - Vol. 223. - P. 1081-1091.

References

1. Gorelik G.S. Kolebaniya i volny [Vibrations and waves]. Moscow-Leningrad: Gosudarstvennoje isdatelstvo tekhniko-teoreticheskoy literatury, 1950, 551 p.

2. Brekhovskikh L.M. Volny v sloistyh sredah [Waves in layered media]. Moscow: Akademia nauk USSR, 1957, 503 p.

3. Ginzburg V.L. Rasprostranenie elektromagnitnyh voln v plazme [Propagation of electromagnetic waves in plasma]. Moscow: Nauka, 1967, 684 p.

4. Brekhovskikh L.M., Godin O.A. Akustika sloistyh sred [Acoustics of layered media]. Moscow: Nauka, 1989, 416 p.

5. Feinberg E.L. Rasprostranenie radio voln vdol’ zemnoi poverhnosti [Radio wave propagation along the Earth's surface]. Moscow: Nauka, 1999, 496 p.

6. Bullen K.E. An introduction to the theory of seismology. Cam-

bridge: Cambr. Univ. Press, 1963, 381 p.

7. Ewing W.M., Jardetzky W.S., Press F. Elastic waves in layered media. New York: McGraw-Hill Book Company, 1957, 380 p.

8. Mead D.J. A general theory of harmonic wave propagation in linear periodic systems with multiple coupling. J. Sound Vibr, 1973, vol. 27, pp. 235-260.

9. Sibiryakov B.P., Maksimov L.A., Tatarnikov M.A. Anizotropiya i dispersiya uprugih voln v sloistyh periodicheskih strukturah [Anisotropy and dispersion of elastic waves in layered periodic structures]. Novisibirsk: Nauka, 1980. 72 p.

10. Shulga N.A. Osnovy mehaniki sloistyh sred periodicheskoi struk-

tury [Fundamentals of mechanics of the periodic structure of layered me-

dia]. Kiev: Naukova dumka, 1981, 200 p.

11. Meyers M.A., Mishra A., Benson D.J. Mechanical properties of nanocrystalline materials. Progress in Materials Science, 2006, vol. 51, pp. 427-556.

12. Lapteva T.V., Tarasenko O.S., Tarasenko S.V. Effekty magnitou-prugogo vzaimodeistviya pri rasprostranenii sdvigovoi volny v odnomer-nom magnitnom akusticheskom fononnom kristalle [Magnetoelastic interaction in a shear wave propagating in a one-dimensional acoustically gyro-tropic magnetic phononic crystal]. Fizika tverdogo tela, 2007, vol. 49, no. 7, pp. 1210 - 1216.

13. Strutt J.W. (Lord Rayleigh). Teoriya zvuka [The theory of sound]. Moscow: Gosudarstvennoje isdatelstvo tekhniko-teoreticheskoy literatury,

1955, vol. II, 476 p.

14. Vinogradova M.B., Rudenko O.V., Sukhorukov A.P. Teoriya voln [Wave Theory]. Moscow: Nauka, 1979, 383 p.

15. Molotov L.A. Matrichnyi metod v teorii rasprostraneniya voln v sloistyh uprugih i gidkih sredah [Matrix method in the theory of wave propagation in layered elastic and liquid media]. Leningrad: Nauka, 1984, 201 p.

16. Abeles F. Recherehes sur la propagation des ondes electromag-netiques sinusoidals dans les milieux stratifiees. Application aux conches minces. Ann. Phys., 1950, vol. 5, pp. 596-640, 706-782.

17. Born M., Wolf E. Osnovy optiki [Principles of optics]. Moscow: Nauka, 1970, 856 p.

18. Berning P.H. Theory and calculations of optical thin films. Physics of thin films. Vol. 1. Ed. G. Hass. San Diego: Academic, 1963, pp. 69-132.

19. Telen A. Design of multilayer interference filters. Physics of thin films. Vol. 5. Ed. G. Hass, R.E. Thun. San Diego: Academic, 1969, pp. 47-86.

20. Jacobson R. Inhomogeneous and coevaporated homogeneous films for optical applications. Physics of thin films. Vol. 8. Ed. G. Hass, M.Y. Francombe, R.W. Hoffman. New York: Academic, 1975, pp. 51-98.

21. Knittl Zd. Optics of thin films (An optical multilayer theory). New York: J. Wiley & Sons, 1975, 548 p.

22. Bondar’ E.A., Shadrina L.P. A Metod opredeleniya opticheskih postoyannyh pogloshayushei plenki v sostave sloistoi sistemy [Method for determining the optical constants of an absorbing film entering into the composition of a layered system]. Optica i spectroskopiya, 2004, vol. 96, no. 1, pp. 128-132.

23. Belyayev Yu.N Rasseyanie voln nepreryvno sloistymi uprugimi sredami. [Wave scattering continuously stratified elastic media]. Vestnik Syktyvkarskogo universiteta. Ser. 1. Mathematika, mekhanika, informatika.

2011, no. 13, pp. 75-88.

24. Belyayev Yu.N. The theory of X-ray diffraction in layered crystals [Teoriya rentgenovskoi difraktsii v sloistyh kristallah]. Summary of the thesis Ph. D. in Physical and Mathematical Sciences. Moscowskiy gosu-darstvenniy universitet, 1982, 17 p.

25. Belyayev Yu. Matrichnyi podhod teorii voln k sloistym sredam [Matrix approach of wave theory to layered media]. Saarbrucken: LAP Lambert Academic Publishing, 2012, 156 p.

26. Tsu R., Esaki L. Tunneling in a finite superlattice. Appl. Phys. Lett, 1973, vol. 22, pp. 562-564.

27. Trzeciakowski W. Boundary conditions and interface states in heterostructures, Phys. Rev. B., 1988, vol. 38, no. 6, pp. 4322-4325.

28. Vladimirov M.R., Kavokin A.V Kraevye elektronnye sostoyaniya v poluprovodnikovyh sverhreshetkah [Boundary electronic states in semiconductor superlattices]. Fisika tverdogo tela, 1995, vol. 37, no. 7, pp. 2163-2188.

29. Bass F.G., Bulgakov A.A., Tetervov A.P. Vysokochastotnye svoistva poluprovodnikov so sverhreshetkami [High-frequency properties of the semiconductor superlattice]. Moscow: Nauka, 1989, 287 p.

30. Baskakov S.I. Radiotehnicheskie tsepi s raspredelennymi paramet-rami [Radio circuits with distributed parameters]. Moscow: Vysshaya shko-la, 1980, 152 p.

31. Elashi Ch. Volny v aktivnyh I passivnyh periodiceskih strukturah [Waves in active and passive periodic structures: a review]. Trudy instituta inzhenerov po elektrotekhnike i radioelektronike, 1976, vol. 64, no. 12, pp. 22-59.

32. Thomson W.T. Transmission of elastic wave through a stratified solid material. J. Appl. Phys. 1950, vol. 21, no. 1, pp. 89-93.

33. Haskell N.A. The dispersion of surface waves on multilayered. Bul. Seism. Soc. Amer. 1953, vol. 43, no. 1, pp. 17-34.

34. Teitler S. Henvis B. Refraction in stratified anisotropic media. J. Opt. Soc. Amer. 1970, vol. 6, pp. 830-834.

35. Berreman D.W. Optics in stratified and anisotropic media: 4 x 4 matrix formulation. J. Opt. Soc. Amer. 1972, vol. 62, pp. 502-510.

36. Berreman D.W. Optics in smoothly varying anisotropic planar structures: application to liquid-crystal twist cells. J. Opt. Soc. Amer. 1973, vol. 63, pp. 1374-1380.

37. Lin-Chung P.J. Teitler S. 4 x 4 Matrix formalisms for optics in stratified anisotropic media. J. Opt. Soc. Amer. A. 1984, vol. 1, pp. 703-705.

38. Abdulhalim I., Benguigui L., Weil R. Selective reflection by heli-coidal liquid crystals. Results of an exact calculation using 4 x 4 characteristic matrix method. J. Phys. (Paris). 1985, vol. 46, pp. 815-825.

39. Abdulhalim I. Analytic propagation matrix method for anisotropic magneto-optic layered media. J. Opt. A: Pure Appl. Opt. 2000, pp. 557-564.

40. Wohler H., Haas G., Fritsch M., Mlynski D.A. Faster 4 x 4 matrix method for uniaxial inhomogeneous media. J. Opt. Soc. Amer. A. 1988, vol. 5, pp. 1554-1557.

41. Arzhannikov A.V., Kuznetsov S. A. Metody rascheta spektralnyh svoistv mnogosloinyh struktur na osnove skreshennyh reshetok-polyariza-torov [Methods for calculating the spectral properties of multilayer structures based on arrays of crossed polarizer]. Gurnal tehnicheskoi fiziki. 2001, vol. 71, vyp. 12. pp. 1-5.

42. Andreeva M.A., Rosete K., Khapachev Yu.P. Matrix analog of Takagi equations for grazing-incidence diffraction. Phys. stat. sol. (a), 1985, vol. 88, pp. 455-462.

43. Andreeva M.A., Rosete K. Teoriya otrageniya ot messbauerovsko-go zerkala. Uchet posloinyh izmenenii parametrov sverhtonkih vzaimodeist-vii vblizi poverhnosti [Theory of reflection from the Mossbauer mirror. Accounting changes fibrewise hyperfine parameters close to the surface]. Vest-nikMoscowskogo universiteta. Ser. 3. Physika. Astronomiya. 1986, vol. 27, no. 3, pp. 57-62.

44. Cowley J. Fizika difraktsii [Diffraction physics]. Moscow: Mir, 1979, 432 p.

45. Monaco S.F. The use of Chebyshev polynomials in the computation of inhomogeneous thin dielectric films with exponentially varying index. Thin Solid Films, 1980, vol. 73, pp. L1-L4.

46. Tolstoy I. Effects of density stratification on sound waves. J. Geo-phys. Res, 1965. vol. 70, pp. 6009-6015.

47. Robins A.J. Reflection of plane acoustic waves from a layer of varying density. J. Acoust. Soc. Amer., 1990, vol. 87, pp. 1546-1552.

48. Huang G.Y., Wang Y.S., Yu S.W. Stress concentration at apenny-shaped cracks in a nonhomogeneous medium under torsion. Acta Mech.,

2005, vol. 180, pp. 107-115.

49. Itou S. Transient dynamic stress intensity factors around two rectangular cracks in nonhomogeneous interfacial layer between two elastic half-spaces under impact load. Acta Mech., 2007, vol. 192, pp. 89-110.

50. Gantmakher F.R. Teoriya matrits [Matrix theory]. Moscow: Nauka, 1988, 552 p.

51. Kurosh A.G. Kurs vysshei algebry [Course in higher algebra]. Moscow: Nauka, 1971, 432 p.

52. Lankaster P. Teoriya matrits [Theory of matrices]. Moscow: Mir, 1969, 272 p.

53. Marcus M., Minc H. Obzor po teorii matrits I matrichnyh neravenstv [Matrix theory and matrix inequalities]. Moscow: Nauka, 1972, 232 p.

54. Bellman R. Vvedenie v teoriyu matrits [Introduction to matrix analysis]. Moscow: Nauka, 1976, 352 p.

55. Liou M.L. A novel method of evaluating transient response. Proc. IEEE, 1966, vol. 54, pp. 20-23.

56. Everling W. On the evaluation of eAt by power series. Proc. IEEE, 1967, vol. 55, p. 413.

57. Bickart T.A. Matrix exponential: Approximation by truncated power series. Proc. IEEE, 1968, vol. 56, pp. 372-373.

58. Angot A. Matematika dlya elektro- i radioingenerov [Mathematics for electrical and radio engineers]. Moscow: Nauka, 1967, 780 p.

59. Cody W.J., Meinardus G., Varga R.S. Chebyshev rational approximation to exp(-x) in [0;ro] and applications to heat conduction problems. J. Approx. Theory, 1969, vol. 2, pp. 50-65.

60. Pade H. Sur la representation approcheee d’une fonction par des fractions rationnelles. These de Doctorat presentee a l’Universite de la Sor-bonne, 1892. Annales Scientifiques de l'Ecole Normale. (3ieme Serie), 1892, vol. 9, pp. 1-93.

61. Gragg W.B. The Pade Table and Its Relation to Certain Algorithms of Numerical Analysis. SIAM Review, 1972, vol. 14, no. 1, pp. 1-62.

62. Fair W., Luke Y. Pade approximations to the operator exponential. Numer. Math, 1970, vol. 14, pp. 379-382.

63. Wragg A., Davies C. Computation of the exponential of a matrix I: Theoretical considerations. J. Inst. Math. Appl, 1973, vol. 11, pp. 369-375.

64. Wragg A., Davies C. Computation of the exponential of a matrix II: Practical considerations. J. Inst. Math. Appl, 1975, vol. 15, pp. 273-278.

65. Arioli M., Codenotti B., Fassino C. The Pade method for computing the matrix exponential. Lin. Alg. Applic, 1996, vol. 240, pp. 111-130.

66. Ward R.C. Numerical computation of the matrix exponential with accuracy estimate. SIAM J. Numer. Anal, 1977, vol. 14, pp. 600-610.

67. Scraton R.E. Comment on rational approximants to the matrix exponential. Electron. Lett, 1971, vol. 7, pp. 260-261.

68. Kamke E. Spravochnik po obyknovennym differentsial’nym urav-neniyam [Handbook on ordinary differential equations]. Moscow: Nauka, 1965, 704 p.

69. Sovremennye metody resheniya obyknovennyh differentsial’nyh uravnenii [Modern methods of solving ordinary differential equations]. Ed. J. Hall, J. Watt. Moscow: Mir, 1979, 312 p.

70. Bahvalov N.S., Zhidkov N.P., Kobelkov G.M. Chislennye metody [Numerical methods]. Moscow: Binom, 2003, 632 p.

71. Moler C., Van Loan C. Nineteen dubious ways to compute the exponential of matrix, Twenty-five years later. SIAM Review, 2003, vol. 45, no. 1, pp. 3-49.

72. Mastascusa E. J. A relation between Liou's method and fourth order Runge-Kutta method for evaluation of transient response. Proc. IEEE, 1969, vol. 57, pp. 803-804.

73. Whitney D. E. More similarities between Runge-Kutta and matrix exponential methods for evaluating transient response. Proc. IEEE, 1969, vol. 57, pp. 2053-2054.

74. Krylov A.N. O chislennom reshenii uravneniya, kotorym v tekhni-cheskikh voprosakh opredelyayutsya chastoty malykh kolebanii material’nykh sistem [One numerical solution of the equations that defines the frequencies of small oscillations of material systems in technical issues]. Isves-tiya Akademii nauk SSSR VIIseriya. Matematika, 1931, no. 4, pp. 491-539.

75. Faddeev D.K., Faddeev V.N. Vychislitel’nye metody lineinoi al-gebry [Computational methods of linear algebra]. Moscow: Nauka, 1963, 656 p.

76. Samuelson P.A. A method of determining explicitly the coefficients of the characteristic equation. Ann. Math. Statistics, 1942, vol. 13, pp. 424-429.

77. Danilevsky A.M. O chislennom reshenii vekovogo uravneniya [On the numerical solution of the secular equation]. Matematicheskii sbornik, 1937, vol. 2, pp. 169-171.

78. Morris J., Heed J.W. Lagrangian frequency equation. Fn “escalator” method for numerical solution. Aircraft Eng., 1942, vol. 14, pp. 312314, 316.

79. Wilkinson J.H. Algebraicheskaya problema sobstvennyh znache-niy [The algebraic eigenvalue problem]. Moscow: Nauka, 1970, 564 p.

80. Smith B.T., Boyle J.M., Dongarra J.J., Garbow B.S., Ikebe Y., Klema V.C., Moler C.B. Matrix Eigensystem Routines: EISPACK Guide, 2nd ed. Lecture Notes in Comput. Sci. 6, New York; Springer-Verlag, 1976, 551 p.

81. Golub G.H., Wilkinson J.H. Ill-conditioned eigensystems and the computation of the Jordan canonical form. SIAM Review, 1976, vol. 18, pp. 578-619.

82. Kagstrom B., Ruhe A. An algorithm for numerical computation of the Jordan normal form of a complex matrix. ACM Transactions on Mathematical Software, 1980, vol. 6, pp. 437-443.

83. Parlett B.N. A recurrence among the elements of functions of triangular matrices. Linear Algebra Appl, 1976, vol. 14, pp. 117-121.

84. Vidysager M. A novel method of evaluating exp(At) in closed form. IEEE Trans. Automatic Control, 1970, vol. AC-15, pp. 600-601.

85. MacDuffee C.C. The Theory of Matrices. New York: Chelsea,

1956, 128 p.

86. Liou M.L. Evaluation of the transition matrix. Proc. IEEE, 1967, vol. 55, pp. 228-229.

87. Bierman G.J. Finite series solutions for the transition matrix and covariance of a time invariant system. IEEE Trans. Automatic Control, 1971, vol. AC-16, pp. 173-175.

88. Chen C.F., Parker R.R. Generalization of Heaviside’s expansion technique to transition matrix evaluation. IEEE Trans. Educ, 1966, vol. E-9, pp. 209-212.

89. Rao K.R., Ahmed N. Heaviside expansion of transition matrices. Proc. IEEE, 1968, vol. 56, pp. 884-886.

90. Zakian V. Solution of homogeneous ordinary linear differential equations by numerical inversion of Laplace transforms. Electron. Lett., 1971, vol. 7, pp. 546-548.

91. Belyayev Yu.N. Algebra tenzorov [Tensor algebra]. Syktyvkarskiy gosudarstvenniy universitet, 2009, 180 p.

92. Belyayev Yu.N. Vektornyi i tenzornyi analiz [Vector and tensor analysis]. Syktyvkarskiy gosudarstvenniy universitet, 2010, 298 p.

93. Belyayev Yu.N. Representation of matrix functions by means of symmetric polynomials. Book of abstracts of the International Conference on Algebra, august 20-26 2012. Kiev: Institute of mathematics NAS of Ukraine, 2012. pp. 30.

94. Belyayev Yu.N. Primenenie simmetricheskikh mnogochlenov v reshenii zadachi Koshi [Application of symmetric polynomials in the solution of the Cauchy problem]. Tezisy Megdunarodnoi konferentsii “Algebra i

lineinaya optimizatsiya” Ekaterinburg: Uralskiy politekhnicheskiy institut,

2012. pp. 20-22.

95. Belyayev Yu.N. Calculations of transfer matrix by means of symmetric polynomials. Days on Diffraction 2012. Proceedings of the International Conference May 28, June 1 2012. Sankt Peterburg, Russia, pp. 36-41.

96. Belyayev Yu.N. Matrichnyi metod rascheta pererasseyaniya voln v periodicheskoi strukture [Matrix method of calculating the rescatteriing wave in a periodic structure]. Vestnik Samarskogo gosudarstvennogo tech-nicheskogo universiteta. Ser. Fiziko-matematicheskie nauki, 2011, no 2(23), pp. 142-148.

97. Belyaev Yu.N., Kolpakov A.V. On the theory of X-ray diffraction in a perfect crystal with distorted surface layer. Phys. stat. sol.(a), 1983, vol. 76, pp. 641-646.

98. Belyayev Yu.N. Harakteristicheskaya matritsa sloisto-periodiches-koi struktury [Characteristic matrix of layered periodic structure]. Vestnik Syktyvkar skogo gosudarstvennogo universiteta. Ser. 1. Matematika, mechanika, informatika, 2010, no. 11, pp. 86-91.

99. Floquet G. Sur les equations differentielles lineaires a coefficients periodiques. Ann. scientifiques de l'Ecole Normale superieure, ser. 2, 1883, vol. 12, no. 1, pp. 47-88.

100. Bloch F. Uber die Quantenmechanik der Elektronen in Kristall-gittern. Z. Physik, 1928, vol. 52, pp. 555-600.

101. Ashcroft N.W., Mermin D. Fizika tverdogo tela [Solid state physics]. Vol. 1. Moscow: Mir, 1979, 400 p.

102. Chen Z., Yu S., Meng L., Lin Y. Effective properties of layered magneto-electro-elastic composites. Compos. Struct, 2002, vol. 57, pp. 177-182.

103. Bagdasaryan G.E., Danoyan Z.N. Elektromagnitouprugie volny [Electromagnetoelasticity wave]. Erevanskiy gosudarstvenniy universitet,

2006, 492 p.

104. Bhangale R.K., Ganesan N. Static analysis of simply supported functionally graded and layered magneto-electro-elastic plates. Int. J. Solids Struct, 2006, vol. 43, pp. 3230-3253.

105. Guo S.H. A fully dynamic theory of piezoelectromagnetic waves. ActaMech., 2010, vol. 215, pp. 335-344.

106. Guo S.H. The thermo-electromagnetic waves in piezoelectric solids. Acta Mech., 2011, vol. 219, pp. 231-240.

107. Xiaoshan Cao, Junping Shi, Feng Jin Lamb wave propagation in the functionally graded piezoelectric-piezomagnetic material plate. Acta Mech, 2012, vol. 223, pp. 1081-1091.

Об авторах

Беляев Юрий Николаевич (Сыктывкар, Россия) - кандидат физико-математических наук, доцент, доцент кафедры математического моделирования и кибернетики Сыктывкарского государственного университета (167001, г. Сыктывкар, Октябрьский пр., 55, e-mail:

ybelyayev@mail.ru).

About the authors

Belyayev Yuriy Nikolaevich (Syktyvkar, Russian Federation) -Ph. D. in Physical and Mathematical Sciences, Ass. Professor, Department of Mathematical Modelling and Cybernetics, Syktyvkarskiy gosudarsrven-niy universitet (55, Oktyabrskiy av., 167001, Syktyvkar, Russian Federation, e-mail: ybelyayev@mail.ru).

Получено 26.02.2013