1611688890-f641c9ec8276824e4686da772eb56520 (826652), страница 60
Текст из файла (страница 60)
По этой причинео степени близости x(k) к x⋆ судят на основании косвенных признаков, важнейшим среди которых является величина невязки решения.Невязка определяется как разность левой и правой частей уравненияпосле подстановки в него приближения к решению, и в нашем случаеэто Ax(k) − b. При этом конкретное применение принципа релаксацииможет заключаться в том, что на каждом шаге итерационного процесса стремятся уменьшить абсолютные значения компонент вектораневязки либо её норму, либо какую-то зависящую от них величину.В этом смысле методы Якоби и Гаусса-Зейделя можно рассматриватькак итерационные процессы, в которых также осуществляется релаксация, поскольку на каждом их шаге компоненты очередного приближения вычисляются из условия зануления соответствующих компонентневязки на основе уже полученной информации о решении.
Правда,это делается «локально», для отдельно взятой компоненты, и без учётавлияния результатов вычисления этой компоненты на другие компоненты невязки.Различают релаксацию полную и неполную, в зависимости от того,добиваемся ли мы на каждом отдельном шаге итерационного процесса(или его подшаге) наибольшего возможного улучшения рассматриваемой функции от погрешности или нет.
Локально полная релаксацияможет казаться наиболее выгодной, но глобально, с точки зрения сходимости процесса в целом, тщательно подобранная неполная релаксациянередко приводит к более эффективным методам.Очень популярной реализацией высказанных выше общих идей является метод решения систем линейных алгебраических уравнений, вкотором для улучшения сходимости берётся «взвешенное среднее» значений компонент предшествующей x(k) и последующей x(k+1) итерацийметода Гаусса-Зейделя.
Более точно, зададимся вещественным числомω, которое будем называть параметром релаксации, и i-ую компоненту26 Отлатинского слова «relaxatio» — уменьшение напряжения, ослабление.3663. Численные методы линейной алгебрыочередного (k + 1)-го приближения положим равной(k+1)ωxi(k)+ (1 − ω)xi ,(k)где xi — i-ая компонента приближения, полученного в результате k(k+1)— i-ая компонента приближения, котороего шага алгоритма, а xi(k)(k)(k+1)(k+1)было бы получено на основе x(k) и x1, .
. . , xi−1 , xi+1 , . . . , xnс помощью метода Гаусса-Зейделя. Псевдокод получающегося итерационного алгоритма, который обычно и называют методом релаксациидля решения систем линейных алгебраических уравнений, представленв Табл. 3.7.Таблица 3.7. Псевдокод метода релаксациидля решения систем линейных уравненийk ← 0;выбираем начальное приближение x(0) ;DO WHILE ( метод не сошёлся )DO FOR i = 1 TO n(k+1)xi(k)← (1 − ω) xi+END DOnXi−1Xω (k)(k+1)aij xj −aij xjbi −aiij=i+1j=1k ← k +1;END DOРасчётные формулы этого метода перепишем в виде(k+1)aii xi+ωi−1Xj=1(k+1)aij xj(k)= (1 − ω) aii xi−ωnX(k)aij xj+ ωbi ,j=i+1i = 1, 2, . .
. , n,3.9. Стационарные итерационные методы367k = 0, 1, 2, . . . . Далее, используя введённые выше в §3.9е матрицы L̃, Dи Ũ , можно придать этим соотношениям более компактный видD + ω L̃ x(k+1) = (1 − ω)D − ω Ũ x(k) + ωb,откуда−1−1x(k+1) = D + ω L̃(1 − ω)D − ω Ũ x(k) + D + ω L̃ωb,k = 0, 1, 2, .
. . .В зависимости от конкретного значения параметра релаксации принято различать три случая:если ω < 1, то говорят о «нижней релаксации»,если ω = 1, то имеем итерации Гаусса-Зейделя,если ω > 1, то говорят о «верхней релаксации».27Последний случай может показаться экзотичным, но во многих ситуациях он действительно обеспечивает улучшение сходимости итерацийв сравнении с методом Гаусса-Зейделя.
Несколько упрощённое объяснение этого явления может состоять в том, что если направление отx(k) к x(k+1) оказывается удачным в том смысле, что приближает к искомому решению, то имеет смысл пройти по нему и дальше, за x(k+1) .Это соответствует случаю ω > 1.Важно отметить, что метод релаксации также укладывается в изложенную в §3.9в схему итерационных процессов, порождаемых расщеплением матрицы решаемой системы уравнений.
Именно, мы берёмA = Gω − Hω с матрицамиGω = D + ω L̃,Hω = (1 − ω)D − ω Ũ .Необходимое и достаточное условие сходимости метода релаксации принимает поэтому вид−1ρ GωHω < 1.Для некоторых специфичных, но очень важных задач математической физики значениерелаксационного параметра ω, при котором величина ρ G−1ω Hω достигает минимума, находится относительно просто. В более сложных задачах для оптимизации ω требуется весьма27 В англоязычной литературе по вычислительной линейной алгебре этот метод обычно обозначают аббревиатурой SOR(ω), которая происходит от термина«Successive OverRelaxation» — последовательная верхняя релаксация.3683. Численные методы линейной алгебры−1трудный анализ спектра матрицы перехода GωHω из представления(3.103). Обзоры состояния дел в этой области читатель может найти в[45, 49, 78, 97, 98].−1Предложение 3.9.3 Если Cω = D + ω L̃(1 − ω)D − ω Ũ — матрица оператора перехода метода релаксации, то ρ(Cω ) ≥ |ω − 1|.
Какследствие, неравенство 0 < ω < 2 на параметр релаксации необходимодля сходимости метода.Доказательство. Прежде всего, преобразуем матрицу Cω для придания ей более удобного для дальнейших выкладок вида:Cω ====−1D + ω L̃(1 − ω)D − ω Ũ−1D(I + ωD−1 L̃)(1 − ω)D − ω Ũ−1 −1(1 − ω)D − ω ŨI + ωD−1 L̃D−1I + ωD−1 L̃(1 − ω)I − ωD−1 Ũ .Желая исследовать расположение собственных чисел λi (Cω ) матрицы Cω , рассмотрим её характеристический полином−1φ(λ) = det(Cω − λI) = det I + ωD−1 L̃(1 − ω)I − ωD−1 Ũ − λI= pn λn + pn−1 λn−1 + . .
. + p1 λ + p0 ,в котором pn = (−1)n по построению. Свободный член p0 характеристического полинома может быть найден как φ(0):−1p0 = det Cω = det I + ωD−1 L̃(1 − ω)I − ωD−1 Ũ= det (I + ωD−1 L̃)−1 · det (1 − ω)I − ωD−1 Ũ= det (1 − ω)I − ωD−1 Ũ = (1 − ω)n ,коль скоро матрица (I + ωD−1 L̃) — нижняя треугольная и диагональными элементами имеет единицы, а (1 − ω)I − ωD−1 Ũ — верхняятреугольная, с элементами (1 − ω) по главной диагонали.С другой стороны, по теореме Виета свободный член характеристического полинома матрицы, делённый на старший коэффициент, равен3.10.
Нестационарные итерационные методы369произведению его корней, т. е. собственных чисел матрицы, умноженному на (−1)n (см., к примеру, [22]), и поэтомуnYi=1λi (Cω ) = (1 − ω)n .Отсюда необходимо следуетmax |λi (Cω )| ≥ |ω − 1|,1≤i≤nчто и доказывает Предложение.Теорема 3.9.5 (теорема Островского-Райха) Если в системе линейных алгебраических уравнений Ax = b матрица A является симметричной положительно определённой, то для всех значений параметраω ∈ ]0, 2[ метод релаксации сходится к решению из любого начальногоприближения.Доказательство опускается. Читатель может найти его, к примеру,в книгах [13, 97, 98]. Обоснование теоремы Островского-Райха будеттакже дано ниже в §3.12 как следствие теоремы Самарского, дающейдостаточные условия сходимости для итерационных методов весьма общего вида.3.10Нестационарные итерационныеметоды для линейных систем3.10аТеоретическое введениеВ этом параграфе для решения систем линейных алгебраическихуравнений мы рассмотрим нестационарные итерационные методы, которые распространены не меньше стационарных.
В основу нестационарных методов могут быть положены различные идеи.В качестве первого примера рассмотрим метод простой итерации(3.104)x(k+1) ← (I − τ A) x(k) + τ b,k = 0, 1, 2, . . . ,исследованный нами в §3.9г. Если переписать его в видеx(k+1) ← x(k) − τ Ax(k) − b ,k = 0, 1, 2, . . . ,(3.115)3703. Численные методы линейной алгебрыто расчёт каждой последующей итерации x(k+1) может трактоватьсякак вычитание из x(k) поправки, пропорциональной вектору текущейневязки (Ax(k) − b). Но при таком взгляде на итерационный процессможно попытаться изменять параметр τ в зависимости от шага, т.
е.взять τ = τk переменным и рассмотреть итерацииk = 0, 1, 2, . . . .(3.116)x(k+1) ← x(k) − τk Ax(k) − b ,Эту нестационарную версию метода простой итерации часто связывают с именем Л.Ф. Ричардсона, предложившего её идею ещё в 1910 году. Он, к сожалению, не смог развить удовлетворительной теории выбора параметров τk , и для решения этого вопроса потребовалось ещёнесколько десятилетий развития вычислительной математики.
Отметим, что задача об оптимальном выборе параметров τk на группе изнескольких шагов приводит к так называемым чебышёвским циклическим итерационным методам (см. [37, 45, 78]).Можно пойти по намеченному выше пути дальше, рассмотрев нестационарное обобщение итерационного процессаx(k+1) ← (I − ΛA) x(k) + Λb,k = 0, 1, 2, . . .
,который получен в результате матричного предобуславливания исходной системы линейных алгебраических уравнений. Переписав его вычислительную схему в видеx(k+1) ← x(k) − Λ Ax(k) − b ,k = 0, 1, 2, . . . ,нетрудно увидеть возможность изменения предобуславливающей матрицы Λ в зависимости от номера шага. Таким образом, приходим квесьма общей схеме нестационарных линейных итерационных процессовk = 0, 1, 2, . . . ,x(k+1) ← x(k) − Λk Ax(k) − b ,где {Λk }∞k=0 — некоторая последовательность матриц, выбор которойзависит, вообще говоря, от начального приближения x(0) .Другой популярный путь построения нестационарных итерационных методов для решения уравнений — использование вариационныхпринципов.Интуитивно понятный термин «вариация» был введён в математикуЖ.-Л.
Лагранжем для обозначения малого изменения («шевеления»)независимой переменной или рассматриваемой функции (функционала). Соответственно, метод исследования экстремумов, основанный на3.10. Нестационарные итерационные методы371изучении зависимости функций от вариаций их аргументов, получилназвание метода вариаций. Но со временем «вариационными» сталиименовать методы решения различных уравнений, которые сводят исходную постановку к тем или иным задачам на нахождение экстремума. Согласно этой терминологии, вариационными принципами теперьназывают переформулировки интересующих нас задач в виде какихлибо оптимизационных задач, т.
е. задач на нахождение минимумовили максимумов. Тогда итерационные методы решения СЛАУ могутконструироваться как итерационные процессы для отыскания этих экстремумов тех или иных функционалов.Вариационные принципы получаются весьма различными способами. Некоторые из них вытекают из содержательного (физического, механического и пр.) смысла решаемой задачи. Например, в классическоймеханике хорошо известны «принцип наименьшего действия Лагранжа», в оптике существует «принцип Ферма» [70].
В последнее столетиеимеется тенденция всё меньше связывать вариационные принципы сконкретным физическим содержанием, они становятся абстрактнымматематическим инструментом решения разнообразных задач.Строго говоря, в вычислительном отношении получающаяся в результате описанного выше сведе́ния оптимизационная задача можетбыть не вполне эквивалентна исходной, так как задача нахожденияустойчивого решения уравнения может превратиться в неустойчивуюзадачу о проверке точного равенства экстремума нулю (этот вопрос более подробно обсуждается далее в §4.2б). Но если существование решения уравнения известно априори, до того, как мы приступаем к его нахождению (например, на основе каких-либо теорем существования), товариационные методы становятся важным подспорьем практическихвычислений.
Именно такова ситуация с системами линейных алгебраических уравнений, разрешимость которых часто обеспечивается различными результатами из линейной алгебры.Как именно можно переформулировать задачу решения СЛАУ в виде оптимизационной задачи? По-видимому, простейший способ можетосновываться на том факте, что точное решение x⋆ зануляет нормуневязки kAx − bk, доставляя ей, таким образом, наименьшее возможное значение.