Full Text

Введение Среда, свойства которой постоянны на каждой плоскости, перпендикулярной к фиксированному направлению, называется плоско-слоистой средой. Если выбрать это специальное направление за ось z декартовой системы координат, то плотность среды, коэффициенты Ламе, диэлектрическая проницаемость и другие параметры среды, характеризующие упругие и электромагнитные свойства вещества, будут функциями одной этой координаты z. Естественными средами, имеющими хорошо выраженную горизонтальную стратификацию, являются океан и атмосфера. Этим объясняется, в частности, сверхдальнее распространение звука под водой и радиоволн в атмосфере [1–5]. Во многих случаях земная кора также успешно может быть аппроксимирована слоистой средой, и это используется в теории распространения сейсмических волн [6, 7]. Широким спектром физических свойств обладают слоистые среды искусственного происхождения. Класс таких объектов постоянно расширяется, и области их применения уже в настоящее время исключительно велики. Значительную роль в разнообразных технических и физических приложениях играют слоисто-периодические структуры. Исследования распространения волн различной физической природы в средах, свойства которых изменяются послойно периодическим образом, весьма обширны по тематике и охватывают большой круг важных вопросов. В механике деформируемого твердого тела интерес к этой проблеме [8, 9] связан в первую очередь с широким внедрением в современную технику слоистых композиционных материалов, в том числе нанокристаллических [10–12]. Постоянно расширяющийся класс слоистых сред и многообразие их возможных применений поддерживают неослабевающий интерес исследователей к проблемам взаимодействия волн с такими структурами. Развитию теории распространения волн в слоистых средах способствуют два равных по силе обстоятельства. Во-первых, дифракция излучений различных диапазонов длин волн лежит в основе методов неразрушающего исследования строения и состава слоев, контроля качества изготавливаемых объектов. Во-вторых, исследования закономерностей распространения волн в слоистых объектах обосновывают создание на их основе новых приборов, устройств и технических конструкций и расширение областей применения существующих. Стандартный подход математической физики к расчету распространения волн в слоистых средах (см., например, [2, 7, 13–15]) состоит в выводе соответствующего волнового уравнения, в решении этого уравнения в каждом из слоев рассматриваемой среды с последующим «сшиванием» полученных решений на границах слоев. Модель одномерного изменения параметров, хорошо аппроксимирующая многие слоистые среды, а также линейность используемых волновых уравнений обусловили широкое использование матричных методов в теории дифракции волн в слоистых средах. Одна из важных особенностей матричного подхода состоит в том, что он позволяет существенно упростить процедуру «сшивания» решений и, как следствие, значительно облегчить процесс моделирования спектров отражения от многослойных систем. 1. Метод характеристической матрицы и его аналоги Суть метода характеристической матрицы (МХМ) состоит в следующем. Компоненты волнового поля, значения которых обуславливаются на границах слоя, составляют вектор-функцию . Элементами вектор-функции в зависимости от решаемой задачи могут быть компоненты тензора напряжений, векторы смещения, составляющие электромагнитного и теплового полей. Значения этой вектор-функции на противоположных границах слоя связаны друг с другом матрицей M, элементы которой определяют величины коэффициентов отражения и пропускания волны слоистой средой, то есть характеризуют рассеивающие свойства последней. Поэтому данную матрицу M называют характеристической матрицей (ХМ) слоя. Пусть пластина занимает область пространства, ограниченного плоскостями и . Тогда (1) Если пластина состоит из N слоев, первый из которых лежит в области второй – N-й слой – а слои характеризуются матрицами ,…, , то ХМ всей пластины равна произведению ХМ отдельных слоев, образующих пластину, (2) Свойство (2) предопределяет наиболее значительные преимущества МХМ. Так как парциальные характеристические матрицы , , зависят каждая только от структуры своего слоя, то изменение толщины и/или физических параметров -го слоя оказывает влияние только на j-ю ХМ , а произведения и остаются неизменными. Поэтому, если при расчете коэффициентов отражения и пропускания волны многослойной средой необходимо учесть изменение структуры j-го слоя, частные произведения и вследствие ассоциативности произведения (2) могут быть вычислены отдельно, и нет необходимости полного перерасчета произведения (2). ХМ второго порядка была впервые предложена в 1950 г. Абелé [16] для описания рассеяния света диэлектрическим слоем (без гиротропии). Данный подход, изложение которого имеется также в книге [17] и обзорах [18–20], стал основным методом при расчете оптических свойств многослойных диэлектрических пленок. ХМ в этом случае еще называют интерференционной матрицей [21, с. 440], а сам метод – методом последовательного наращивания слоев [22]. Второй порядок ХМ достаточен также для описания распространения горизонтальной волны сдвига (так называемая SH волна [7]) и упругих волн в средах, в которых модуль сдвига = 0 [23]. Развитие МХМ применительно к двухволновой рентгеновской дифракции в слоистых кристаллах было выполнено как в геометрии Брэгга [24], так и Лауэ [25]. В исследованиях электронных состояний полупроводниковых сверхрешеток (см., например, [26–28]) очень плодотворным оказался метод матрицы переноса T. Эта матрица является обратной по отношению к характеристической матрице M: Т = М –1. Еще одно название, используемое для МХМ как в электронике, так и в радиотехнике, – это метод матрицы передачи [29, 30]. Наконец, отметим, что метод, заключающийся в получении характеристической матрицы второго порядка, связывающей выходные и входные параметры единичного элемента дискретной периодической структуры, успешно используется при разнообразных технических расчетах. Ссылки на соответствующие работы имеются в обзоре [31]. Матрицы четвертого и более высоких порядков. Характеристическая матрица четвертого порядка была введена для расчета распространения волн упругих напряжений в слоисто-неоднородных средах с плоскопараллельными границами в работах Томсона [32] и Хаскела [33]. Описание данного метода и развитие его в подход, использующий матрицы пятого и шестого порядка, дано в монографии [18]. Развитие МХМ до матриц четвертого порядка в оптике слоистых анизотропных сред было сделано Тейтлером и Хенвисом [34] и применено Берреманом к исследованию жидких кристаллов [35, 36]. Матрицы четвертого порядка были использованы, в частности, для расчета перерассеяния электромагнитных волн слоистыми анизотропными средами в оптическом [37–40] и СВЧ [41] диапазонах, для рассмотрения рентгеновской и мессбауэровской дифракции в условиях полного отражения [42, 43]. При дифракции излучения рентгеновского диапазона частот в кристаллах возможно возникновение трех и большего числа сильных волн, а при рассеянии электронов на тонких кристаллических пленках число дифракционных пучков, появляющихся одновременно, может достигать нескольких сотен. Соответственно, n-волновой дифракции соответствует матрица переноса n-го порядка, которая в теории электронографии носит название матрицы рассеяния [44]. Определяющее уравнение матрицы переноса. Теоретические исследования распространения волн, например упругих деформаций, основывается на решении соответствующих волновых уравнений. Обычно из исходной системы дифференциальных уравнений в частных производных первого порядка относительно компонент поля (компоненты тензора напряжений, вектора смещений, напряженностей электрического и магнитного полей, температуры и др.) выводится уравнение второго или более высокого порядка относительно одной из компонент поля или соответствующих потенциалов [2, 10, 15, 17]. Из решений полученных уравнений для каждого слоя структуры составляются элементы матрицы , удовлетворяющей условию (1). Аналитические решения в таком подходе были получены для характеристической матрицы второго порядка для случаев однородных по составу слоев [2, 16, 17], для слоев с линейным [21] и экспоненциальным [45] законами изменения параметров среды по толщине слоя. Кроме перечисленных существует еще небольшое число случаев (см., например, [2, 3, 15, 46, 47]), в которых соответствующие волновые уравнения имеют аналитические решения и, соответственно, может быть получен явный вид ХМ. Для нахождения ХМ второго порядка слоистой среды с непрерывно изменяющимися параметрами по толщине в работе [23] использован метод, основанный на решении соответствующих интегральных уравнений Вольтерра методом последовательных приближений. Аппроксимация неоднородного слоя набором слоев с линейным изменением упругих свойств была предложена в статье [48]. Основным теоретическим подходом к расчету ХМ непрерывно слоистой среды является аппроксимация последней наборами однородных слоев [2, 15, 17, 49]. Альтернативный подход к расчету характеристической матрицы слоя состоит в решении системы дифференциальных уравнений первого порядка (3) которая может быть получена из исходных уравнений, например с помощью Фурье-преобразования компонент поля по времени и одной (или двух) координат. Порядок n системы (3) равен числу независимых компонент вектор-функции и, соответственно, порядку матрицы W. Решение уравнения (3) выражается через начальные значения с помощью матрицы T: (4) В теории дифференциальных уравнений матрица T, выражающая значение неизвестной вектор-функции через начальные данные , называется фундаментальной. В теории дифракции матрица T, «транслирующая» вектор на расстояние d, известна как матрица переноса. Значения элементов этой матрицы определяются свойствами слоя рассматриваемой среды, то есть матрицей коэффициентов W системы дифференциальных уравнений (3). Поэтому матрицу W, как и систему (3), будем называть определяющей. Из сравнения соотношений (1) и (4) очевидно, что матрица, обратная T, является характеристической: (5) Наиболее простой связью, определяющей матрицы W с матрицей переноса T и, следовательно, с характеристической матрицей M, является их связь для однородного слоя. В этом случае [51, с. 115] (6) Один из способов решения системы уравнений (3) для неоднородного по толщине слоя состоит в аппроксимации последнего достаточно большим числом однородных слоев и использовании для матриц переноса каждого из этих слоев соотношения (6), так что где – толщины слоев. Существуют десятки методов численного расчета матричной экспоненты. Их обзор представлен в следующем разделе. 2. Методы вычислений матричной экспоненты Эффективность некоторых алгоритмов вычислений экспоненты матрицы сильно зависит от вида матрицы. Существенную роль в некоторых подходах играют собственные значения и собственные векторы матрицы, даже если они не вовлечены в конкретный алгоритм. Напомним их определения, а также некоторые другие понятия, обозначения и теоремы, связанные с рассматриваемым вопросом. Некоторые сведения из теории матриц. Для обозначения элементов матрицы (таблицы чисел), скажем A, используются соответствующие строчные буквы с двумя индексами. Сама матрица изображается в виде таблицы, заключенной в прямые двойные скобки и обозначается символами . Первый индекс элемента матрицы указывает номер строки, а второй – номер столбца, на пересечении которых в таблице находится данное число. Элементы квадратной матрицы называются диагональными. Сумма диагональных элементов – это след матрицы. Порядок квадратной матрицы n – это число ее строк (столбцов). С алгеброй матриц можно познакомиться по любому учебнику, например [52]. Определитель матрицы A обозначается . Матрица A – невырожденная, если Квадратная матрица I, у которой диагональные элементы равны единице, а остальные нули, называется единичной. Матрица A–1, удовлетворяющая условию называется обратной по отношению к матрице А. Если поменять в матрице строки и столбцы местами, то в результате получается транспонированная матрица . Матрица А* называется сопряженной по отношению к A, если (здесь чертой обозначен переход к комплексно сопряженной величине). Норма квадратной матрицы A в данной работе определяется равенством . Пусть дана квадратная матрица n-го порядка. Ненулевой вектор, представляемый матрицей-столбцом , называется собственным вектором матрицы A, если имеют место равенства (7) где под повторяющимся индексом подразумевается суммирование. Числа называются собственными значениями матрицы A. Данная система однородных линейных уравнений имеет ненулевое решение, если определитель, составленный из коэффициентов при неизвестных , равняется нулю: . Это соотношение называется характеристическим уравнением матрицы A. Раскрыв определитель в последнем равенстве и расположив полученное выражение по убывающим степеням , получаем характеристическое уравнение в виде (8) Коэффициенты в этом уравнении равны суммам главных миноров j-го порядка определителя , т.е. С другой стороны, как известно [52], коэффициенты являются элементарными симметрическими многочленами относительно корней , характеристического уравнения: Собственные значения являются решениями характеристического уравнения (8), поэтому их еще называют характеристическими числами матрицы A. Теорема 1 (Гамильтона–Кэли). Каждая квадратная матрица удовлетворяет своему характеристическому уравнению [51, 52]: (9) где I – единичная матрица. Здесь использованы обозначения (10) Теорема 2. Если функция разлагается в степенной ряд в круге, (11) то это разложение сохраняет силу, если скалярный аргумент l заменить любой матрицей A, характеристические числа которой лежат внутри круга сходимости. То есть если собственные значения матрицы A удовлетворяют условию , где r – радиус круга сходимости ряда (11), то функция определена и удовлетворяет следующему соотношению: . Следствие. Для любой квадратной матрицы A и скаляра z экспонента существует и может быть представлена сходящимся степенным рядом [50, 52–54]: (12) Методы, основанные на аппроксимации ряда конечной суммой. Эта категория численных методов основана на прямом приложении к функциям матриц методов аппроксимации, разработанных первоначально для аналогичных функций скалярного аргумента. В этих методах ни порядок матрицы, ни ее собственные значения непосредственно не влияют на вычисления. Усечение ряда Тейлора (13) Этот алгоритм при некоторых значениях аргумента не является удовлетворительным даже в скалярном случае. Ошибки округления при выполнении арифметических операций по формуле (13) могут приводить к значениям элементов матрицы , отличающихся от истинных значений не только по величине, но и по знаку. Оценки числа слагаемых N, при которых ошибка, обусловленная усечением ряда, не превышает допустимых пределов, сделаны в работах [55–57]. Рациональная аппроксимация по Чебышеву. Общая задача приближения функции в интервале (a,b) заключается в том, чтобы найти функцию , наилучшим образом приближающую в этом интервале. Приближение по Чебышеву определено условием сделать как можно меньшим наибольшее отклонение между и (см., например, [58, с. 714]). Пусть – отношение двух полиномов порядка q. Рассмотрим . Коэффициенты частного для различных значений q, минимизирующие максимум, были определены в работе [59]. Эти результаты могут быть непосредственно распространены на оценку , если матрица A эрмитова (т.е. выполняется равенство ) и ее собственные числа . Аналогичные оценки для других видов матриц неизвестны. Аппроксимация Паде. Изложение метода, разработанного Падé [60] для приближения функций скалярного аргумента с помощью рациональных функций, имеется в обзоре [61]. Матричная экспонента аппроксимируется выражением [62–65] где Невырожденность матрицы обеспечивается при достаточно больших значениях p и q или если собственные значения матрицы A отрицательные. Как и в предыдущем методе, ошибки округления делают аппроксимацию Падé ненадежной. При больших значениях q матрица аппроксимирует ряд, в который разлагается , тогда как – Поэтому аннулирующие ошибки могут помешать точному определению этих матриц. Аналогичные замечания относятся к общим (p, q) аппроксимациям. Кроме того, матрица может оказаться плохо обусловленной к операции инверсии. Опасность этого особенно велика при большом разбросе собственных значений матрицы A. Масштабирование. Трудности, связанные с ошибками округления, и затраты машинного времени при вычислении матричной экспоненты по предыдущим трем методам возрастают по мере увеличения нормы матрицы A или с возрастанием разброса собственных значений матрицы A. Обе эти трудности можно контролировать, если использовать фундаментальное свойство экспоненциальной функции: Идея состоит в том, чтобы выбрать такое целое k, для которого матрица может быть надежно и эффективно найдена, а затем вычислить . Один из часто используемых критериев для выбора k состоит в том, чтобы сделать его наименьшей степенью двойки, такой, что . При этом условии матрица может быть удовлетворительно вычислена по методу Тейлора или Паде. Среди тех, кто предложил некоторые аналитические оценки данного метода или его усовершенствования, были Уорд [66] и Скрэтон [67]. Методы, основанные на численном интегрировании обыкновенных дифференциальных уравнений. Величина является решением системы обыкновенных дифференциальных уравнений (3), удовлетворяющим начальным условиям , при условии, что роль матрицы W в (3) играет матрица A с постоянными коэффициентами. Поэтому методы, основанные на численном интегрировании дифференциальных уравнений (см., например, [68–70]), являются естественным подходом к задаче вычисления матричной экспоненты. Все эти методы просты в использовании при вычислении и требуют очень небольшого дополнительного программирования. Основным недостатком данной группы методов является относительно большое время компьютерных вычислений. Примеры вычисления матричной экспоненты с использованием программ, предназначенных для решения обыкновенных дифференциальных уравнений общего вида, рассмотрены в обзоре [71]. По оценке авторов указанной работы, трудоемкость такой процедуры примерно в 10 раз больше той, что требуется по методу масштабирования. Более специализированные подходы оказываются более эффективными. Два классических подхода к решению дифференциальных уравнений – это методы Тейлора и Рунге–Кутты четвертого порядка с фиксированным шагом. Вся область изменения переменной z разбивается на m равных участков, так что и . Применительно к уравнению метод Тейлора дает а метод Рунге–Кутты Здесь использованы следующие обозначения: , и Несложные преобразования показывают, что в этом случае оба метода приводят к формально одинаковому результату, но с разными ошибками округления. При фиксированном шаге h матрицу необходимо вычислять лишь один раз, и затем находится в результате одного умножения матрицы на вектор столбец . Стандартный метод Рунге–Кутты требует 4 таких умножения на каждом шаге. Если вычислить для одного частного значения z, скажем , то при находим Это соотношение указывает на тесную связь рассматриваемых алгоритмов с методом масштабирования [72, 73]. «Масштабированной» матрицей в данном случае является hA, а ее экспонента аппроксимируется матрицей . Методы имеют одинаковые свойства с точки зрения ошибок округления, поэтому метод Рунге–Кутты с фиксированной величиной шага не имеет преимуществ перед методом масштабирования. Обсуждение возможности использования алгоритмов с переменным шагом интегрирования и многошаговых методов численного интегрирования к вычислению матричной экспоненты дано в обзоре [71]. Методы преобразования матриц основаны на преобразовании подобия вида При каждом преобразовании подобия сохраняется результат сложения матриц, умножения матриц и умножения матрицы на скаляр. Поэтому из разложения (12) экспоненты в ряд следует Идея состоит в том, чтобы найти такую матрицу B, которая приводит к легко вычисляемой экспоненте . Метод собственных векторов. Наивный подход заключается в выборе такой матрицы S, столбцами которой являются собственные векторы матрицы A, так что n уравнений (7) могут быть записаны в виде где матрица D – диагональная, элементы которой . Вычисление экспоненты матрицы D выполняется элементарно, если имеется удовлетворительный метод вычисления экспонент от собственных значений. – это диагональная матрица, элементы которой равны Поскольку матрица V неособенная, Этот метод хорошо работает для любой симметричной матрицы (). Теоретическая трудность возникает, когда A не имеет полного набора линейно независимых собственных векторов и, таким образом является дефектной. В этом случае матрица не существует и алгоритм ломается. Для расчетов собственных значений и собственных векторов матриц применяются специальные методы, среди которых наибольшую известность получили методы Крылова [74], Хессенберга [75, с. 273], Самуэльсона [76], Данилевского [77], эскалаторный метод [78]. Подробное изложение этого и других методов решения задачи о собственных значениях имеется в книге [75]. QR-алгоритм вычисления собственных векторов. Улучшение предыдущего метода с точки зрения эффективности и надежности может быть достигнуто, если собственные векторы вычисляются по QR-алгоритму [79]. Предположим на время, что все собственные значения матрицы A действительны. Идея состоит в том, чтобы найти ортогональную матрицу Q (такую, что ) и треугольную , (у нее если ) такие, что Это является преобразованием подобия, так как . Поэтому искомые собственные значения матрицы A находятся на главной диагонали матрицы . После этого нужно найти собственные векторы матрицы . Следующим шагом является вычисление матрицы R и диагональной D, которые удовлетворяют условию После этого собственные векторы матрицы A находятся простым перемножением матриц Q и R: . Ключевым здесь является то, что R – это верхняя треугольная матрица. Одна из возможных реализаций указанного метода состоит в использовании специальных компьютерных программ, например ORTHES и HQR2 пакета EISPACK [80]. Если матрица A имеет комплексные собственные значения, то матрица R будет квазидиагональной (на главной диагонали будут стоять блоки размерностью соответствующие каждой паре комплексно сопряженных собственных значений). В этом случае при вычислениях необходимо использовать комплексную арифметику. Жорданова каноническая форма [50, с. 144], [75, с. 71]. Для любой квадратной матрицы A существует такая невырожденная матрица C, которая преобразует A к блочно-диагональному виду , где – квадратная матрица, порядок которой равен кратности собственного значения матрицы A. При таком преобразовании экспонента может быть вычислена по формуле в которой матрицы определяются согласно следующему примеру: если , то . Главная трудность такого подхода состоит в том, что при близких величинах собственных значений матрицы кратность собственных значений может быть найдена неверно, что приведет к изменению всей структуры матриц J и C [82, 83]. Преобразование Шура состоит в нахождении треугольной (квазитреугольной) матрицы с помощью ортогональной (унитарной) матрицы Q. Сложности вычисления связаны с проблемами вычисления экспоненты треугольной матрицы [84]. Треугольная блочно-диагональная матрица. В этом методе используется не ортогональная, но хорошо обусловленная к операции инверсии матрица S. В результате преобразования получается одновременно матрица треугольная и блочно-диагональная матрица B. С точки зрения нахождения преобразующей матрицы данный метод точнее метода «Жорданова каноническая форма», но уступает подходу «Преобразование Шура». Зато вычисляется легче, чем в предыдущем методе и, конечно, сложнее, чем в методе канонической жордановой формы. Полиномиальные методы. Из теоремы Гамильтона–Кэли (9) следует, что любая целая степень матрицы A может быть выражена через n степеней, например , A,, ..., : . (14) Это представление и разложение (12) позволяют, в свою очередь, выразить матричную экспоненту в виде полинома . (15) Методы вычислений , которые представлены в данном разделе, основаны на использовании представления (15). Если собственные значения l1, l2, …, ln матрицы A известны, для вычисления могут быть использованы следующие три метода. Формула Лагранжа–Сильвестра Эта формула применима, если матрица A не имеет кратных собственных значений. В иных случаях формула Лагранжа–Сильвестра имеет не столь простую форму [50]. Более того, если матрица имеет не равные, но близкие по величине собственные значения, расчет даже по специальным формулам не дает удовлетворительных результатов. Это замечание справедливо и по отношению к формулам двух следующих методов. Формула Бэкера [58, с. 196] выражает функции матриц через собственные значения матрицы с помощью определителя Вандермонда Применительно к матричной экспоненте данный метод дает , где это определитель, получаемый из заменой , , …, соответственно на , , …, . Применение данного метода к матрицам с кратными собственными значениями рассмотрено в работе [84]. Формула Ньютона , где разностные отношения зависят от z и определяются рекуррентными формулами Обсуждение данного метода для случая близких собственных значений дано в книге [85]. П р е о б р а з о в а н и е Л а п л а с а. Изображение матричной экспоненты , получаемое в результате преобразования Лапласа и обозначаемое символами , – это матрица , элементы которой являются рациональными функциями переменной а именно . (16) Здесь матрицы определяются рекуррентными формулами , а коэффициенты могут быть вычислены по методу Д.К. Фаддеева (см., например, [50, с. 90] или [75, с. 295]): Обратное преобразование Лапласа соотношения (16) дает искомую функцию . Обратное преобразование может быть разложено в степенной ряд по z [86]. Иные подходы к использованию преобразования Лапласа в задаче вычисления матричной экспоненты рассмотрены в статьях [87–90]. Общим для последних четырех методов является то, что все они требуют для своей реализации проведения вычислений с точностью . Это делает их применение, за исключением случаев матриц, для которых известны собственные значения, запредельно дорогими. Объем памяти, требуемый для хранения матриц, на порядок больше, чем требуется в методах, явно не зависящих от собственных значений матриц. Оценка различных методов нахождения может быть проведена по таким параметрам, как общность, надежность, стабильность, точность, простота, легкость использования и эффективность численного алгоритма, объем памяти компьютера, требуемый для реализации алгоритма, и «аналитичность» метода. При этом надо иметь в виду, что на практике каждый из методов может быть реализован в виде различных компьютерных программ, отличающихся деталями и эффективностью использования. Чем шире класс матриц, к которым применим конкретный алгоритм, тем выше общность этого метода вычислений. Именно этот параметр сыграл свою роль в отборе методов для данного обзора. При определении таких терминов, как надежность, стабильность и точность, важно проводить различие между ошибками, обусловленными самой задачей, и ошибками, свойственными конкретному алгоритму решения этой проблемы. Например, обращение почти вырожденной матрицы по своей природе является задачей очень чувствительной к малейшим изменениям параметров. Такие проблемы, как считают, относятся к разряду плохо обусловленных задач. Нет алгоритма работы с конечной арифметической точностью, который бы позволил провести вычисления на компьютере обратной матрицы, не отягощенные большими ошибками. Ошибки округления могут быть увеличены за счет больших множителей, что сделает компьютерный результат полностью ошибочным. Алгоритм считается надежным, если он позволяет избежать чрезмерных ошибок вычислений или предупредить их возникновение. Алгоритм устойчив, если не задает большую чувствительность к малым изменениям параметров, чем заложено самой задачей. Устойчивый алгоритм выдает ответ, который является точным для задачи с близкими данными. Метод может быть неустойчивым и при этом надежным, если нестабильность может быть обнаружена. Точность алгоритма относится прежде всего к ошибке, которая возникает в результате усечения рядов. Эта ошибка является одним из компонентов, но не единственным, влияющим на точность полученного ответа. Часто при наличии большого компьютерного времени повышают точность вычислений, что делает метод устойчивым. Например, точностью итерационного метода решения системы уравнений можно управлять, изменяя число итераций. Сравнение различных методов по погрешностям, вызванным округлением результатов арифметических операций, можно проводить по числу используемых элементарных операций умножений (чисел). Эффективность измеряется количеством машинного времени, необходимым для решения конкретной задачи. Применительно к вычислению экспоненты следует различать расчеты и для нескольких значений z. При оценке времени, необходимого на матричные вычисления, традиционно оценивается число используемых элементарных умножений (либо суммарное число элементарных сложений и умножений) и вводится некоторый коэффициент для учета остальных операций. Наиболее конкурентоспособным по таким параметрам, как общность, надежность, стабильность и точность, среди представленных выше подходов является метод масштабирования. Под «аналитичностью» метода понимается возможность получения аналитических оценок при проведении вычислений или нахождения аналитических решений с помощью данного метода. Принципиальная возможность анализа , как функции собственных значений матрицы A и параметра z заложена в формулах Лагранжа–Сильвестра, Бэкера и Ньютона. Еще один метод вычисления матричной экспоненты, который, с одной стороны, не уступает методу масштабирования по общности, надежности и стабильности и имеет более высокую точность, а с другой – превосходит полиномиальные методы по аналитическим возможностям, рассмотрим более подробно. 3. Метод симметрических многочленов (МСМ) Этот метод относится к числу полиномиальных, но принципиально отличается от вышеперечисленных тем, что для вычисления целых функций матрицы не требует нахождения решений характеристического уравнения, а основан лишь на использовании коэффициентов этого уравнения. Центральным понятием МСМ являются симметрические многочлены n-го порядка матрицы A. Они определяются рекуррентными формулами [91, с. 63] (17) в которых коэффициенты выражаются через элементарные симметрические многочлены матрицы A соотношениями (10). Там, где необходимо уточнение, что речь идет о симметрических многочленах конкретной матрицы, скажем A, вместо обозначений , , соответственно используются , , . Представление функций матриц. Теорема 3 [92, с. 274]. Любая целочисленная степень j невырожденной матрицы A n-го порядка может быть представлена с помощью ее симметрических многочленов через n первых степеней A, …, формулой (18) Для значений , формула (18) легко проверяется, а для других j – доказывается по индукции. Если соотношение (18) выполняется для Соотношение (18) играет важную роль не только с точки зрения представления аналитических функций матриц, но и с точки зрения численных расчетов характеристических матриц слоисто периодических сред. Дело в том, что начиная с вычисление по формуле (18) становится существенно менее затратным в сравнении с обычным перемножением матриц как по числу элементарных операций умножения, так и по числу сложений. Очевидным следствием Теорем 2 и 3 является представление целых функций матриц с помощью симметрических многочленов [93]. В частности [94], . (19) Если выразить симметрические многочлены n-го порядка явным образом через собственные значения матрицы A, то суммирование по j в формуле (19) может быть выполнено аналитически без особого труда. Результат суммирования в этом случае выражается через собственные значения матрицы. Его конвертация к функциональной зависимости от элементарных симметрических многочленов гарантированно выполняется для матриц второго, третьего и четвертого порядков. Экспонента матрицы второго порядка. Симметрические многочлены второго порядка определяются соотношениями (20) где и – соответственно след и определитель матрицы. Решением уравнений (20) являются функции где , i – мнимая единица. Подстановка этих функций в формулу (19) дает [95] (21) Пример 1. Горизонтальная волна сдвига. Рассмотрим распространение плоской волны упругих деформаций в однородном слое твердого тела, если колебания происходят перпендикулярно плоскости падения волны (вдоль оси y). Определяющее уравнение и матрица переноса для данного случая выглядят следующим образом: . Здесь компонентами вектор-функции являются смещение и компонента тензора напряжений ; , k – волновое число; θ – угол падения; ω – циклическая частота; ρ – плотность; μ – модуль сдвига. Столь же просто находится матрица переноса продольных колебаний в жидкости [23]. Экспонента матрицы четвертого порядка. Представим здесь выражение экспоненты матрицы четвертого порядка A в важном для практических приложений случае, когда ее элементарные симметрические многочлены удовлетворяют соотношениям , , . (22) При этих условиях все многочлены четвертого порядка с четными индексами равны нулю, а многочлены с нечетными индексами определяются равенствами . Решение этих уравнений можно представить в виде [96] , где . Использование этих функций в формуле (19) дает следующий результат: , (23) где (24) и . Пример 2. Волна P – SV типа. Если колебания происходят в плоскости падения, то определяющее уравнение при том же выборе системы координат, что и в предыдущем примере, имеет вид Здесь и – коэффициенты Ламе; Нетрудно убедиться, что определяющая матрица W в этом примере удовлетворяет условиям (22). Поэтому матрица переноса легко вычисляется в соответствии с формулами (23)–(24) [25, с. 78]. Масштабирование. Для экспоненты матрицы общего вида, порядок которой , получение аналитической зависимости от элементарных симметрических многочленов представляется невозможным. При сравнении различных методов численного расчета матричной экспоненты следует различать погрешность, вызванную усечением ряда (12), и погрешности, обусловленные методом вычисления каждого из слагаемых ряда и связанные с округлением результатов расчетов. Относительная погрешность при аппроксимации ряда (12) первыми слагаемыми, а именно характеризуется величиной (25) Метод симметрических многочленов позволяет оценить погрешность усечения ряда аналитически. Покажем это на примере матрицы переноса. Представим матрицу переноса однородного слоя толщиной d по формуле (19) с применением масштабирующего коэффициента k и вспомогательной матрицы : , (26) где , (27) (28) Величина , определяемая равенством (28), имеет простой смысл. Легко заметить, что Последнее равенство выполняется тем точнее, чем меньше относительная погрешность Теорема 4. Относительная погрешность ε, обусловленная заменой точных формул (27)–(28) на приближенную (29) не превосходит величины при условии, что величина параметра меньше единицы. Пример 3. Пусть коэффициент масштабирования k матрицы порядка равен наименьшему целому, удовлетворяющему условию . В этом случае и вычисление по формуле (29) со значением дает матрицу K с относительной погрешностью ε менее . Замечание. После того, как вспомогательная матрица K найдена, наиболее точная процедура вычисления матрицы переноса , если , состоит в использовании теоремы (18). Пример на составление определяющего уравнения и нахождение матрицы переноса. Найдем матрицу, характеризующую распространение термоупругих волн в изотропном слое. Влияние тепловых возмущений на механические движения проявляются при ()-поляризации вектора упругих перемещений. Выберем ту же ориентацию декартовой системы координат, что и в предыдущих примерах. Ненулевыми компонентами вектора упругих перемещений будут и . Уравнения движения в этом случае имеют вид (30) где , , – компоненты тензора напряжений; ρ – плотность слоя; t – время. Обозначим через и коэффициенты упругости Ламе, – коэффициент линейного расширения, , , – отклонение температуры тела от его температуры при нулевых деформациях и напряжениях. Тогда, используя определение линейного тензора деформаций и закон Дюамеля–Неймана, получаем следующие три уравнения: (31) Уравнения (30) и (31) будем решать совместно с уравнениями (32) выражающими закон Фурье, и уравнением теплопроводности (33) Здесь и – компоненты вектора плотности теплового потока; – коэффициент теплопроводности; – теплоемкость. Если деформации в слое порождаются падающей на него плоской волной вида , то функции имеют одинаковую функциональную зависимость от координаты и времени : Из восьми неизвестных независимыми являются шесть. Функции , находятся в результате решения уравнения вида (3), определяющая матрица которого выглядит следующим образом: . (34) А функции и выражаются через , и равенствами Нетрудно получить значения элементарных симметрических многочленов матрицы (35) Из определения (17) и равенств (35) следует, что все симметрические многочлены -го порядка с четными индексами равны нулю, а с нечетными индексами выражаются равенствами Дальнейшее вычисление характеристической матрицы по формуле (19) или методом масштабирования (26)–(28) не представляет труда. Рассмотрим частный случай. Обозначим – длину волны и будем считать, что выполняется условие (36) Элементарный симметрический многочлен матрицы n-го порядка равен сумме главных миноров j-го порядка определителя матрицы . Таких миноров . Минор j-го порядка содержит слагаемых, каждое из которых является произведением j элементов матрицы . Обозначим наибольшее по модулю значение элемента матрицы через . Очевидно, что модуль элементарного симметрического многочлена не превосходит . Поэтому в соответствии с определением (10) коэффициентов получаем следующие соотношения: Используя эту оценку и определение (17) симметрических многочленов, можно показать, что Очевидно, что при выполнении условия (36) модули симметрических многочленов с индексами являются величинами -го порядка малости по безразмерному параметру , определенному соотношением (36). Поэтому и с точностью до членов третьего порядка малости по величинам матрица переноса термоупругих волн может быть аппроксимирована формулой . Заключение Метод численного расчета целых степеней матрицы A n-го порядка общего вида, основанный на применении симметрических многочленов n-го порядка использует наименьшее число элементарных операций умножений и поэтому является наиболее точным по ошибкам округления в сравнении с другими известными подходами. Этот показатель при численном расчете характеристической матрицы многослойной периодической структуры является решающим при выборе метода вычислений. МСМ применительно к вычислению экспоненты матрицы общего вида по надежности и точности вычислений, простоте программирования и использования имеет существенные преимущества перед другими подходами. Эти преимущества обусловлены тем, что: a) не требуется вычисление (оценка) собственных значений матрицы; б) надежно контролируются ошибки усечения рядов; в) аналитически определяется коэффициент масштабирования k; г) вычисление при проводится с использованием симметрических многочленов . Для матрицы A второго–четвертого порядков МСМ позволяет получить формулы и другие аналитические функции , определяемые через экспоненту, аналитически выражающие с помощью элементарных симметрических многочленов . Для слоисто-периодических сред, описываемых характеристическими матрицами второго порядка, МХМ дает точные аналитические решения, выражающие, в частности, коэффициенты отражения и пропускания через элементы характеристической матрицы одного слоя [23, 25, 96–98] с помощью ортогональных полиномов. В отличие от теоремы Флоке [99] (аналогом которой в физике твердого тела является теорема Блоха [100, 101]), доказанной для бесконечных периодических сред, МХМ применим к среде произвольной толщины, состоящей из одного, двух и т.д., вплоть до бесконечно большого числа одинаковых слоев. Большой интерес представляют эффекты, связанные с распространением упругих волн в слоистых структурах, состоящих из слоев с различными упругими, тепловыми, электрическими и магнитными свойствами. Обширная библиография, посвященная этим вопросам, имеется, например, в работах [102–107]. Применение метода матрицы переноса (характеристической матрицы) в сочетании с методом симметрических многочленов представляется весьма перспективным к расчетам распространения электроупругих, магнитоупругих и термоэлектромагнитоупругих волн в слоистых средах.

About the authors

Yuriy Nikolaevich Belyayev

Syktyvkarskiy gosudarsrvenniy universitet

Email: ybelyayev@mail.ru
55, Oktyabrskiy av., 167001, Syktyvkar, Russian Federation Ph. D. in Physical and Mathematical Sciences, Ass. Professor, Department of Mathematical Modelling and Cybernetics, Syktyvkarskiy gosudarsrvenniy universitet

References

Statistics

Views

Abstract - 445

PDF (Russian) - 343