CC BY
16
19. Felippa С. A. Optimization of finite element grids by direct energy search // Appl. Math. Model. 1976. 1. N 1. P. 93-96.
20. Delfour M.,Payre G.,Zolesio J.-P. An optimal triangulation for second-order elliptic problems// Comput. Meth. Appl. Mech. Engrg. 1985. 50. N 3. P. 231-261.
21. To u rig n у Y., Hiilsemann F. A new moving mesh algorithm for the finite element solution of variational problems // SIAM J. Numer. Anal. 1998. 35. N 4. P. 1416-1438.
22. Babuska I., Aziz A. K. On the angle condition in the finite element method // SIAM J. Numer. Anal. 1976. 13. N 2. P. 214-226.
23. Fried I. Condition of finite element matrices generated from nonuniform meshes // AIAA J. 1972. 10. N 2. P. 219-221.
24. Verfiirth R. A posteriori error estimators for the Stokes equation non-conforming discretization // Numer. Math. 1991. 60. N 2. P. 235-249.
25. Rivara M.-C. Mesh refinement processes based on the generalized bisection of simplices // SIAM J. Numer. Anal. 1984. 21. N 3. P. 604-613.
Поступила в редакцию 30.03.05
УДК 519.6
M. Дана, Х. Д. Икрамов
О РЕШЕНИИ СИСТЕМ ЛИНЕЙНЫХ УРАВНЕНИЙ, МАТРИЦЫ КОТОРЫХ ЯВЛЯЮТСЯ МАЛОРАНГОВЫМИ ВОЗМУЩЕНИЯМИ ЭРМИТОВЫХ МАТРИЦ
(кафедра общей математики факультета ВМиК)
1. Введение. Пусть нужно решить систему линейных уравнений
Ах = Ь (1)
с нормальной га X га-матрицей А. Нормальность в данном контексте означает выполнение условия
АА* = А* А. (2)
Для разреженной неэрмитовой матрицы А стандартным выбором алгоритма для решения системы (1) является обобщенный метод минимальных невязок GMRES (см. [1, раздел 6.6.6]). Однако в этом методе не удается использовать специфику матрицы А, выражаемую равенством (2): алгоритм Арнольди, лежащий в основе GMRES, приводит нормальную матрицу А к той же форме Хессенберга, что и матрицу общего вида. В результате сохраняются два главных недостатка, присущие GMRES по сравнению с методами типа Ланцоша: большой расход памяти и зависимость количества арифметической работы, выполняемой на отдельном шаге, от его номера к.
При дополнительном предположении, что спектр матрицы А принадлежит алгебраической кривой Г невысокого порядка (примером реалистической ситуации, где это предположение выполнено, могут служить нормальные матрицы с почти-вещественным спектром), в [2] предложен иной подход к решению системы (1), который вкратце состоит в следующем. Методы типа Ланцоша и GMRES ищут приближенные решения в крыловских подпространствах возрастающей размерности, т.е. в линейных оболочках конечных отрезков степенной последовательности
х, Вх, В2х,... . (3)
Здесь В — это или сама матрица А, или некоторая матрица, тесно с А связанная; в качестве х обычно (при отсутствии хорошего начального приближения) берут правую часть Ь системы (1) или полученной из нее "предобусловленной" системы.
Метод, предложенный в [2], является алгоритмом минимальных невязок для обобщенных крылов-ских подпространств, натянутых на конечные отрезки последовательности
т 4т А*т А^т А А*т А*^т А^т А^А*т (А\
В отличие от СМНЕЭ этот метод, называемый MINRES-N, обладает фиксированной длиной рекурсии, зависящей от порядка то кривой Г. Так, для случая то = 2, наиболее подробно рассмотренного в [1], длина рекурсии равна 6. То же самое можно сказать иначе: вместо унитарного приведения матрицы А к форме Хессенберга, неявно выполняемого в СМНЕЭ, А теперь приводится к специальной компактной форме, которую можно считать блочно-трехдиагональной или ленточной матрицей с шириной ленты, зависящей от то. Тем самым устраняются оба отмеченных выше недостатка СМНЕЭ.
В настоящей статье мы исследуем поведение метода МШНЕЭ^ для двух классов матриц А, которые могут рассматриваться как малоранговые возмущения эрмитовых матриц. Напомним вначале, что всякая квадратная комплексная матрица А может быть единственным образом представлена в виде
А = В + 1С, В = В*, С = С*. (5)
Эрмитовы матрицы
В=1-(А + А*), С=±(А-А*) (6)
называются соответственно вещественной и мнимой частями матрицы А. Нормальные матрицы А могут быть охарактеризованы тем свойством, что соответствующие им матрицы (6) перестановочны. В разделе 2 рассматривается класс нормальных матриц А, выделяемый условием
к = гапкС <С п, (7)
где п — порядок А. Спектр такой матрицы А принадлежит объединению вещественной оси и самое большее к прямых, параллельных мнимой оси, т.е. вырожденной алгебраической кривой порядка не выше к1. Следовательно, метод МШНЕЭ^ применим к матрицам этого типа. Мы показываем, что для матриц со свойством (7) МШНЕЭ^ имеет следующую привлекательную особенность: начиная с некоторого шага рекурсия длины 3(к + 1), обычно выполняемая этим методом, может быть заменена трехчленной рекурсией, после чего метод фактически совпадает с эрмитовым процессом Ланцоша. Достигаемое этим путем ускорение сходимости иллюстрируется несколькими примерами. В разделе 3 рассматриваются системы с матрицами (5), для которых
гапкС =1. (8)
В отличие от раздела 2 не требуется, чтобы А была нормальной матрицей, т.е. ее вещественная и мнимая части В и С не обязаны коммутировать. Мы показываем, что к системам этого типа применим метод МШНЕЭ-^, представляющий собой специализацию метода МШНЕЭ^ для нормальных матриц со спектром, лежащим на кривой второго порядка. Возможность применения МШНЕЭ-^ к анормальным матрицам со свойством (8) тем более удивительна, что их спектры могут быть совершенно не похожи на те, для которых метод был разработан. Эффективность метода иллюстрируется несколькими примерами решения ленточных систем; на этих примерах производительность МШНЕЭ-^ сопоставляется с производительностью библиотечной программы Ма^аЬ'а, реализующей СМНЕЭ.
2. Нормальные малоранговые возмущения эрмитовых матриц. Рассмотрим матричную последовательность
/„, А, А*, А2,АА*, А*2, А3, А2 А*, АА*2,А*3, А4,..., (9)
соответствующую векторной последовательности (4). Члены последовательности (9), называемые словами из А и А*, получаются из одночленов
1, х, у, ж2, ху, У2, Ж3, х2у, ху2, у3, х4,... (10)
подстановкой
X А, у А*.
Степенью слова Аг А*3 называется число I + ]. Совокупность слов одной и той же степени то образует то-й слой последовательности (9). По определению единичная матрица /„ составляет нулевой слой.
Обозначим через Ст линейную оболочку матриц, входящих в слои последовательности (9) с номерами 0,1, 2,..., т. Положим
£т=дш1Ст (11)
и
= - Ап-1 (12)
Число ит назовем шириной т-го слоя. По определению и0 = 1. Заметим, что в блочно-трехдиагональ-ной форме, к которой МШНЕЭ^ (неявно) приводит матрицу А, размеры блоков определяются как раз числами ит.
Числа £т можно вычислить иначе, чем указывает определение (11). Используя представление (5) и вытекающие из него равенства
А* = В- ¿С, А2 = В2 - С2 + 2 ¡ВС,
(13)
АА* = В2+С2, А*2 = В2 - С2 - 21ВС, ...,
заменим (9) последовательностью
1п,В,С,В2,ВС,С2,В3,В2С,ВС2,С3,В4,... . (14)
Вместо слов из А и А* мы имеем в (14) слова из В и С. Степенью слова ВгС3 по-прежнему называется число Как и (9), последовательность (14) разбивается на слои 0,1,2,3,4,.... Из формул (13)
видно, что введенное выше подпространство Ст можно рассматривать как линейную оболочку слоев 0,1, 2,..., т последовательности (14).
Теорема 1. Пусть нормальная п X п-матрица А удовлетворяет условию (7). Тогда
ит ^ 1 для т > к. (15)
Доказательство. Наряду с будем рассматривать подпространства
Мг = эрап{В*~8С8, 0 < в ^ Ь ^ г}, ¿ = 1,2,.... (16)
Содержательно М{ есть линейная оболочка тех слов из первых г + 1 слоев последовательности (14), в которые входит хотя бы один множитель С.
Из перестановочности эрмитовых матриц В ж С следует, что все члены последовательности (14) могут быть приведены к диагональному виду одним и тем же подобием. При этом диагональные матрицы, соответствующие словам из подпространств М,, содержат нули в тех диагональных позициях, где стоят нули в диагональной форме матрицы С. Отсюда следует, что
с1 \тМг^к Мг. (17)
Покажем, что для последовательности
М0 = {0} С Мг С М2 С ... С Мг С ... (18)
справедлива импликация
Мг = Мг+1 Мг+1=Мг+2. (19)
Для этого прежде всего заметим, что подпространства, выписанные в (18), связаны рекурсией
Мг = эрап{Мг-иСМг-ъВ'^С}, г = 1,2,... .
Если М, = MiJГl, то
Отсюда
С'МгСМг И В'С'еМг-СМ{+1 = CMi С Mi
и
Вг + 1С = В(ВгС) 6 ВМг С Мг +1 = Мг.
Поэтому
Мг+2 = зрап{М^1,СМ^1,Вг+1С} С Мг Учитывая очевидные включения Mi С М{+\ С Mi+2, получаем (19).
Доказанная импликация вместе с оценкой (17) означает, что строгое возрастание размерностей в последовательности (18) возможно лишь для первых к значений индекса I. Для то, больших, чем к, приращение размерности подпространства Ст по сравнению с Ст-\ возможно лишь за счет матрицы Вт. Это доказывает утверждение (15).
Теорему 1 можно использовать для ускорения сходимости метода МШНЕЭ-М. Объясним это с помощью следующего примера.
Пусть А — диагональная матрица порядка 2000, все диагональные элементы которой, за исключением четырех, расположены в двух симметричных сегментах вещественной оси [—14, —10] и [10,14]. Остальные четыре элемента образуют две сопряженные пары ^ ±0,1 % и ^ ±0,5г. Итак, спектр сосредоточен на координатных осях, т.е. принадлежит (вырожденной) кривой второго порядка
Поэтому к системе с матрицей А можно применить метод МШНЕЭ-^. То обстоятельство, что диагональная система уравнений решается тривиально, не должно вызывать недоразумений. Действительно, всякая нормальная матрица А приводится к диагональному виду И в некотором ортонормированием базисе {е}, а метод МШНЕЭ-^, проводимый с матрицей А, изометричен такому же процессу, применяемому к И в базисе {е}. Поэтому с точки зрения исследования сходимости метода ограничение диагональными матрицами не влечет за собой никакой потери общности.
Правую часть Ь системы (1) мы выбирали в виде
Таким образом, точным решением системы является вектор, составленный из единиц.
На рис. 1, а показаны графики евклидовых длин невязок для методов МШНЕЭ-^ и СМНЕЭ, применяемых к описанной выше системе. В этих расчетах СМНЕЭ был представлен библиотечной функцией gmres системы Ма^аЬ, а МШНЕЭ-^ — составленной нами процедурой Ма^аЬ. Для обоих методов использовался один и тот же критерий остановки:
где г — текущая невязка, а положительное число е задается пользователем. Для обсуждаемого примера
Результаты сравнения двух методов для данной системы можно считать благоприятными для МЩЦЕЭ^2: при приблизительно одинаковом числе итераций (50 против 47 итераций для СМНЕЭ) выигрыш в общем времени счета происходит вследствие меньшего количества выполняемой арифметической работы. Однако обращает на себя внимание тот факт, что график ||г|| для МШНЕЭ-^ состоит из ступенек ширины 4. Это означает, что после спуска на очередную ступеньку три последующие итерации не уменьшают норму невязки.
Объяснение состоит в следующем: МШНЕЭ-^ исходит из предположения, что все числа ит равны двум. Между тем согласно теореме 1 для рассматриваемой системы ит ^ 1 при то > 4. Пусть ^1,^2) • • • — векторы ортонормированной последовательности, которую строит МЩЦЕЭ^2. (Они называются обобщенными векторами Ланцоша.) При этом в качестве д! берется нормированный вектор правой части Ь, а каждая последующая пара векторов Ланцоша теоретически получается посредством ортонормализации, применяемой к векторам очередного слоя последовательности (4). В частности, векторы образуют ортонормированный базис подпространства £4. После того как на 10-м
шаге процесса построен вектор дю, МШНЕЭ-^ пытается вычислить вектор дц, оставаясь в том же (пятом) слое последовательности (4). Однако = 1, т.е. в пятом слое нет вектора, независимого с дю и ортогонального к £4. В методе МШНЕЭ-^ это обстоятельство проявляется в том, что вектор г, полученный из второго вектора пятого слоя ортогонализацией к предыдущим векторам Ланцоша, имеет малые компоненты, фактически являющиеся накопленными ошибками округлений. Несмотря на это, МШНЕЭ-^ нормирует г и принимает его в качестве дц. Естественно, что такой вектор, случайный по отношению к подпространствам Ст, мало помогает в уменьшении длины невязки. Описанная ситуация повторяется на каждом последующем слое.
Исправить отмеченный недостаток можно, вводя в процесс ортогонализации проверку
ху = 0.
Ъ=Ае, где е= (1,1,..., 1)т.
(20)
(21)
£ = Ю"8.
Рис. 1
выполняемую перед нормированием, дающим очередной вектор Ланцоша. Число 6 задается пользователем; в описываемых экспериментах мы полагали 5 = Ю-4. Если неравенство (22) выполнено, то текущий шаг завершается стандартным образом; в противном случае уменьшаем значение и ширины текущего слоя на единицу, отбрасываем вектор г и переходим к вычислениям, соответствующим следующему шагу процесса ортогонализации.
Рис. 1, б показывает поведение норм невязок для СМИЕв и модифицированного алгоритма МШИЕЭ-^. Эти графики соответствуют системе, построенной как ранее: матрица А — диагональная со спектром, сосредоточенным на сегментах [—14, —10] и [10,14], кроме двух сопряженных пар ~ ±0,5г и рй ±0,9г. Некоторое отличие этой матрицы от предыдущей связано с тем, что при заданном характере спектра (заданное число чисто мнимых сопряженных пар + равномерное распределение прочих собственных значений на заданных сегментах вещественной оси) конкретные значения диагональных элементов выбирались с помощью датчика псевдослучайных чисел. Мы видим, что картина сходимости для СМИЕЭ почти полностью сохранилась; в то же время модифицированный МЩЦЕЭ^2 сходится значительно быстрее первоначального метода: критерий (21) теперь выполняется (при е = Ю-8) после 30 итераций. Время счета в три раза меньше, чем для СМИЕЭ; ступеньки ширины четыре сменились ступеньками ширины два.
Почему ступеньки лишь уменьшили ширину, а не исчезли вовсе, — авторам в настоящее время не вполне ясно. Возможно, это связано с симметричным расположением спектра рассмотренных выше матриц относительно нуля. Во всяком случае, ступенек нет (после девятого шага ортогонализации) на рис. 1, в. Показанные здесь графики норм невязок соответствуют матрице А с двумя сопряженными парами собственных значений ^ ±0,2 % И ^ ±0,4г и остальным спектром, равномерно распределенным на сегменте [1,4]. Для выполнения критерия (21) с тем же е = Ю-8 здесь достаточно 24 итераций. Заметим, что и СМИЕЭ сходится для этой системы быстрее (32 итерации), чем в предыдущих случаях.
*\\Ъ-Ахт\\
♦ |Ь — Лж«тгев||
# >№ V »« « • >0»< ♦ «к >*« \ "Ч. V "м % > >«•< * >«, V *"*- -
\ " ш ж«< ***
О 40 80 120
Номер итерации
в
Номер итерации
г
Рис. 2
3. Одноранговые возмущения эрмитовых матриц. Если матрица А не является нормальной, то вместо (4) нужно рассматривать последовательность
х^ ¿4.xу А. х^ А. Ху Ху А. ¿Ах^ А. х^ А. х^ А. А. х^....
В к-м слое этой последовательности содержится 2к векторов, а не к + 1 векторов, как в (4).
Матричные множители, стоящие перед вектором х в (23), по-прежнему можно получать из одночленов от ж и у подстановкой х —> А, у —> А*; однако теперь переменные ж и у не коммутируют.
В данном разделе символ Ст обозначает линейную оболочку векторов (а не матриц, как в разделе 2), входящих в слои последовательности (23) с номерами 0,1, 2,..., т. При этом условии символы £т и ит сохраняют свои определения (11) и (12).
Как и в разделе 2, при построении подпространств Ст последовательность (23) можно заменить последовательностью
X, Вх, Сх, В2х, ВСх, СВх, С2х, В3х, В2Сх, ВСВх, СВ2х,.... (24)
Теорема 2. Пусть квадратная матрица А = В + 1С удовлетворяет условию (8). Тогда
ьзт ^ 2, т = 1,2,________(25)
Доказательство. При т = 1 утверждение (25) очевидно: первый слой последовательности (23) (как и последовательности (24)) состоит из двух векторов.
Пусть неравенство (25) выполнено для т= 1,2,..., г. Будем различать два возможных случая: 1) в слоях 1,2,..., г последовательности (24) найдется вектор, дающий приращение размерности соответствующего подпространства Ст, причем в записи этого вектора в виде w(B,C)x самым левым множителем слова w является С; 2) такого вектора нет.
В первом случае при переходе от С{ к £¿+1 слова, начинающиеся с С, не могут дать приращения размерности, поскольку rank С = 1. Новые линейно независимые векторы можно получить лишь из независимых векторов, найденных в г-м слое, умножая последние на матрицу В. Таким образом, в этом случае o^+i ^ Ui ^ 2.
Во втором случае слои 1,2,..., г имеют ширину 1. Из единственного независимого вектора z, найденного в г-м слое, можно получить самое большее два новых независимых вектора Bz и Cz в слое г + 1. Итак, и в этом случае o^+i ^ 2.
Следствие. Квадратная матрица А, удовлетворяющая условию (8), может быть приведена посредством унитарного подобия к блочно-трехдиагональной форме, где порядок каждого диагонального блока не превосходит двух.
Из теоремы 2 и ее следствия вытекает, что метод MINRES-N2 может быть применен к решению систем с анормальными матрицами коэффициентов, удовлетворяющими условию (8). На рис. 2, а, б, е, г мы показываем графики норм невязок этого метода (и аналогичные графики для GMRES) для ленточных систем порядка 2000 с шириной ленты т = 3,5, 7, 9. Трехдиагональная матрица А, определяющая первую из этих систем, была получена так: диагональные элементы ее вещественной части В суть псевдослучайные числа, равномерно распределенные на паре сегментов [-350,-300] и [275,350], а элементы первой наддиагонали (и первой поддиагонали) равномерно распределены на паре сегментов [—70, —40] и [65,80]. Мнимая часть С имеет единственный ненулевой элемент с псевдослучайным значением из сегмента [-60,-20], стоящий в диагональной позиции (г, г) со случайным индексом г. Остальные ленточные матрицы были построены сходным образом.
Несмотря на большее число итераций, MINRES-N2 сходился значительно быстрее, чем GMRES, для всех четырех систем. Так, для т = 9 решение было найдено за 10,297 с, тогда как GMRES потребовал 34,188 с. Все вычисления проводились на персональном компьютере Pentium IV с тактовой частотой 256 Мгц.
СПИСОК ЛИТЕРАТУРЫ
1. Деммель Дж. Вычислительная линейная алгебра. Теория и приложения. М.: Мир, 2001.
2. Дана М., Зыков А. Г., Икрамов Х.Д. Метод минимальных невязок для специального класса линейных систем с нормальными матрицами коэффициентов // ЖВМиМФ. 2005. 45. № 11.
Поступила в редакцию 20.05.05
УДК 537.8
H. В. Гришина, Ю. А. Еремин, А. Г. Свешников МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ЛОКАЛЬНОГО БИОСЕНСОРА1
(лаборатория вычислительной электродинамики факультета ВМиК)
I. Введение. Развитие био- и нанотехнологий становится определяющим фактором научно-технического прогресса во многих областях. Существует настоятельная потребность в разработке эффективных средств экспресс-диагностики растворов для проведения иммунологического анализа крови, определения наличия вирусов или изменений в ДНК. При разработке современных средств
1 Работа выполнена при финансовой поддержке РФФИ-ННИО, № 04-02-04009.
