CC BY
46
Вычислительные технологии
Том 11, № 6, 2006
МЕТОД ТОЧНЫХ РЕЛАКСАЦИЙ
С. Е. МИХЕЕВ
Санкт-Петербургский государственный университет, Россия
e-mail: him2@mail.ru
Exact relaxations employing some additional information about location of the desired solution are able to improve the convergence of iterative methods, which can be presented in the simple iteration method form. Formulae for the relaxations are obtained using minimization of the maximum estimation of the error arising in the subsequent iteration.
1. Описание
Большая доля r-точечных итеративных методов поиска приближения к некоторому элементу а нормированного пространства B может быть представлена в виде
xk+1 = A(Xk), k = 0,1, 2,..., (1)
где A — базовый алгоритм, связанный с исходной задачей f из некоторого семейства F, которая входит в него неявным параметром; Xk := {xk,... , xk-r+1}; (хг}^-г — итеративные точки, принадлежащие B; Хо = {х0,... ,x1-r} — набор начальных точек, выбираемый из каких-то внешних для данного метода соображений.
Под сходимостью итераций из (1) понимается существование предела: а = xk.
При r =1 итерации (1) именуются методом простых итераций, а элемент а, если в нем непрерывен алгоритм A, — неподвижной точкой алгоритма A.
Когда алгоритм A корректно связан с исходной задачей f, элемент а (если существует) должен быть решением задачи f.
Согласно (1) алгоритм A для построения итеративной точки xk+1 использует информацию только о наборе Xk. Однако доступной и применяемой часто оказывается еще дополнительная информация Ik, в явном виде не вошедшая в список аргументов алгоритма A, а вычисляемая вспомогательным алгоритмом:
Ik+1 = Ai(Xk, Ik), k = 0 1...
В качестве примера можно привести оценку погрешности dk k-го приближения
dk > ||xk — а||, (2)
которая используется для критерия остановки итеративного процесса (dk < е), где е — заданный уровень допустимой погрешности приближенного решения. Если помимо dk не используется другой информации, то Ik = {dk}.
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.
Вместе с тем и базовый алгоритм может зависеть от некоторой, вообще говоря, иной информации. После пополнения ею /д. итеративному методу можно придать вид
Будем также считать, что результат работы базового алгоритма А не меняется при описанном выше расширении итеративной информации. Потребуем еще, чтобы все /д., к = 0,1,..., обладали единым строением 2.
Поясним сказанное примером. Если известна информация о сжимающих свойствах базового алгоритма метода простой итерации:
то из нее легко извлечь рекуррентную формулу для оценки погрешности к-й итерации
Опишем методику, позволяющую в некоторых случаях предложить на основе алгоритма А иные итерации, сходящиеся лучше, чем процесс (3). Часто объем дополнительных к базовому алгоритму вычислений для подобной модификации легко оценивается.
Принцип минимальности. Если /к не содержит оценку погрешности , введем ее туда с сохранением обозначений (/д, 2). Если начальная оценка неизвестна, назначаем
4 := +то. Положим = А(Хд, /д).
Если итерации (3) имеют теоретическое обоснование сходимости, то набору Хд, /д должно соответствовать содержащее решение а некоторое ограниченное множество
, эд,/& следующим образом. Итеративной точкой жк+1 назначить минимайзер величины (у) := тах^ехк ||у — г||, иными словами, величину
где Х£ — множество из (6). Очередная итеративная информация /д.+1 должна обладать
к = 0,1,...
(3)
Уж Є В ||А(ж) — а|| < с||ж — а||р,
(4)
(5)
а расширенному набору Хд, и^/д — суженное, вообще говоря, множество
ІД := 1(Хк,ш,4) С В.
(6)
Для построения очередного итеративного набора жк+1, /&+1 нужно использовать набор
уЄ£ гЄХк
где / обладает строением 2. Если минимайзер величины Д*(у) или очередная итеративная информация не определены единственным образом согласно вышеизложенному, то это следует сделать с привлечением каких-нибудь дополнительных соображений об их выборе.
Замечание 1. Использование принципа минимальности отменяет применение вспомогательного алгоритма А; в итеративном процессе.
Замечание 2. Принцип минимальности применим также и к итеративным процессам, не имеющим изначально отличной от X* итеративной информации (/* = 0).
Модификация алгоритма А с помощью принципа минимальности есть новый алгоритм М(Х,/, А(Х, /)). Его естественно именовать точной релаксацией (ТР) алгоритма А, ибо с имеющейся информацией он обеспечивает минимум оценки погрешности следующего приближения.
Точные релаксации порождают итеративный процесс
{жк+1Л+1} = М (X*,4, А(Х*, /*)), к = 0,1,...
Здесь начальные точки X, г =1 — г, 0, и начальная итеративная информация /° — те же, что и для исходного итеративного метода (3) 1. Для прочих индексов к = 1, 2,... итеративные точки и итеративные информации, вообще говоря, отличны от тех, что появляются в методе (3).
Итеративный процесс приближения к искомому элементу а с помощью ТР назовем методом точных релаксаций.
2. Применение принципа минимальности к методу простых итераций
Естественно ожидать, что различия в нормах и в строениях итеративной информации приведут к различиям в вычислительных формулах метода точных релаксаций любого итеративного метода и в частности метода простых итераций. Добавим конкретности.
К1. В качестве В возьмем предгильбертово (в частности, и евклидово) пространство, а под || • || будем понимать естественную норму в нем.
К2. Итеративная информация — оценка (4) и оценка текущей погрешности ^. Оценка начальной погрешности ^° > ||ж° — а|| считается заданной.
Таким образом, дополнительная информация /* состоит из постоянной части — оценки (4) и переменного числа — оценки текущей погрешности (см. соотношение (2)). Тогда, если опустить постоянный компонент итеративной информации, (1) принимает [1] вид
(ж*+1,4+1) = М(ж*, ^, А(ж*,4)), к = 0,1,...
По своему духу предлагаемая модификация М имеет отдаленное сходство с релаксациями в итеративных методах решения систем линейных уравнений (в основном в модификации методов типа метода Гаусса — Зейделя) как в “чистом виде” (см. например, [2]), так и в составе итеративных методов, требующих решения линейных систем, например, в методе Ньютона [3] или в методе сеток [4].
Когда с < 1 и р =1 либо когда Ы° 1 < 1 и р > 1, итерации (1) согласно (4) порождают сходящуюся к а последовательность {ж*}^°.
1 Если оценка начальной погрешности отсутствовала в исходном итеративном методе, то в /о она входит как с!о = +то.
Далее, где возможно, индекс текущей итерации у оценок погрешностей и итеративных точек будет опускаться, а вместо индекса последующей итерации (к + 1) будет “*”, т. е. ж := ж*, I := 4, ж* := ж*+1, 4 := 4+ь
Оказывается, условия К1 и К2 однозначно определяют очередные итеративные точку ж* и информацию 4 для точной релаксации. Покажем, как это делается.
Определим зону Аполлония и оценочный шар:
:= {а ||А(ж) — а|| < с||ж — а||р}, 5Х := {у ||у — ж|| < I}.
Поскольку решение а должно находиться как в 5Х)д, так ив 5Х и это единственное условие для формирования Х* 2, имеем Х* = Р| б^.
Согласно принципу минимальности, (к + 1)-й итеративной точкой следует назначить центр минимального шара, т. е. шара минимального радиуса, содержащего множество I*. Лучшая оценка погрешности 4 итеративной точки ж* тогда будет равна этому радиусу.
Так как минимальный шар единствен для любого множества, привлечения дополнительных соображений для устранения неединственности в принципе минимальности не требуется.
Сходимость, скорость сходимости, текущие погрешности, удобство использования метода точной релаксации естественно сравнивать с такими же характеристиками породившего его метода простой итерации.
Распространенный подход к исследованию погрешности текущей итерации ж в методе простых итераций с использованием оценки (4) заключается в получении оценки 1° погрешности начального приближения (||ж° — а|| < 1°) и последовательном применении формулы 4+1 < с1*, вытекающей из (4): 4 < с1+р+р2+---+рк 11° = (С4)рк/С, где С = р-/с. В линейном случае (р = 1 ):
4 < с4-1 < ... < с*4. (7)
Использование такого типа оценок для остановки итеративного процесса по достижении заданной точности фактически превращает метод простых итераций в метод вида (3).
Поведение оценок погрешностей в методе ТР носит совсем иной характер в силу их зависимости от последовательности результатов применения ТР на каждом шаге. То есть априори эти оценки носят вероятностный характер, и корректное сравнение их поведения следует проводить, привлекая ту или иную статистику. Вместе с тем, согласно принципу минимальности, оценки погрешностей итераций метода ТР всегда не больше оценок модифицируемого метода, а в некоторых случаях для погрешности следующей итерации по ТР возможно получить оценку сверху, которая для любой текущей итерации будет строго меньше оценки следующей итерации по исходному алгоритму.
2.1. Вывод расчетных формул для линейной сходимости (р = 1)
На шаге к по построению известно, что искомое а лежит в шаре 5^. Вычисление А(у) дает дополнительную информацию о расположении а. При фиксированном у оценка
||А(ж) — а|| < с||ж — а|| (8)
1 Добавление к двум определяющим множество Х'к условиям, оценки ||А(х) — а|| < с<1р не сузит множе-
ства Х’к, так как она из них вытекает.
описывает множество возможных расположений искомого а: 5(у, А) — известное, если В — двумерное вещественное евклидово пространство, как круг Аполлония. В пространствах произвольной размерности будем именовать его шаром Аполлония.
Теорема 1. Шар (круг) Аполлония есть действительно шар (круг), и его центр находится на прямой I, проходящей через у и А(у) [5].
Лемма. Если итерация ж не неподвижная точка алгоритма А, то ж € 5х,д.
Доказательство. Положив в (8) переменную а равной ж € 5х,д, получили бы || А(ж) — ж|| < 0, т. е. А(ж) = ж, иначе говоря, итерация ж — неподвижная точка алгоритма А. □
Из леммы, очевидно, следует, что если а — искомое решение, то
ж = а =^ 5Х,д ф 5^. (9)
Найдем диаметр минимального по величине диаметра шара 5, содержащего пересечение оценочных множеств, т. е. пересечение шара 5^ := {г ||г — ж|| < I} с шаром Аполлония
5х,д. (Оно, конечно, не пусто, если обе оценки истинны.)
В силу симметрии ситуации относительно прямой I центр шара 5 должен лежать на этой прямой. Для нахождения этого центра и радиуса (т. е. величин ж* и I*) рассмотрим содержащую I произвольную двумерную плоскость. В ней сечения исследуемых трех шаров имеют вид кругов, диаметры которых равны соответственно диаметрам этих трех шаров. В силу (9), если А(ж) = ж, то либо Бх,д С 5-, либо круги Бх,д, 5^ имеют общую хорду Q ± I. Ясно, что := Q П I является серединой (центром) хорды ^.
Проведем в круге Аполлония диаметр 5 перпендикулярно к прямой I. Затем выполним следующие действия:
— когда этот диаметр принадлежит кругу 5-, следует положить 5 = 5х,д, ибо нет меньшего круга, который бы содержал 5;
— когда 5 С 5-, найдем общую хорду Q кругов ,А и 5-.
Теорема 2. Если I) С то:
1) центры шаров 5-, Бх,д лежат на прямой I по разные стороны от Ql;
2) шар, где Q — диаметр, является минимальным шаром, содержащим |^| Бх>д.
Доказательство. Достаточно проверить утверждение в двумерном сечении.
1. Пусть точка Сд € I есть центр круга 5х,д. Параметризуем прямую I:
ССА +
Восстановим перпендикуляр к I в С(£). Он пересечет границу круга Аполлония в некоторой точке У (£). По теореме Пифагора
IIж — У (£) II2 = (У5?2)^ ) + (11Сд — ж)! + £)2 = 52/4 + (Сд — ж)2 + 211Сд — ж11£.
Следовательно, ||ж — У (£)|| монотонно увеличивается с ростом £. Но 5 С означает, что концы диаметра 5 как самые удаленные от ж не принадлежат кругу 5-, т. е. ||ж — У (0) || > I, что при £ > 0 влечет ||ж — У (£) || > I. Поэтому У (£) при £ > 0 не может быть концом хорды Q, а следовательно, (> 0^1 = С(£), т. е. Сд, равный С(0), не может находиться между ж = С(—||Сд — ж||) и Ql. По лемме ж € 5х,д, следовательно, ж не может быть между Сд и Ql (тогда бы выполнялось условие ||Сд — ж|| < ||Сд — Ql|| < 5/2, что влечет ж € 5х,д). Остается одна возможность — Ql находится между Сд и ж.
2. Поскольку центры кругов и БХуд находятся по разные стороны от их общей хорды Q, их пересечение состоит из объединения двух сегментов, отсекаемых хордой Q от этих кругов. Причем ни один из таких сегментов не содержит центра круга, от которого он отсечен, поэтому каждый из сегментов находится в круге, для которого Q — диаметр. Следовательно, и объединение этих сегментов также лежит в этом круге. □
Пусть г := А(у) — у, г := ||Г||. Используя шаговую оценку погрешности I и оценку (8), запишем
г < ||у — а|| + ||А(у) — а|| < I + !с. (10)
Найдем параметры круга Аполлония — величину его диаметра 5 и центр О. Согласно теореме 1 поиск сводится к одномерному на прямой I. Введем на ней координаты так, что элементу у будет соответствовать О, а элементу А(у) — г. Тогда для координат £1>2 пересечения круга Аполлония с прямой I согласно (8) будем иметь уравнение |г — £| = с|£|. Отсюда
£1 = 1 I , £2 = ----.
1 + с 1 — с
Следовательно,
О = у + ^ = у +-Ч. 5 = £2 — £, = т2^. (11)
2 г 1 — с2 1 — с2
На перпендикуляре к прямой, проходящей через у и А(у), восстановленном в центре О, точки, более близкие к центру, находятся ближе также и к итерации у. Значит, чтобы проверить, будет ли диаметр 5 принадлежать шару 5-, достаточно проверить, будут ли его концы в шаре 5-, что эквивалентно соотношению
^(О — у)2 + 52/4 < й. (12)
Найдем г, для которого в (12) выполняется равенство
п2 г2 с2г2
(О - у)2 + Т - (Г—^ + (Т—^¥ =
Отсюда искомое значение г есть
А - ^лтж ■ ■“>
В силу формул (11) выражения для (О — ж)2 и 52 растут с увеличением г, поэтому 5 С
5- г < г.
Случай г < Г. Точная релаксация и оценка ее погрешности определяются по формулам
у* = О = у + —Ц; (14)
1 с2
сг
1с
2
(15)
Случай г > Р. Вначале найдем центр к = у + кг/г общей хорды где к — расстояние от этого центра до у — удовлетворяет системе
к + г = - у||,
Т2/4 - г2 = ?2 - к2 = ?2
(г — расстояние между центрами круга Аполлония и общей хорды). Преобразуем эту систему с помощью формул (11):
г
к + г
1 — с2 ’
г ТТ2 с2 г2
(к — г)---------------- = ?2--------------— = ?2 —
1 - с2 4 (1 - с2)2
Отсюда
у2 = к = у + Ц 1 + і?2*!-і2» І. (17)
Половина длины хорды Q есть
I* = 2 ||Q| = Р. (18)
Итак, расчетные формулы в первом случае будут (14) и (15), во втором — (16)—(18).
2.2. Сводные формулы
Объединяя (14) с (17) и (15) с (16), (18), получим формулы для точной релаксации М: г = А(у) — у, г = ||г||, г = !(1 — с2)/л/ 1 + с2; (19)
'у + г(1+ I2(1 — с2)/г2)/2 при г < г,
у* = ^ г (20)
у + 1 _ с2 при г > г;
-у/?2 - (г + ?2(1 - с2)/г)2/4 при г < г,
^ сг ^ (21)
= --------о при г > г.
1 с2
2.3. Оценки погрешности и вычислительной трудоемкости ТР
Определим оценку погрешности I* от точной релаксации через текущую оценку I.
В случае г < Г наибольшее значение I* достигается при г = Г и равно !с/\/1 + с2.
Следовательно, I* < !с/^ 1 + с2.
В случае г > Г наибольшая оценка I* соответствует наименьшему значению Л. С помощью формулы (16) видно, что Л(г) строго выпукла и стремится к при г \ 0 и
г ^ +то. Поэтому существует единственный абсолютный минимум, который соответствует корню производной
7/ 1( !2(1 — с2)'
0 = Л = ^1 г2"
находимому без труда: Г := !\/1 — с2. Очевидно, что Г > Г и
I + с! = IV1 + с V1 + с > (1\/1 + с V1 — с = г.
Следовательно, г принадлежит области задания Л, и
Лт1п = Г = IV 1 — с2, I* < ^сР — Лт;п = сй.
Таким образом, когда ||А(х) — х|| = ^\/1 — с2 , точная релаксация не заменяет значение А(х) иным и, естественно, дает оценку погрешности следующей итерации — ту же (Ы), что и базовый алгоритм. Во всех прочих случаях она обеспечивает I* < ^ при в < с. Величина в на каждой итерации своя. При желании ее можно извлечь как функцию некоторого параметра 5 € (0,1) из формулы (15) или из формул (16), (18), полагая г = ^
с5 1 — с2
при 5 <
1 - с2 " - v/тT^’
I 1 _ с2
V1 - (5 + (1 - с2)/^)2 /4 при 5 > ^Гс2 •
Как и следовало ожидать, при 5 = v/T— с2 параметр 5 равен 1.
Вычислив к(Р) в силу (16), а затем в силу (18), также получим ?с/\/1 + с2. Следовательно, формулы (15), (16), (18) определяют непрерывную функцию ?2(г), г Є [0,? + ?с].
С вычислительной точки зрения экономнее использовать последовательность } вместо {?&}. Тогда вычислительные затраты на один шаг точной релаксации есть:
— вычисление разности г = А(у) - у (применение алгоритма А не входит в состав собственно ТР);
22
— вычисление скалярного квадрата г2 = г 2;
— умножение на (1 - с2)2(1 + с2)-1, т. е. вычисление Р2;
— умножение элемента пространства на число и сложение элементов в формулах (14) либо (17);
— одно умножение в (15) либо одно вычитание в (18) для вычисления ?|+1.
В случае г > Р для вычисления значений к и к2 потребуются еще пять арифметических операций и извлечение квадратного корня из вещественного числа.
Итого: 3... 8 арифметических операций и (что существеннее) вычисление скалярного произведения + два сложения элементов + умножение элемента на число.
2.4. Точные релаксации простых итераций в скалярном пространстве
Наиболее наглядно проявление метода точной релаксации в скалярном пространстве. Шар 5^ превращается в сегмент [х—х+^]. Для описания зоны Аполлония введем обозначения
г := А(х) — х, г := |г |. (22)
Получим расчетные формулы для линейной сходимости (р = 1). Таковой сходимостью обладает, например, модифицированный метод Ньютона решения нелинейного уравнения $(х) = 0 по схеме ж* = х — 7-1(х0)$(х), где $ : Д С ^ Дга; 7 — матрица Якоби для функции $.
Традиционно считается, что оценка (4) при р =1 имеет ценность, когда с < 1. Будем и здесь пока полагать с < 1. Тогда в линейном случае зона Аполлония есть сегмент
Я
А : =
х + ----, X +
1 + с 1 — с
(23)
(Здесь допущена небольшая вольность речи: когда г < 0, левый конец последнего сегмента становится больше правого.)
Пересечение двух сегментов 5^ и есть сегмент. Его центр, согласно ТР, назначается
следующей итеративной точкой, а половина его длины — оценкой ^+1 погрешности точки
Х^+1
Вывод расчетных формул для х* и ^* несложен:
х + (і sgn г +-----------) /2 при і <
V 1 + сг
г
х* = { _1 + с) 1 - с (24)
і г г
х + -----------2 при і > -------;
1 — с2 1 — с
(і - /2 при й <
й. = Г ^ 1 - с (25)
----2 при і > -------.
1 — с2 1 — с
Если неизвестна оценка погрешности начального приближения то она полагается равной +то, что приводит на первом шаге к использованию вторых строк в формулах (24), (25): х1 = х + г/(1 — с2), іі = гс/(1 — с2).
Дополнительные вычислительные затраты на итерацию точной релаксации определяются из формул (24), (25). Это — одна операция сравнения и пять либо восемь (вторые либо первые строки формул) арифметических операций.
Легко установить, что всегда
І. < Т+-- (26)
1 + с
г
Равенство в оценке (26) возможно лишь в редких случаях, когда || А(у)— у|| = г = ^—^с. Поэтому практическая ценность модификации выше, чем коэффициент с/(1 + с) в (26) вместо с в (7).
При возрастании с к 1 скорость сходимости метода точной релаксации уменьшается. Тем не менее всегда гарантируется уменьшение оценки погрешности более чем в два раза на шаге.
Примечательно, что если оценка начальной погрешности известна, то метод ТР сходится к решению и в случае с =1. Тогда расчетные формулы ТР упрощаются до
х* = х + ^ sgn г + г/2^/2 при ^* = ^ — г/2^/2.
При этом на каждой итерации будет выполняться условие г < 2^ и погрешность продолжает уменьшаться более чем в два раза.
3. Численные эксперименты
Подвергнем точной релаксации модифицированный метод Ньютона (ММН) решения скалярного уравнения $(х) = 0:
хк+1 = хк — 7-1(х0)$(хк) =: А(хк), (27)
где 7 := $. Это метод простой итерации с линейной сходимостью. Поскольку речь идет о практической реализации, в вычислительном процессе должен присутствовать критерий остановки, например, по малости невязки ||$(хк)|| < е или по малости оценки погрешности ^ < е, ||хк — а|| < ^, где а — искомое решение уравнения, которое, очевидно, является неподвижной точкой базового алгоритма А. Для точной релаксации наиболее удобен последний критерий. Его и возьмем для исследования. В этом случае о сходимости и скорости сходимости ММН имеется теорема Мысовских (см. теорему 4 в [7]), однако ее заключение о скорости сходимости имеет форму
||хк — а|| < дк-1|х1 — а||, д := Рм + РМ/4 £ [0,1), (28)
не пригодную для использования в ТР. Но непосредственно из доказательства этой теоремы можно извлечь достаточно информации нужного вида. Синтез посылок из теоремы 4 в работе [7] с упомянутой информацией без сужения с банахова пространства на скалярное таков:
Теорема 3. Если уравнение $(х) = 0 имеет решение а в шаре 5* и выполнены условия:
1) существует оператор 7-1(х0) и |7-1(х0)| < г0, 7 :=
2) ||/(х)|| < Ь Ух € 55+Рм"о/2;
3) Рм := Т0Ь4 < 2\[2 — 2,
то решение а единственно в шаре 5* и к нему сходятся последовательные приближения (27), а быстрота сходимости характеризуется рекуррентными формулами с0 := Рм/2, и для к = 1, 2,...
IIх а11 < ск—1^й-1 : ) ск : Рм + 2 ^к— 1. (29)
Отметим, что в этой теореме условия 2 и 3 взаимозависимы. Это представляет заметное неудобство в ее непосредственном применении. Чтобы уйти от дополнительной проблемы,
огрубим результаты, заменив в радиусе шара из условия 2 величину Рм на какую-то ее оценку сверху, например на 2\/2 — 2 < 0.83 или даже на 1: 2/) ||^//(х)| < ЬУж Є $1оМо.
Первое приближение по ММН совпадает с первым приближением по обычному методу Ньютона и имеет квадратичное изменение оценки погрешности:
ж1 — а = х0 — а — 7-1(ж0)д(ж0) =
= ж0 — а — 7-1(ж0) (#(а) + 7(ж0)(ж0 — а) + о(ж0 — а)) = 7-1(ж0)о(ж0 — а),
где ||о(ж0 — а)|| < Ь||ж0 — а||2/2. Отсюда
Цх1 — а|| < г0Ь||ж0 — а||2/2. (30)
Квадратичная (30) и линейная (28) оценки приводят к существенно различным формулам
точной релаксации, что повышает сложность программирования ТР. Поэтому огрубим (30) до линейной
Р
1 Рм 0 0
||ж — а|| < — ||ж — а|| = С0|ж — а||, (31)
что позволит на всех итерациях использовать формулы (24), (25) с переменным с.
Были проведены численные эксперименты с несколькими скалярными уравнениями по единой схеме.
Ш!аг 1. Выбиралась скалярная функция f с известным корнем а. На некотором расстоянии от него выбиралась начальная точка ж0, т. е. начальная погрешность (не ее оценка) была ^0 = |ж0 — а|.
Ш!аг 2. Определялся сегмент
а := [ж0 — 1.5^0, ж0 + 1.5^0] 3 [ж0 — л/2^0, ж0 + л/2^0].
Шаг 3. Находилась оценка Ь > эир^,ео. f //(ж).
Ш!аг 4. Выяснялось, удовлетворяют ли найденные параметры условию теоремы 3 Рм = Ьг0^0 < 2\/2 — 2 или немного более сильному условию Рм < 0.83. Если нет, то возврат на шаг 1.
Ш!аг 5. Программой на С производилось десять итераций по ММН, и оценка погрешности велась с переменным с по формулам из (29). Из той же начальной точки у0 := ж0 и с той же начальной оценкой погрешности ^0° := ^0 производилось десять итераций по методу ТР с базовым алгоритмом А из (27) согласно формулам (24) и (25), к которым были добавлены формулы из (29) для вычисления {ск0)}.
Шаг 6. Из той же начальной точки у02 := у0^ := ж0 производилось по десять итераций по методу ТР в двух вариантах ухудшения оценки характерного параметра Рм совместно с ухудшением оценки начальной погрешности:
РМ ) := рм + (2^2 — 2 — рм)/2,
^1 := (РІ1}/Рм — і) /2 + 1, ^ := ^1, Ь^ := Р^^. (32)
Переменный коэффициент с, используемый в методе ТР, определялся по формулам
с01} := Р^/2, сЦ := Р^ + ЬГ1)^і-)1. (33)
Второе огрубление сильнее. Вначале := Рм + (2\^2 —2 — Рм)3/4. Затем применяются
формулы (32), (33) с заменой индекса 1 на 2.
Пояснение к шагу 6. Ухудшение параметра Рм назначается таким, чтобы сохранилось условие 3 теоремы 3, а оценки начальных погрешностей назначаются так, чтобы они не обеспечивали ухудшение параметра Рм.
В табл. 1-4 хк — итерации по модифицированному методу Ньютона; у^ — итерации по точной релаксации с коэффициентами {с£0)}; У1 — итерации по ТР с коэффициентами {с^}; у£2) — итерации по методу точной релаксации с коэффициентами {ск2)}, им соответствуют оценки погрешностей ^, ^к,0), ^к_2). Отметим, что хк не зависят от ко-
эффициента с.
X
Пример 1. Уравнение д(х) := —------------= 0.
х2 + 6х + 5
Шаг 1. х0 = 0.15 , = 0.15.
Шаг 2. а = [—0.075, 0.375].
Шаг 3. Для удобства оценки константы Липшица для д' разложим д в сумму простейших:
5 1 ' —5 1
д(х) = 77—— 77—“ГТ, д (х) = 77—ГТТ7 +
4(х + 5) 4(х + 1) ’ 4(х + 5)2 4(х + 1)2
"( \ 5 1 д (х) =
2(х + 5)3 2(х + 1)3’
Вторая производная имеет единственный корень х'' в (—1, +то). Оценим его снизу:
$5 1 5 — ^5 5 — 1 ,
' ‘ х = ^--------> ------= 4.
х + 5 х + 1 1^5 — 1 2 — 1
Выясним, будет ли вторая производная монотонно расти на сегменте а. Ее производная
Ш, \ —15 3
д (х) = ^ +
2(х + 5)4 2(х + 1)4
также имеет единственный корень х''' в (—1, +то). Так как д''(1 + 0) = — то, очевидно, что д'''(х) > 0 Ух € (—1, 4). Следовательно, если а € (—1, 4), то тахд'' достигается на левом конце сегмента а, т. е. д''(—0.075) < 0.611 =: Ь.
Шаг 4. Рм = 0.6456.
Шаги 5,6. Расчеты показали, что помимо значительно большей скорости уменьшения оценки погрешности (сравнительно с ММН) и сами итерации по методу ТР существенно быстрее сходились к решению как при хороших оценках на Ь и ^0, так и при сильном огрублении (табл. 1).
Отметим, что выпуклость д в сочетании с д(х0) < 0 порождает “перескок” итераций по ММН через решение, нарушая их монотонность.
Пример 2. Уравнение д(х) := ех/3 — 1 = 0.
Шаг 1. х0 = —1, ^0 = 1.
Шаг 2. а = [—2.5, 0.5].
Шаг 3. Ь = е°'5/3/9.
Шаг 4. Рм = Ье1/33 = е/3 < 0.550.
Шаги 5, 6. Расчеты показывали сходное с предыдущим примером поведение итераций по методу точной релаксации и модифицированному методу Ньютона. Как и в первом примере, происходил “перескок” итераций по ММН через решение (табл. 2).
Таблица 1.
к к Ук / (хк) / (Ук) 4 ^(0) “к уГ / (У(1)) 4:) (2) Ук / (ук2)) “(2) “к
0 1.500е-01 1.500е—01 2.533е—02 2.533е—02 1.500е—01 1.500е—01 1.500е—01 2.533е—02 1.606е—01 1.500е—01 2.533е—02 1.659е—01
1 —2.848е-02 7.539е—03 —5.896е—03 1.494е—03 4.842е—02 7.539е—03 4.484е—03 8.920е—04 1.510е—02 2.901е—03 5.782е—04 1.882е—02
2 1.307е—02 6.015е—04 2.574е—03 1.202е—04 3.631е—02 6.015е—04 —4.805е—03 —9.667е—04 5.810е—03 —7.599е—03 —1.534е—03 8.323е—03
3 —5.066е—03 4.357е—05 —1.019е—03 8.713е—06 2.628е—02 4.357е—05 3.050е—05 6.100е—06 9.736е—04 —4.715е—04 —9.435е—05 1.195е—03
4 2.118е—03 3.130е—06 4.225е—04 6.259е—07 1.846е—02 3.130е—06 —6.499е—05 — 1.300е—05 7.081е—05 3.121е—04 6.240е—05 4.118е—04
5 —8.594е—04 2.247е—07 —1.720е—04 4.493е—08 1.265е—02 2.247е—07 —3.224е—06 —6.447е—07 9.041е—06 —1.699е—05 —3.398е—06 8.272е—05
6 3.531е—04 1.613е—08 7.058е—05 3.226е—09 8.512е—03 1.613е—08 2.605е—06 5.210е—07 3.213е—06 3.108е—05 6.217е—06 3.464е—05
7 —1.443е—04 1.158е—09 —2.887е—05 2.316е—10 5.652е—03 1.158е—09 —5.836е—08 — 1.167е—08 5.497е—07 1.476е—06 2.952е—07 5.037е—06
8 5.912е—05 8.312е—11 1.182е—05 1.662е—11 3.718е—03 8.312е—11 1.217е—07 2.434е—08 1.327е—07 —1.626е—06 —3.252е—07 1.935е—06
9 —2.420е—05 5.967е—12 —4.839е—06 1.193е—12 2.430е—03 5.967е—12 5.975е—09 1.195е—09 1.698е—08 — 1.567е—08 —3.133е—09 3.248е—07
10 9.907е—06 4.284е—13 1.981е—06 8.567е—14 1.582е—03 4.284е—13 —4.941е—09 —9.881е—10 6.068е—09 4.134е—08 8.267е—09 4.462е—08
Таблица 2.
к к Ук / (хк) / (Ук) <4 “(0) “к у^ (/ у 4^ (2) Ук (/ 2) “(2) “к
0 —1.000е+00 — 1.000е+00 —2.835е—01 —2.835е—01 1.000е+00 1.000е+00 — 1.000е+00 —2.835е—01 1.127е+00 —1.000е+00 —2.835е—01 1.190е+00
1 1.868е—01 —3.450е—02 6.42е—02 — 1.143е—02 2.748е—01 3.450е—02 4.792е—03 1.599е—03 1.221е—01 2.535е—02 8.486е—03 1.649е—01
2 —8.221е—02 —1.897е—03 —2.703е—02 —6.320е—04 1.718е—01 1.897е—03 — 1.126е—02 —3.748е—03 1.226е—02 —6.664е—02 —2.197е—02 7.293е—02
3 3.096е—02 —9.478е—05 1.037е—02 —3.159е—05 1.025е—01 9.478е—05 —5.091е—04 —1.697е—04 1.506е—03 —4.700е—03 — 1.566е—03 1.099е—02
4 —1.247е—02 —4.710е—06 —4.149е—03 — 1.570е—06 5.922е—02 4.710е—06 4.541е—04 1.514е—04 5.428е—04 2.650е—03 8.837е—04 3.638е—03
5 4.899е—03 —2.340е—07 1.634е—03 —7.799е—08 3.351е—02 2.340е—07 —4.884е—06 —1.628е—06 8.382е—05 —2.196е—04 —7.319е—05 7.685е—04
6 —1.944е—03 —1.162е—08 —6.476е—04 —3.875е—09 1.872е—02 1.162е—08 8.094е—06 2.698е—06 8.942е—06 2.518е—04 8.393е—05 2.972е—04
7 7.680е—04 —5.774е—10 2.560е—04 —1.925е—10 1.039е—02 5.774е—10 2.787е—07 9.291е—08 1.127е—06 3.306е—06 1.102е—06 4.869е—05
8 —3.040е—04 —2.869е—11 —1.013е—04 —9.562е—12 5.738е—03 2.869е—11 —4.000е—07 —1.333е—07 4.485е—07 —7.568е—06 —2.523е—06 8.251е—06
9 1.202е—04 —1.425е—12 4.008е—05 —4.750е—13 3.163е—03 1.425е—12 —1.053е—08 —3.512е—09 5.895е—08 —4.400е—07 — 1.467е—07 1.123е—06
10 —4.757е—05 —7.084е—14 —1.586е—05 —2.354е—14 1.741е—03 7.084е—14 1.746о—08 5.818е—09 1.929е—08 2.959е—07 9.8630—08 3.867е—07
Метод точных релаксаций 83
Таблица 3.
к к Ук / (хк) / (Ук) 4 ^(0) “к У^ (/ У 4:) (2) Ук (2) (У (/ “к
0 1.000е+00 1.000е+00 3.956е—01 3.956е—01 1.000е+00 1.000е+00 1.000е+00 3.956е—01 1.127е+00 1.000е+00 3.956е—01 1.190е+00
1 1.496е-01 1.665е—01 5.113е—02 5.705е—02 2.748е—01 1.665е—01 1.203е—01 4.092е—02 2.472е—01 9.660е—02 3.272е—02 2.869е—01
2 3.969е-02 4.479е—02 1.332е—02 1.504е—02 1.718е—01 4.479е—02 —2.717е—02 —9.015е—03 9.968е—02 —6.495е—02 —2.142е—02 1.253е—01
3 1.106е-02 1.204е—02 3.694е—03 4.023е—03 1.025е—01 1.204е—02 1.712е—02 5.722е—03 3.321е—02 1.023е—02 3.416е—03 5.014е—02
4 3.121е-03 3.238е—03 1.041е—03 1.080е—03 5.922е—02 3.238е—03 —3.086е—03 — 1.028е—03 1.301е—02 —9.359е—03 —3.115е—03 1.549е—02
5 8.835е-04 8.703е—04 2.945е—04 2.901е—04 3.351е—02 8.703е—04 1.212е—03 4.039е—04 2.995е—03 2.781е—04 9.272е—05 5.851е—03
6 2.504е-04 2.339е—04 8.346е—05 7.798е—05 1.872е—02 2.339е—04 —4.494е—04 — 1.498е—04 1.147е—03 — 1.979е—04 —6.596е—05 3.630е—04
7 7.096е-05 6.288е—05 2.365е—05 2.096е—05 1.039е—02 6.288е—05 1.647е—04 5.491е—05 4.236е—04 2.389е—05 7.963е—06 1.412е—04
8 2.011е-05 1.690е—05 6.705е—06 5.634е—06 5.738е—03 1.690е—05 —6.013е—05 —2.004е—05 1.550е—04 —1.6466—05 —5.487е—06 3.062е—05
9 5.702е-06 4.543е—06 1.901е—06 1.514е—06 3.163е—03 4.543е—06 2.191е—05 7.305е—06 5.654е—05 2.201е—06 7.336е—07 1.196е—05
10 1.616е—06 1.221е—06 5.388е—07 4.071е—07 1.741е—03 1.221е—06 —7.982е—06 —2.661е—06 2.060е—05 — 1.515е—06 —5.051е—07 2.820е—06
Таблица 4.
к к Ук / (хк) / (Ук) 4 “к У^ / (У(1)) 4^ (2) Ук / (у!2)) “к
0 1.047е+00 1.047е+00 1.913е+00 1.913е+00 1.047е+00 1.047е+00 1.047е+00 1.913е+00 1.096е+00 1.047е+00 1.913е+00 1.121е+00
1 —2.283е—01 5.087е—02 —4.546е—01 1.017е—01 3.655е—01 5.087е—02 3.758е—02 7.514е—02 8.645е—02 3.074е—02 6.147е—02 1.040е—01
2 7.478е—02 5.665е—03 1.495е—01 1.133е—02 2.997е—01 5.665е—03 —1.939е—02 —3.877е—02 2.949е—02 —3.225е—02 —6.449е—02 4.107е—02
3 —2.488е—02 6.109е—04 —4.976е—02 1.222е—03 2.392е—01 6.109е—04 2.604е—03 5.207е—03 7.500е—03 6.506е—05 1.301е—04 8.755е—03
4 8.291е—03 6.564е—05 1.658е—02 1.313е—04 1.860е—01 6.564е—05 —2.128е—03 —4.255е—03 2.768е—03 —1.781е—04 —3.563е—04 1.951е—04
5 —2.764е—03 7.051е—06 —5.527е—03 1.410е—05 1.414е—01 7.051е—06 6.006е—05 1.201е—04 5.806е—04 —1.448е—05 —2.897е—05 3.141е—05
6 9.212е—04 7.573е—07 1.842е—03 1.515е—06 1.054е—01 7.573е—07 —1.321е—04 —2.641е—04 1.467е—04 6.597е—06 1.319е—05 1.033е—05
7 —3.071е—04 8.133е—08 —6.142е—04 1.627е—07 7.726е—02 8.133е—08 —8.776е—06 —1.755е—05 2.343е—05 —1.016е—06 —2.031е—06 2.715е—06
8 1.024е—04 8.736е—09 2.047е—04 1.747е—08 5.593е—02 8.736е—09 6.258е—06 1.252е—05 8.398е—06 7.187е—07 1.437е—06 9.803е—07
9 —3.412е—05 9.382е—10 —6.824е—05 1.876е—09 4.009е—02 9.382е—10 —3.071е—07 —6.142е—07 1.833е—06 —3.825е—08 —7.650е—08 2.233е—07
10 1.137е—05 1.008е—10 2.275е—05 2.015е—10 2.852е—02 1.008е—10 6.738е—07 1.348е—06 7.487е—07 8.762е—08 1.752е—07 9.747е—08
84 С. Е. Михеев
Пример 3. Уравнение g(x) := ex/3 — 1 = 0.
Шаг 1. x0 = 1, d0 = 1.
Шаг 2. а = [0.5, 2.5].
Шаг 3. L = e2'5/3/9.
Шаг 4. Pm = Le-1/33 = е/3 < 0.55.
Шаги 5, 6. Расчеты показали, что при хороших оценках величин L и d0 метод точной релаксации помимо значительно большей скорости уменьшения оценки погрешности (сравнительно с ММН) обеспечивает и более быструю сходимость. При грубых оценках этих величин метод ТР также обеспечивает значительную скорость уменьшения оценки погрешности по сравнению с ММН, но сами итерации ММН сходятся быстрее, чем в методе точной релаксации. Что, правда, не гарантируется теорией (табл. 3).
Отметим, что выпуклость g в сочетании с g(x0) > 0 благоприятна для ММН. Она обеспечивает монотонность итераций.
Пример 4. Уравнение g(x) := x + sin x = 0.
Шаг 1. x0 = n/3, d0 = п/3.
Шаг 2. а = [—п/6,п5/6].
Шаг 3. L = 1.
Шаг 4. Pm =------L----0d0 = 2апП5/12 п < 0.698 < 2^2 — 2.
1 + cos x0 v 3 6
Шаги 5, 6. Картина сходна с той, что была в примерах 1 и 2 (табл. 4).
Ни выпуклости, ни вогнутости в корне уравнения нет, но “перескок” через корень происходил, и расчеты подтвердили эффективность применения метода ТР для этого случая.
Заключение
Метод точных релаксаций основан на очень простом не вызывающем сомнений принципе минимальности. Для некоторых типов базовых алгоритмов и пространств (как в разд. 2, например) можно получить простые расчетные формулы, вычислительная трудоемкость которых легко оценивается. Оценка погрешности ТР по худшему варианту может совпадать с оценкой погрешности итерации по базовому алгоритму, и тогда ТР совпадает с этой итерацией. Поскольку худший вариант встречается далеко не всегда, реальное улучшение сходимости по методу ТР, как показывают численные эксперименты, оказывается значительным.
Описание вывода расчетных формул для метода точной релаксации существенно длиннее, чем изложение собственно принципа минимальности, поэтому они здесь были получены для весьма обширного, но лишь одного типа итеративных методов. Без должного внимания остались доступные для метода ТР не менее интересные типы, как-то: пространства с иными нормами (например, ||у||д = ^/У1 |УгК), с иной оценкой скорости сходимости (р > 1), с другим строением итеративной информации, г-точечные методы.
Численные эксперименты в одномерном вещественном пространстве показали, что текущая оценка погрешности, полученная с использованием метода точной релаксации, и в худшем случае существенно меньше оценки базовых простых итераций. Однако сама погрешность (не ее оценка) в некоторых случаях может в простой итерации быть меньше, чем в ее ТР.
Список литературы
[1] Михеев С.Е. Нелинейные методы в оптимизации. СПб.: Изд-во Санкт-Петербург. гос. ун-та, 2001. 27б с.
[2] Фаддеев Д.К., Фаддеева В.Н. Вычислительные методы линейной алгебры. М.: Физматгиз, 19б0. бБб c.
[3] Ортега Дж., РЕйнволдт В. Итерационные методы решения нелинейных систем уравнений со многими неизвестными: Пер. с англ. М.: Мир, 197Б. ББВ с.
[4] Хейгеман Л., Янг Д. Прикладные итерационные методы. М.: Мир, 19Вб. 44б c.
[Б] Miheev S.E. Contraction of attainability domains in a game of pursuit jj Nova J. Math.: Game
Theory and Algebra. N.Y.: Nova Sci. Publ., 1997. Vol. б, N 2j3. P. 147-1б1.
[6] Михеев С.Е. Релаксационное ускорение на основе областей достижимости jj Вест. Санкт-
Петербург. гос. ун-та. Сер. 1. 1999. Вып. 3, № 1Б. C. 29-35.
[7] Мысовских И.П. О сходимости метода Л.В. Канторовича решения функциональных уравнений и его применениях j j Докл. АН СССР. 1950. Т. 70, № 4. С. БбБ-БбВ.
Поступила в редакцию їв августа 2QQ4 г., в переработанном виде — 13 сентября 2QQ6 г.
