ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2007, том 47, < 11, с. 1819-1829
УДК 519.614
О НЕКОТОРЫХ ДВУСТОРОННИХ АНАЛОГАХ МЕТОДА НЬЮТОНА РЕШЕНИЯ НЕЛИНЕЙНОЙ СПЕКТРАЛЬНОЙ ЗАДАЧИ
© 2007 г. Б. М. Подлевский
(79000 Львов, ул. Наукова, 3-6, ИППММ НАН Украины) e-mail: podlev@iapmm.lviv.ua
Поступила в редакцию 20.03.2007 г.
Переработанный вариант 01.06.2007 г.
Рассматриваются итерационные алгоритмы нахождения двусторонних приближений к собственным значениям нелинейных спектральных задач для матриц, которые используют эффективную численную процедуру вычисления первой и второй производных от детерминанта матрицы. Рассматриваются также вычислительные аспекты использования этой процедуры в алгоритме отыскания всех собственных значений нелинейной спектральной задачи из заданной области комплексной плоскости изменения спектрального параметра. На модельных задачах иллюстрируется эффективность алгоритмов. Библ. 14. Табл. 4.
Ключевые слова: нелинейные алгебраические спектральные задачи, итерационный алгоритм, двусторонний аналог метода Ньютона.
ВВЕДЕНИЕ
В работе предлагается численный метод решения следующей спектральной задачи. Пусть D(X) - заданная матрица размерности n х n, нелинейно зависящая от спектрального параметра X. Требуется найти те значения параметра X е C, называемые собственными значениями, при которых существуют нетривиальные решения x, y е Cn уравнений
x *D(X) = 0, D(X) y = 0, (1)
где индекс звездочка означает сопряженное транспонирование. Искомые значения X одни и те же для обеих задач (1) и являются решениями уравнения
f(X) = detD (X) = 0. (2)
Далее предположим, что элементы матрицы D(X) - достаточно гладкие функции скалярного параметра X (дважды непрерывно дифференцируемые в действительном случае или аналитические в комплексном) в некоторой области.
В отличие от существующих численных методов решения нелинейных спектральных задач (см., например, [1]-[5]), в данной работе для нелинейной спектральной задачи (1) мы предлагаем итерационный процесс двустороннего аналога метода Ньютона отыскания простого действительного собственного значения как корня соответствующего нелинейного скалярного уравнения (2), не раскрывая определителя. Это означает, что левая часть уравнения (2) в явном виде не задается, а предлагается алгоритм нахождения функций f(X),f(X) и/1'^) при фиксированном значении параметра X, используя для этого LU-разложение матрицы D(X). Кроме того, предложенный алгоритм позволяет, используя принцип аргумента аналитической функции, определить как число собственных значений, расположенных в заданной области G комплексной X-плоскости, так и найти начальные приближения ко всем собственным значениям задачи (1) в области G с последующим их уточнением одним из известных итерационных методов, в частности и двусторонним аналогом метода Ньютона.
1. АЛГОРИТМ ВЫЧИСЛЕНИЯ f(X), f(X) И f"(X)
Известно, что матрица D(X) при любом фиксированном X может быть представлена в виде
D (X) = L (X) U (X), (3)
где L(X) - нижняя треугольная матрица с единичными диагональными элементами, а U(X) - верх-
1819
няя треугольная матрица. Тогда
п
/(X) = ае1Ь (X) аеш (X) = П ^(А) •
г = 1
Поскольку элементы квадратной матрицы D(X) (а следовательно, и и(А)) являются дифференцируемыми функциями X, то для любых X получаем, что
пп
/ (X) = X V ГГ(А)П МА),
г=1 г =1
г ф г
/" (X) = X ^(А)П и"(А) + X V кк(А)
к = 1
г = 1 г ф к
к = 1
X V^П uгг(X)
1 = 1 1 ф к
г = 1 г ф к
г ф1
где V,,(X) = иг'г (X), а (X) = VУи (X). Для определения значений vii(X) продифференцируем (3) по X. Получим
В( X) = М(А)и(А) + L(X)V(X), (4)
где
B(X) = ад), M(X) = L'(X), V(X) = U'(X),
а V(X) - элементы матрицы V(А). Теперь, дифференцируя (4) по X, получаем
ОД = N(А)U(А) + 2M(А)V(А) + L(X)W(X), (5)
где
С(А) = В'(А), N(X) = М'(А), W(X) = V'(X), а ^Й(А) - элементы матрицы W(X).
Таким образом, для вычисления/(X),/'(X) и/"(X) необходимо при фиксированном X = Ат вы-
числить
D = LU, В = Ми + LV, С = NU + 2MV + LW,
(6)
откуда
пп
/(Ат) = П иг г , ^ (Ат) = X Vгг П иг '
г =1 г = 1 гФ г
пп
/'(Ат) = £ №ккП и г+ X V
кк
к =1 г = 1 к =1
г ф к
пп
X V11 П и г
1 = 1 г= 1
1 ф к г Ф к
гФ 1 У
(7)
Для вычисления элементов матриц в разложении (6) могут использоваться рекуррентные соотношения
г = 1, 2, ..., п,
п
п
п
п
г = 1
п
игк =
г - 1
1гк = ыгк - X к = г, _, п,
1=1
/ г -1 Ч
=
гг - X 1
1=1
г]и]г
игг, г = г + 1, ..., п•
Vгк = Ьгк - X (тг1и1к + 1г] V]k), к = г> n,
1=1
тгг =
Ьгг - X (+ VVjг) - 1гVг
1=1
игг, г = г + 1,___, п.
п=
™гк =
г-1
Сгк - X (пг]'и}-к + 2т^ 1к + /^д), к = г, _, п,
1=1
Сгг - X (пг1и1г + 2тгУ ^г + А^Л - 2тгVгг - 1гг™г
1=1
г = г + 1, _, п.
Нужно отметить, что этот алгоритм может оказаться неустойчивым и даже недействительным, если некоторый элемент игг = 0. Чтобы исключить такие случаи, применяют при Ци-разло-жении матрицы ряд перестановок строк (и/или столбцов) матрицы D с одновременным выбором главного элемента, как в методе Гаусса. В этом случае разложение (6) запишется в виде
PD = LU,
РВ = ми + ЦУ,
РС = NU + 2МУ + LW,
где Р - матрица перестановок, причем detP = (-1)?, где q - число перестановок (например, строк). В таком случае соотношения (7) примут вид
/(Ат ) = (-1 )<? П иг, / (Ат ) = (-1 ^ X V кк П
к= 1 г = 1 г ф к
пп
/"(Ат) = (-1 Г X ™кк П игг + (-1 Г X V кк
к =1 г = 1 г ф к
к = 1
пп
X VЛ П иг
1=1 г=1 1 ф к г ф к
г ф1 у
Далее в работе предлагается использовать описанный алгоритм вычисления производных на основе Ци-разложения матрицы D(А).
г-1
г-1
г - 1
п
п
п
г=1
п
2. ВЫЧИСЛИТЕЛЬНЫЕ АСПЕКТЫ АЛГОРИТМА ОТЫСКАНИЯ НАЧАЛЬНЫХ ПРИБЛИЖЕНИЙ СОБСТВЕННЫХ ЗНАЧЕНИЙ
Принцип аргумента неоднократно применялся при решении различных задач, связанных с определением количества собственных значений, лежащих в заданной области (см., например, [4], [6]-[8]). Так, в [6] был предложен алгоритм нахождения всех нулей аналитической функции в заданной области на комплексной плоскости. В [4], [7] и [8] были предложены модификации этого метода и его обобщение для решения нелинейных задач на собственные значения. В настоящей работе, в частности, также рассматриваются некоторые вычислительные аспекты его использования.
Пусть характеристическая функция (2), которая по условию является аналитической, имеет в О (с учетом кратности) т нулей Х1, ..., Xm и не имеет нулей на границе Г области О. Тогда, как известно, число нулей т функции/, лежащих в О, определяется в соответствии с принципом аргумента формулой
т - 'о - 2Ыжл (8)
г
Введем обозначения
" -шSXk/mdX, к 2--; (9)
г
тогда (см. [9, с. 151]) имеют место соотношения
т
^¡Xk] - 'к, к - 1, 2, ..., т. (10)
1 -1
Таким образом, зная т и 'к, к = 1, 2, ..., т (см. (8), (9)), из системы (10) можно найти нули функции (2), т.е. все собственные значения задачи (1), лежащие в заданной области О.
Остановимся подробнее на вычислительных аспектах алгоритма, а именно на той его части, где вычисляются величины 'к, к = 1, 2, ..., т.
В [4] предлагалось в качестве О выбирать круг и использовать квадратурную формулу прямоугольников для нахождения 'к. В [8] в качестве О предлагалось выбирать круг или прямоугольник, а подынтегральную функцию ф(Х) = /'(X)//(X) аппроксимировать линейным или квадратичным многочленом на каждом отрезке контура Г, предварительно разбив его на части. При любом способе вычисления величин 'к необходимо вычислять значения функции ф(Х) в заданных точках. Для этого в [6] предлагалось использовать LU-разложение матрицы D(X) для определения функции /(X) в заданной точке, а производную определять при помощи соотношений
т -1
/ '(Х)~£ 1а} (X - X о)1 -1,
1 -о
где
а, - 1 Г -^^Л,
1 ^'га - x„)'*l
или воспользоваться теоремой о следе матрицы и заменить отношение /XX)//(X) на Тгф-1^^^)}, а нахождение Тгф-1^^^)} свести к решению систем линейных уравнений (см. [8]). Здесь предлагается использовать LU-разложения матрицы D(X) для вычисления как /(X), так и /XX).
Не ограничивая общности, рассмотрим в качестве области О, ограниченной контуром Г, круг О^*, р) с центром в точке X* и радиусом р. Используя замену X(t) = X* + р ехр(2л'0, преобразуем интеграл (9) к виду
1
'к - ]Х(0'рехр(2пи)^/■(XX(t))Л. (11)
о
Разобьем отрезок [0, 1] на N равных частей и заменим интеграл (11) какой-нибудь квадратурной формулой (например, прямоугольников, как в [4]). Получим
'к- N )'р ехр (' ^/ш - (12)
1 -1 1
где X = X* + р ехр('N1, 1 = 1, 2, ., N.
Таким образом, в (12) требуется вычислить значение функции /(X) и ее производной только на границе области О. Это можно осуществить с помощью разложения (6). Тогда отношение /(Ау)//(А,), используя представления (7), приведем к виду
/ (X1)//(X1) = X ^
(13)
г = 1
Таким образом, учитывая (13), для вычисления величин sk, к = 0, 1, ..., т, получаем соотноше-
ния
N
^ =
¿X <А/Рехр(.т£ V*
1=1
:2Ш. N
иг
(14)
Следовательно, при помощи соотношений (14) можно вычислить количество т = 50 собственных значений, находящихся в области О, а также правую часть системы (10), для решения которой можно применить, как в [4], метод Ньютона. В результате получим некоторые, вообще говоря, грубые приближения ко всем собственным значениям задачи (1), лежащим в области О.
п
гг
3. ДВУСТОРОННИЕ АНАЛОГИ МЕТОДА НЬЮТОНА
Для уточнения некоторых грубых приближений к собственным значениям, полученным по тем или иным соображениям, в частности и при помощи алгоритма, описанного выше, предлагается использовать итерационные алгоритмы двусторонних приближений. Они обладают тем важным преимуществом по сравнению с другими приближенными методами, что дают возможность на каждом шаге итерационного процесса искомые решения заключать в "вилку" и тем самым получать удобную апостериорную оценку погрешности последовательных приближений.
Пусть известно приближение А0 к собственному значению X*, достаточное для применения двустороннего аналога метода Ньютона (см. [10])
А = А__/ ( А 2 т )/'( А2т )
А2 т +1 А2т . 2 - - ,
Г (А2т) - /(А2т) У" (А2 т) п л
т = 0, 1, _ , (15)
X _ X /( А2 т + 1 )
А2т +2 = А2т +1- /-1 / -
I (А2т+1)
для п
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.