ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2015, том 55, № 6, с. 933-946
УДК 519.642.8
О НОВЫХ РЕАЛИЗАЦИЯХ 2-ФАКТОР-МЕТОДА1)
© 2015 г. А. Ф. Измаилов
(119992 Москва, Ленинские горы, МГУ) e-mail: izmaf@ccas.ru Поступила в редакцию 25.09.2014 г.
Предлагаются новые способы реализации так называемого 2-фактор-метода, предназначенного для поиска особых решений нелинейных уравнений. Известные до сегодняшнего дня варианты метода имели весьма трудоемкую итерацию, реализация которой требовала вычисления сингулярного разложения производной решаемого уравнения. Новая экономичная реализация основана на методе Гаусса с выбором главного элемента. Кроме того, в работе рассматриваются возможности глобализации сходимости метода. В совокупности рассматриваемые средства превращают принципиальную схему 2-фактор-метода в действительно практический алгоритм. Библ. 16. Табл. 2.
Ключевые слова: нелинейное уравнение, особое решение, ньютоновские методы, идентификация структуры особенности.
DOI: 10.7868/S004446691506006X
1. ВВЕДЕНИЕ
Будем рассматривать систему нелинейных уравнений
F(x) = 0 (1.1)
с дифференцируемым отображением F: —»- R". Решение X е системы (1.1) называется особым, если
detF'(X) = 0;
в противном случае решение X называется регулярным. Эффективные методы отыскания особых решений — важная проблема численного анализа, интерес к которой объясняется как ее математической нетривиальностью, так и наличием важных приложений в теории бифуркаций (см., например, обзор в [1]), в оптимизации, включая комплементарные и связанные с ними задачи (см. [2], [3]), и в других областях. Известно, что классический метод Ньютона при сходимости к особым решениям утрачивает присущую ему сверхлинейную скорость сходимости, не говоря уже о том, что его итерации могут не быть определены в сколь угодно малой окрестности особого решения.
Информацию о существующих подходах к поиску особых решений нелинейных уравнений можно найти в [4]—[6]. В данной работе речь пойдет о подходе, который развивался в [7], [8], [9, глава 3], а идея которого восходит к так называемому 2-фактор-методу (см. [10], [11]). При этом основное внимание будет уделено следующим двум вопросам: экономичной практической реализации данного подхода, при которой каждая итерация 2-фактор-метода не требовала бы недопустимо больших вычислительных затрат, а также глобализации сходимости этого метода. Известные на сегодняшний день реализации рассматриваемого подхода используют весьма затратные средства идентификации структуры особенности в решении, и в результате сверхлинейная скорость сходимости итераций метода часто не компенсирует высокой дополнительной (по сравнению с методом Ньютона) стоимости каждой итерации. Что же касается глобализации сходимости, данный круг вопросов для методов поиска особых решений нелинейных уравнений вообще не обсуждался в известной автору литературе.
1) Работа выполнена при финансовой поддержке РФФИ (код проекта 14-01-00113) и CNPq (Бразилия; грант PVE 401119/2014-9).
2. РЕГУЛЯРИЗОВАННОЕ УРАВНЕНИЕ И 2-ФАКТОР-МЕТОД
Коротко опишем подход, разработанный в [5], [7], [8], [9, глава 3]. Всюду далее считаем выполненным следующее предположение:
(F) отображение Fдважды дифференцируемо в окрестности точки X, причем его вторая производная непрерывна в точке X.
Иногда будет использоваться более сильное предположение:
(F*) отображение F дважды дифференцируемое в окрестности точки X, причем его вторая производная липшицева относительно точки X , т.е.
F'(x) - F' (X) = O01 x - XI)
при x —- X.
Рассматриваемый подход подразумевает, что определены отображения P : R" —»- R" х n и p : R" —- R", удовлетворяющие следующим требованиям:
(P1) P = Р( X) удовлетворяет включению imF'(X) с ker P;
(P2) p = p( X) принадлежит ker F'(X);
(РЗ) отображения P и p непрерывны в точке X, причем
IIP(X) - Pllllp(X) -p\\ = o(||x - XI)
при x —► X;
(P4) отображение F2-регулярно в точке X по направлению p , т.е.
det( F' (X) + PF " (X )[p ]) = 0.
Наряду с (Р3) будет также использоваться более сильное требование
(Р3*) ||P(x) - P || = O(||x - X ||), |P(x) - p || = O(||x - X ||) при x — X.
Введем отображение Ф : R" —»- R",
Ф( x) = F(x) + P( x) F' (x) p (x). (2.1)
Идея рассматриваемой регуляризующей конструкции состоит в замене уравнения (1.1) уравнением
Ф( x) = 0. (2.2)
Основные свойства этой регуляризации суммированы в следующей теореме.
Теорема 1. Пусть X е R" - решение уравнения (1.1), причем выполнены предположения (F) и (Р1)-(Р3).
Тогда X является решением уравнения (2.2) с введенным в (2.1) отображением Ф, и это отображение дифференцируемо в точке X, причем
Ф' (X) = F (X) + PF " (X )[p ]. (2.3)
Соответственно, если дополнительно выполнено предположение (Р4), то detФ'(X) Ф 0, т.е. X -регулярное решение уравнения (2.2).
Кроме того, если вместо (F) и (Р3) выполнены предположения (F*) и (Р3*) соответственно, то имеет место оценка
Ф(x) - Ф^) - Ф'(X)(x - X) = O(||x - X||2)
при x —- X.
Заметим, что отображения P и p здесь не предполагаются гладкими и даже непрерывными в окрестности точки X. Соответственно, отображение Ф может не быть гладким и даже непрерывным в точках x Ф X, сколь угодно близких к X. Поэтому в сделанных предположениях стандартный метод Ньютона к регуляризованному уравнению (2.2) не применим. С другой стороны, фор-
мула (2.3) дает информацию о структуре производной Ф в точке X, что подсказывает идею специальных итерационных методов решения уравнения (2.2), использующих прямую аппроксимацию Ф'( X). Самый очевидный из таких методов, называемый 2-фактор-методом, задается итерационной формулой
Xk +1 = Xk- (f(Xk) + P(Xk)F"(Xk)[p(Xk)])-1Ф(Xk), k =0, 1, .... (2.4)
Из теоремы 2 легко выводится следующий результат о локальной сходимости и скорости сходимости данного метода.
Теорема 2. Пусть X е I — решение уравнения (1.1), причем выполнены предположения (F) и (Р1)—(Р4).
Тогда для любого начального приближения x0, достаточно близкого к X, формула (2.4) корректно определяет последовательность {xk}, которая сходится к X со сверхлинейной скоростью.
Если же вместо (F) и (Р3) выполнены предположения (F*) и (Р3*) соответственно, то скорость сходимости квадратичная.
Приведенная регуляризующая конструкция весьма элегантна, но для ее практического использования необходимо указать конструктивные способы задания отображений P ир, удовлетворяющих предположениям (Р1)—(Р4). Способы такого рода, предлагавшиеся до сих пор, были основаны либо на знании величины rank F'(X) (см. [2]), либо на дорогостоящем вычислении сингулярного разложения оператора F'(x) в имеющемся приближении x е I к решению X (см. [9, п. 3.1.3]). B следующем разделе предлагается значительно более практичный способ задания этих отображений, использование которого, соответственно, делает значительно более практичными и 2-фактор-метод, и саму рассматриваемую регуляризующую конструкцию, а также более общие регуляризующие конструкции из [4], [6], [12]. Данный способ использует разработанный в [13] алгоритм аппроксимации так называемого подпространства вырожденности, который, в свою очередь, в своей основе имеет не что иное, как общеизвестный метод Гаусса с выбором главного элемента (см., например, [14, § 28]).
3. ОТОБРАЖЕНИЯ P и р
Рассмотрим следующую задачу. Пусть A е R" х " — некоторая матрица ранга r е {0, 1, ..., "}. Необходимо построить такое отображение U некоторой окрестности A в R"х" в пространство I х (" -r), чтобы для любой предельной точки U отображения U(-) в точке A были выполнены соотношения rank U = " — r и AT U = 0. Ранг r известным не предполагается, и его определение — часть задачи.
Следующий алгоритм, который для заданной матрицы A е R" х " и известной величины t = = t(A) > 0, некоторым образом оценивающей расстояние от A до A, вычисляет матрицу U(A) е е R" х (" - r) таким образом, чтобы получаемое отображение U(-) обладало требуемыми свойствами, был предложен в [13, алгоритм 2].
Алгоритм 1
Полагаем = 0 , и0) = I е К" х п, А(0) = А, Я(0)(А, 0 = А, 5 = 0. 1. Если 5 = п или
|ДМ(А, ОН ^ t,
полагаем г = 5 и переходим к шагу 4. В противном случае находим индексы 4 е {1, ..., п}\$5 ие е {1, ..., п}, для которых справедливо
Кл = тах К1-
I е { 1,..., п}\ У е { 1,-, п}
Полагаем +1 =
2. Определяем матрицу преобразования Ts) е
Ts) =
1 ij
A(s)
-j, i = ^ j Й +1,
hJs i,j = 1, n.
1, i = j,
0, в противном случае,
Полагаем Us + 1) = Us)Ts), A(s +1) = (7(s))TA(s) и определяем R(s + 1)(A, t) как матрицу, составленную из столбцов A(s + 1) с номерами i е {1, ..., n}\$s + 1.
3. Увеличиваем s на 1 и переходим к шагу 1.
4. Определяем U(A, t) как матрицу, составленную из столбцов Ur) с номерами i е {1, ..., n}\$r. Возвращаясь к определению отображения P, предположим, что известна функция а : R" —»- R такая, что она непрерывна в точке X, а(X) = 0, и имеет место оценка
||x - XI = O (|а( x)|)
при x —- X. Если, например, отображение Fдважды дифференцируемо в точке X, и выполняется равенство
{^ е R"|f(X)£, = 0, F'(X)[^,^]e imF'(X)} = {0}, (3.1)
то нужными свойствами обладает функция a(x) = ||F(x)||1/2 (см. [9, следствие 1.3.1, замечание 1.3.8]). Определим отображение P следующим образом:
P( x) = U( x)(( U(x))TU(x))-1 ( U( x ))T, (3.2)
где матрица U(x) = U(F '(x), ||F(x)||)1/2 определяется с помощью алгоритма 1, применяемого к A = F '(x) и t = ||F(x)||1/2.
Нетрудно убедиться, что если матрица U е Rn х (n - r) удовлетворяет соотношениям rank U = n - r и (F'(x))T U = 0, то ортогональный проектор P на (imF'(X ))1 имеет вид
P = U( UT U)U
(вне зависимости от конкретной реализации U с указанными свойствами). Поэтому из [13, предложение 2] и предположений на функцию а следует, что отображение P, определяемое согласно (3.2), непрерывно в точке X, причем P(X) = P, а значит, выполняется (Р1). Более того, как показано в [13, разд. 4], это отображение является липшицевым относительно X, т.е. выполняется первая оценка в (Р3*).
Далее, чтобы задать отображение p, определим отображение П : R" —► Rn х ",
П( x) = V( x)(( V(x))TV(x))-1 ( V( x ))T, (3.3)
где матрица V(x) = U((F '(x))T, ||
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.