Содержание к диссертации
Введение
Глава 1. Обзор существующих методов и постановка задачи 16
Системы линейных алгебраических уравнений 16
Численные методы решения СЛАУ 16
Метод минимальных невязок 19
Многошаговый метод минимальных невязок 21
Использование предобуславливателей 23
Экономия вычислительных ресурсов в итерациях 25
Решение интегральных уравнений 28
Рассматриваемые интегральные уравнения 28
Уравнения с вырожденным ядром и резольвента 31
Дискретизация методом Галеркина 34
Метод коллокации для кусочно-постоянной аппроксимации 35
Постановка задачи 39
Глава 2. Изложение предлагаемых методов 41
Разреженный предобуславливатель для одномерных интегральных уравнений 42
Кусочно-постоянный предобуславливатель 42
Предобуславливатель для уравнения в пространстве непрерывных функций 48
Предобуславливатель в пространстве функций, непрерывных с производной 52
Предобуславливатель для уравнений с периодическими функциями 55
Обобщение разреженного предобуславливателя на многомерные задачи 63
Структура предобуславливателя в многомерном случае 63
Предобуславливатель в пространстве многомерных непрерывных функций 66
Многомерный предобуславливатель для периодических функций 69
Предобуславливатель для систем интегральных уравнений 72
Плотный предобуславливатель для интегральных уравнений 73
Предобуславливатель на основе резольвенты 73
Плотный предобуславливатель с непрерывной резольвентой 76
Плотный предобуславливатель для систем интегральных уравнений 78
Формирование плотного предобуславливателя с внешним параметром 79
Сокращение количества операций при умножении на вектор S2
Некоторые вопросы сходимости рассматриваемых итерационных методов при решении интегральных уравнений 83
Обоснование возможности подбора предобуславливателя 83
Обход отсутствия сходимости метода минимальных невязок 88
Глава 3. Программная реализация и численные эксперименты 92
Платформа для проведения экспериментов 92
Реализованное программное обеспечение 92
Примеры интегральных уравнений с плохой сходимостью 93
Численные эксперименты 97
Формирование прототипа матрицы СЛАУ 97
Сходимость для одномерного интегрального уравнения 99
Сходимость для одномерного уравнения с периодическими функциями 109
Сходимость для двумерного интегрального уравнения 117
Соответствие сходимости для сеток разных размерностей 120
Эффективность модификации ММН и МММН 123
Выводы 124
Глава 4. Применение к задачам математического моделирования 126
Решение задачи рассеяния электромагнитных волн 126
Теоретическое описание математической модели 126
Аналитическая подготовка к вычислениям 129
Численное решение и результаты 131
Расчет звукового излучения в пространстве 138
Описание задачи математического моделирования 138
Результаты численного решения 140
Заключение 144
Список литературы 148
- Метод минимальных невязок
- Предобуславливатель для уравнения в пространстве непрерывных функций
- Реализованное программное обеспечение
- Аналитическая подготовка к вычислениям
Введение к работе
Настоящая диссертационная работа посвящена разработке новых способов и алгоритмов повышения сходимости итерационных методов при численном решении многомерных линейных интегральных уравнений, возникающих в корректно поставленных задачах математического моделирования.
Метод минимальных невязок
Метод минимальных невязок (ММН) [72, 109] относится к нестационарным итерационным методам. Он широко используется при численном решении интегральных уравнений с оператором без положительной определенности.
На каждом шаге итерационного процесса производится вычитание из вектора текущего приближения йп вектора текущей невязки hn с некоторым коэффициентом тв, причем коэффициент выбирается так, чтобы длина вектора невязки следующего шага была минимальной из всех возможных при таком подходе вариантов: Ч,+і= 7„-г,Д,; К=Аии-/; (1.6) т.. " (А.А)"
Угловыми скобками здесь обозначено скалярное произведение. Как показано в [72], минимум достигается в случае, когда следующая невязка /ги+] будет ортогональна вектору Ahu. Соответственно, коэффициент тп в приведенной формуле рассчитывается из условия их ортогональности: hll+i=Aul„l-f = A(uH,fr,)-f = h,,lAh,; (17) {h Ah -rJh Ah Ah -r Ah Q.
Итерационный процесс начинается с некоторого заданного начального приближения й0 и продолжается до тех пор, пока относительная точность вектора текущего приближения йа не достигнет некоторого заданного значения є (1.4). Заданная точность может быть никогда не достигнута, несмотря на то, что йа стабилизируется. Это может произойти в случае, если йи приблизится к значению, при котором (hn, Ahn \ и 0. Таким образом, ММН сходится для любого й0 в случае СЛАУ с такой матрицей А, при которой для любого ненулевого вектора v выполняется условие: Mv} 0, И 0. (1.8) Несмотря на то, что это не столь жесткое условие, нежели положительная определенность ((v, Av) 0, v Ф 0), возможность применения ММН к решению конкретной задачи все же требует доказательства [72].
Многошаговый метод минимальных невязок При решении СЛАУ (1.2) с помощью ММН на каждом шаге невязка последующего шага hiHl образуется, как видно из формул (1.6), путем коррекции невязки \ вектором Ahn: К А,- Л Е-тпА)К (1-9)
При этом корректирующий коэффициент тп на каждом шаге выбирается так, чтобы последующая невязка была ортогональна корректирующему вектору Ahn. Однако при такой коррекции после нескольких последовательных итераций невязка в общем случае не будет минимально возможной. Выразим невязку после М последовательных шагов ММН (M N): Кш =(Я-w K» -] =(-Wi4"(-v4R =K-f/«A%. (ЇЛО) /=і
Видно, что при образовании невязки hii+M невязка hn корректируется циклическим базисом из М векторов, порожденным оператором А. При такой коррекции, как показано в [72], невязка h)l+M будет минимальной тогда, когда корректирующие коэффициенты выбираются исходя из условия ортогональности hH+M циклическому базису AhH,...,AMhir
Многошаговый метод минимальных невязок (МММН) является усовершенствованием ММН [72]. Он основан на группировке нескольких последовательных шагов ММН и единовременном вычислении корректирующих коэффициентов с целью уменьшения невязки по сравнению с ММН.
Предобуславливатель для уравнения в пространстве непрерывных функций
Приведенный предобуславливатель на основе скалярных блоков, будучи примененным к уравнению (1.25), описанному в пространстве непрерывных функций, переводит уравнение в пространство кусочно-непрерывных функций, т.е. его решение в общем случае будет функцией, кусочно-непрерывной на границах интервалов (2.11). Однако известно, что решение уравнения (1.25) при непрерывном по совокупности аргументов ядре К\х, у) и непрерывной функции неоднородности fix) будет также непрерывной функцией. Таким образом, можно сделать предположение, что предобуславливатель, порождающий из непрерывной сеточной функции также непрерывную сеточную функцию, будет обеспечивать лучшую сходимость для уравнения в пространстве непрерывных функций. Сеточной функцией будем называть результат конечномерной аппроксимации некоторой функции. Под непрерывностью сеточной функции подразумевается непрерывность функции, которую данная сеточная функция аппроксимирует.
Опишем способ формирования на основе матрицы прототипа такого предобуславливателя, который сохранял бы непрерывность отображенной им сеточной функции. Представим обобщенную структуру приведенного ранее предобуславливателя (рис. 1).
Матрица предобуславливателя состоит из квадратных диагональных блоков. В описанном выше случае диагональные блоки являются скалярными матрицами, т.е. их диагональные элементы не изменяются в пределах блока. На границах блоков в диагоналях происходит резкое изменение значений элементов, что и порождает разрывы при отображении предобуславливателем непрерывной сеточной функции. Чтобы приблизить результат такого отображения также к непрерывной сеточной функции, можно добиться, чтобы по всем ненулевым диагоналям элементы не претерпевали резких изменений, т.е. их значения также изменялись в соответствии с некоторыми непрерывными функциями. При этом, в соответствии с использованным методом дискретизации, в центрах блоков значения должны быть приближены к соответствующим значениям малой матрицы прототипа. Таким образом, значения на диагоналях предобуславливателя могут быть получены путем дискретизации некоторых интерполяционных кривых, построенных для узловых значений, взятых с соответствующих диагоналей прототипа предобуславливателя.
Однако даже при подобном способе формирования матрицы предобуславливателя результат его умножения на непрерывную сеточную функцию в общем случае не будет непрерывной сеточной функцией, в ней все равно будут разрывы на границах интервалов (2.11). Это происходит из-за отличий произведений значений диагоналей с правого и левого края на, соответственно, элементы в начале и конце вектора. На рис. 2 обведены точки, о которых идет речь. Далее будем называть их критическими точками. Разрывы происходят вследствие отличий между значениями произведений вблизи каждых двух точек на одной горизонтали.
Чтобы устранить порождаемые этим обстоятельством разрывы, следует значения произведений элементов диагоналей и вектора на краях слева и справа также приблизить друг к другу. К примеру, можно приблизить все интерполяционные кривые, описывающие значения соответствующих диагоналей, вблизи обозначенных точек к нулю.
Интерполяционные кривые можно строить множеством способов. Поскольку сейчас рассматривается уравнение в пространстве непрерывных функций, используем кусочно-линейную интерполяцию.
При кусочно-линейной интерполяции значения на диагоналях меняются линейно между двумя центрами соседних по диагонали блоков (рис. 3). Как уже сказано выше, в центре блока значение диагонали должно быть приближено к соответствующему значению матрицы прототипа
предобуславливателя. Для крайних правых или левых блоков предобуславливателя в том направлении, где нет соседнего блока, значения диагонали линейно опускаются до нуля так, чтобы обращаться в ноль в критической точке. Это не выполняется для крайних блоков, расположенных на главной диагонали, иначе матрица предобуславливателя будет близка к вырожденной. Для крайних блоков сверху и снизу в направлении, в котором для них нет соседних блоков по диагонали, соседним значением принимается значение центра блока, ближайшего к отсутствующему. Т.е. для верхних блоков берется значение центра соседнего блока слева, для нижних -соседнего справа. Для крайних блоков на главной диагонали соседним значением принимается значение в центре этого же блока, т.е. значения главной диагонали на интервале от центра крайнего блока до границы матрицы не меняются.
Реализованное программное обеспечение
Для проведения численных экспериментов были выбраны несколько интегральных уравнений с плохой сходимостью. Для большей наглядности рассмотрение предобуславливателей ведется на примере одномерных и двумерных интегральных уравнений. Одномерные предобуславливатели рассма фиваются на примере следующего интегрального уравнения Фредгольма второго рода: ы( )- JK(x,y)u(y)cty = f(x\ а К(х,у)=еМ; (3-1) [а;Ь]=[0;х] Известно, что оператор этого интегрального уравнения не имеет нулевых собственных значений, значит, уравнение можно решать при любой функции неоднородности. На рис. 8 представлен график расположения значений функции и{х) решения (ЗЛ) на комплексной плоскости.
На примере (ЗЛ) исследуется сходимость итерационных методов при использовании разреженных и плотных предобуславливателей для уравнений с непериодическими функциями.
Для исследования сходимости итерационных методов при использовании предобуславливателей на основе ДПФ требуется интегральное уравнение с периодическими по всем аргументам функциями: ь и(х) JK(x,y)u(y)dy = f{x\ / ч simc + sinj ч Ч У)=\ Г7Г У (3,2) }posx + cosy\-[\ + i) [а;Ь] [0;2тг] Оператор уравнения (3.2) также не имеет нулевых собственных значений. График расположения значений искомой функции и(х) при указанной функции неоднородности представлен на рис. 9.
Расположение значений функции решения интегрального уравнения (3.2)
В качестве примера многомерных уравнений будем рассматривать следующее двумерное интегральное уравнение Фредгольма: Kfabw eb-»" -3 - ДадЬА (3.3) [a,: 6, ]х [а2; Ь2 ] - [О; лг]х [О; л]
Двумерный случай вполне отражает обобщение одномерных предобуславливателей на многомерные, при этом имеются достаточно широкие возможности выбора размерностей каждого измерения для полноценного описания поведения входящих функций, учитывая доступные вычислительные ресурсы.
Эксперименты будут проводиться путем численного решения уравнений (3.1), (3.2) и (3.3) при различных условиях. При этом из исходного уравнения путем дискретизации будет формироваться СЛАУ размерности N: Ай = й-Км=/. (3.4)
Эту СЛАУ будем решать итерационными методами до достижения некоторой точности (1.4) как без предобуславливателя, так и с применением предобуславливателей. При этом будут использоваться разреженные кусочно-постоянный предобуславливатель Рчс, кусочно-линейный предобуславливатель Рк!, кусочно-кубический предобуславливатель Psb и предобуславливатель на основе ДПФ .Р,, а также плотные кусочно постоянный предобуславливатель Pik, кусочно-линейный предобуславливатель Pd!, кусочно-кубический предобуславливатель Pdh и предобуславливатель на основе ДПФ PdJ. Эти обозначения введены для удобства дальнейшего описания результатов экспериментов, при этом верхний индекс будет прямо или косвенно характеризовать размерность исходного прототипа, сформированного для построения предобуславливателя.
В численных экспериментах будет производиться решение СЛАУ с использованием МММН с количеством шагов в группе, равным 2, 3 и 5. Большое количество на практике нечасто используется, поскольку с ростом количества шагов существенно увеличиваются вычислительные затраты на проведение ортогонализации методом Грама-Шмидта. Использование же меньшего количества шагов в группе (ММН) нецелесообразно, поскольку МММН с малыми значениями входного параметра не требователен к ресурсам и при этом, как правило, значительно повышает сходимость.
Результаты, полученные при проведении численных экспериментов, будут представляться в количестве произведенных операций умножения матрицы решаемой СЛАУ на вектор, поскольку в итерационных методах решения СЛАУ это, как правило, самая ресурсоемкая операция. При решении СЛАУ с помощью МММН в процессе одной итерации выполняется столько операций умножения матрицы на вектор, каков параметр МММН -количество шагов в группе.
Аналитическая подготовка к вычислениям
Рассмотрим следующий класс задач электродинамики. Пусть в конечной области Q среда характеризуется тензорами диэлектрической и магнитной проницаемости є и /и соответственно. Вне области Q параметры среды постоянны и изотропны, т.е. = Q = const и // - jua const. Требуется определить параметры электромагнитного поля, возбуждаемого в этой среде внешним полем. Описанный класс задач носит название задач рассеяния и имеет большое значение в электродинамике.
Здесь мы рассмотрим частный случай описанного класса задач, когда среда характеризуется всюду постоянным значением магнитной проницаемости /ла и тензор-функцией диэлектрической проницаемости є, отличной от постоянного тензора є J только в области Q. При этом область Q среда также однородна и изотропна, т.е. s-єі. Внешнее поле характеризуется временной зависимостью в виде множителя е ,а". Ставится задача найти векторные поля иЯ, возбуждаемые в среде Q внешними полями Е и Н . При таких исходных данных векторное поле Е определяется системой объемных интегральных уравнений, а векторное поле Й вычисляется исходя из найденного поля Ё [72].
На основе уравнений Максвелла, а также исходя из условия непрерывности тангенциальных компонент поля на границе области Q и условия излучения на бесконечности, в [72] выводится следующее интегральное уравнение относительно векторного поля Е:
Уравнение (4.2) должно быть решено относительно J. Для этого производится дискретизация методом коллокации, в процессе которой вся область Q делится на элементарные объемы П . За узловую точку в каждом элементарном объеме принимается центр П t. После дискретизации полученная СЛАУ выглядит следующим образом [72]: І/ Л--ІХ5 =- = 1АЗ. (4.4)
Здесь суммирование по q = (ql,q7 q ) ведется по всем узловым точкам Q; Jpl, Еpi - значение компоненты / соответствующей функции в узловой точке р = (р1гр2,ръ); fipllii - компоненты тензор-функции р{х)=\є{х) е0і) в узловой точке р. Мы будем рассматривать случай рассеяния волн в однородной изотропной области Q, поэтому в нашем случае тензор Р\х) будет скаляром: (AJ)pl , -±± в ,,= Е и J 1 з й т=\ д (4.5)
Значения оператора R в точках коллокации описываются следующим ооразом: в - ДОод З 2ik0 "її Х1 Iх рт : Х/ 1 Х1 АХрш Хт ІШ -8, г г \\ dx dx,. (4.6) Поскольку уравнение (4.1) описано на основе сингулярных интегралов, для случаев, когда р - q, в [72] приводится другая формула: Зр,, ,,, = А, + + а L-G(r) fa-x/fr -yJV Уф(г)р(у- 1 - ) 1 г2 ) \ж г2 dx.dx dx,. (4.7)
В формулах (4.6) и (4.7) г - расстояние от точки коллокации, т.е. от центра П , до текущего параметра интегрирования х = (х1,х2 х3): r = yj(x] xpj+(x1-xll2f+[x3-xpif. В (4.7) коэффициенты а, определяются формой границы П . Для куба a, =-,1 = 1,2,3, что совпадает со значениями этих коэффициентов для шара. Также в (4.7) используется следующая всюду дифференцируемая функция: o(r)=i+(iwv)g !7v kQr Полученный путем решения СЛАУ (4.5) вектор будет приближенно описывать значения компонент электрического тока поляризации J (4.3) в точках коллокации. Исходя из этих значений могут быть найдены значения векторного поля Е в точках коллокации, на основе которых, в свою очередь, могут быть приближенно вычислены значения компонент векторного ПОЛЯ Н по формулам, приведенным в [72], Поскольку нас интересует только решение интегрального уравнения, мы не будем останавливаться на вычислении векторных полей Е и Н, а ограничимся только поиском значений электрического тока поляризации J.
Аналитическая подготовка к вычислениям
Для использования в численных расчетах приведенные формулы определения коэффициентов следует преобразовать к приемлемому виду, а именно избавиться от интегралов.
Будем проводить дискретизацию с равномерным разбиением вдоль координатных осей прямоугольной декартовой системы координат и возьмем равный шаг дискретизации по всем направлениям. В результате вся область Q будет состоять из N TV, х N2 х Ny элементарных ячеек объемом П , каждая из которых представляет собой куб со стороной h. Значение (4.6) для малых кубов П( можем взять приближенно равным произведению значения подынтегральной функции в узловой точке на объем куба:
