Введение в вычислительную физику - Федоренко Р.П.
ISBN 5-7417-0002-0
Скачать (прямая ссылка):
Pt=l + xXt, Byk = % + XAipk = + tXt%.
Подберем т таким, чтобы max | | = IP1I- Очевидно, это достига-
ется при х = I/L. В самом деле,
Pt* I -IxXtI > I-LZL = O,
max I pt I = max Pt = max (I + tXt) = I + t max Xt = I + XX1.
Спектр В устроен так: P1 может быть больше или меньше единицы в зависимости от знака X1, Pt Є [О, P1].
Теперь можно применить степенной метод к оператору В:
ui+i = But,
іі+1 = м'+1/||м'+1
P^‘+l) = (м' + 1, Ui).
Как было выяснено, и‘ -» Cip1, P^ -* P1. Можно оценить и скорость сходимости: погрешность убывает, как
CS GQ. I H-TX2
Pi I-I-TX1
I - -CtX1 — X2) И
1
Xi X2
так как обычно | X11 -c^=: L, | X21 «L.
Для того чтобы составить себе представление о скорости сходимости, обратимся опять-таки к разностному оператору Лапласа на сетке NxN. В этом ,случае, как показывают простые вычисления,
(Xi — X2)/L «к IX2HN2. Скорость сходимости невысока. Легко понять, что ее можно повысить примерно вдвое, взяв X як 2/L.
Метод обратной итерации. Запишем уравнение (1) в форме Bu = рм (где В — А~1, р = X-1) и применим степенной метод к оператору В, определяя его максимальное по модулю собственное зна-
§ 16]
ГЛАВНАЯ СПЕКТРАЛЬНАЯ ЗАДАЧА ДЛЯ КРАЕВЫХ ЗАДАЧ
197
чение. Пусть, как это часто бывает, | X11 = min | Xk \, тогда В‘и° выделяет именно то, что нас интересует. Используя для иллюстрации оператор Лапласа, находим, что спектр B = А-1 расположен на интервале [1/Х,, 0], X1 < 0, а скорость сходимости степенного метода для В определяется величиной P2Z^1. Для лапласиана это дает (в терминах § 14) X1 JX1 2«0.4, т.е. скорость сходимости практически на зависит от числа узлов сетки N.
Итак, метод обратной итерации очень эффективен (погрешность убывает, как (0.4)'), но сама стандартная итерация весьма громоздка. Поскольку явного выражения для оператора В мы не имеем, метод реализуется в форме Au1+1 = и1 и каждая итерация требует приближенного решения системы линейных уравнений с матрицей А. Если min I Xk I =*= I X11, то предварительно следует сдвинуть спектр, т.е. перейти к оператору А' = А — аЕ, подобрав соответствующее значение а. Легко понять, что итерации будут тем эффективнее, чем ближе а к X1. Однако при этом заметно осложняется решение уравнения
(А — aE)ul+l = и1, так как оператор А — аЕ приближается к вырожденному оператору А — X1E.
Вычисление второй собственной функции. В некоторых задачах представляет интерес вторая собственная функция (соответствующая точке спектра X2), а иногда и последующие. Степенной
метод позволяет вычислить ее, правда, после того, как уже вычислена первая собственная функция. Используется тот же алгоритм, но в подпространстве, ортогональном найденной собственной функции гр,. В выражении Ui = B1U0 = J ск{Pjk)^* слагаемое C2(P2)iTp2
к
будет играть главную роль в случае, если начальное приближение и0 выбрано так, что C1 = 0. Другими словами, начальное приближение должно быть ортогонально первой собственной функции Ip1 (считаем, что она уже найдена), так как C1 = (и0, Ip1).
Процесс итераций организуется так. Возьмем какую-то функцию її0. Спроецируем ее на подпространство, ортогональное функции Ip1: и0 J= и0—(и0, ^p1) -Vp1. Выполняем стандартные итерации степенного метода: ul+1 = Bu1. В результате погрешностей округления (и неточности определения Tp1) функция и1 содержит пусть маленькую, HO все же ненулевую компоненту C1Tp1, поэтому нужно проводить систематическую ортогонализацию и нормировку:
~/+1 := ~/+1 _ (-/+Ij ui+1 = ~;+1/||~«+1||.
198
ПРИБЛИЖЕННЫЕ МЕТОДЫ ВЫЧИСЛИТЕЛЬНОЙ ФИЗИКИ
[Ч. II
(Иногда ортогонализацию и нормировку проводят не на каждой итерации, а периодически, например через пять-десять итераций.) Собственное значение X^ (/-е приближение) вычисляется по той же формуле: X^ = (Aui, и‘)/(и\ и1).
Ускорение степенного метода. Один из наиболее часто употребляемых приемов ускорения степенного метода состоит в варьировании итерационного параметра т. Вычисления ведутся практически по той же самой схеме, но с переменным т:
ui + 1 = (Е + xj + lA)u'.
Разлагая и0 в ряд Фурье по собственным векторам, получаем для и1 формулу
І
и‘ = 2 С*П (1 + tA)?-
к j = 1
Здесь, конечно, тоже следует использовать идею чебышевского метода, но с некоторой модификацией. Мы заинтересованы в том, чтобы полином PJ (I + XjXk) был как можно меньше на всем спектре, за исключением точки X1.
Для реализации этой идеи необходимо располагать некоторой информацией. Прежде всего нужно иметь оценку (сверху) для максимального (по модулю) собственного числа L и оценку для второго собственного числа. Достаточно знать число /, разделяющее X1 и X2, т.е. должно быть Xk Є [—L, /] (к = 2, 3, 4, ...), X1 > /. Параметры Xj следует брать обратными к корням полинома Чебышева (некоторой назначенной степени і) на интервале {—L, I]. Эффект ускорения тем выше, чем ближе I к X2.
Есть еще один простой прием, позволяющий ускорить сходимость: сдвиг спектра. Пусть имеется некоторая оценка первого собственного числа X*. Она может быть получена на первом этапе решения задачи, т.е. является не очень точной. Уточнение осуществляется теми же итерациями, но с небольшим изменением: используется оператор со сдвинутым спектром А — XtE, т.е.
