Введение в вычислительную физику - Федоренко Р.П.
ISBN 5-7417-0002-0
Скачать (прямая ссылка):
Что же реально мы имеем? Рассмотрим для простоты двумерный случай (функции зависят от х и у). Введем в зоне реактора сетку с шагом h, узлами xk = kh, ут — mh (к, т = 0, 1, ..., N) и сеточные функции ФІ^ m, Ф2к т. Тогда вместо дифференциального уравнения можно рассматривать аппроксимирующее конечно-разностное уравнение:
(второе уравнение запишется точно так же).
Обычно в памяти ЭВМ мы имеем вектор {ФІ^ т, Ф2к т} и коэффициенты Dlk+l/2 т, ..., А21к . Иногда для них нет места в памяти и приходится использовать подпрограммы, которые по индексам к, т и какой-то относительно небольшой информации о структуре реактора позволяют получить Dlk+l/2 т и остальные коэффициенты схемы. Таким образом, имея точку и, можно вычислить точку (той же размерности) Au. Это сравнительно «дешевая» операция: она требует 0(JV2) арифметических действий.
Степенной метод определения границ спектра матрицы. Ограничимся сравнительно простым, но важным в приложениях случаем, когда оператор А самосопряженный: A = A*. В этом случае все собственные значения Xk вещественны. Следующий алгоритм позволяет вычислять максимальное по модулю собственное значение и соответствующий собственный вектор. Выбирается некоторый более
7 — 1833
194
ПРИБЛИЖЕННЫЕ МЕТОДЫ ВЫЧИСЛИТЕЛЬНОЙ ФИЗИКИ
[Ч. U
или менее произвольный вектор и0 (начальное приближение). Затем производятся итерации (г — номер итерации):
Нетрудно убедиться, ЧТО при І-* 00 вектор и1 стремится к собственной функции, соответствующей max | Xk |. В самом деле,
Ui = AiU0. Пусть — собственные векторы матрицы А. Разложим
и0 в ряд по базису яр: и0 = Ji Тогда
Очевидно, что компоненты суммы, у которых I Xk I < L, стремятся к нулю и в конце концов остается только тот член, у которого
I Xk I = L (для простоты считаем, что такой член только один).
Оценим скорость сходимости. Пусть- К — номер максимального (по модулю) Xk, К — 1 — номер следующего за ним (по модулю) собственного значения. Тогда
Таким образом, и‘ состоит из слагаемого, пропорционального Ipysr, и убывающей при / -*¦ оо погрешности. Ее можно оценить по норме
Итак, погрешность убывает, как |Я^/Я;с_1|'.
Разумеется, мы неявно предполагаем, что коэффициент ск не слишком мал, т.е. что начальное приближение не слишком плохое. Плохое оно при ск = 0, однако и в этом случае метод дает правильный результат: за счет погрешностей округления в каком-то приближении обязательно появится ненулевой коэффициент ск. Ho если он слишком мал, процесс придется доводить до слишком боль-
ui+i — Au1, і = 0, 1,2,...
к
Ui = AiJ Ck^k = J CkAi^k = J Cll(Xk)iVv
к
к,
к
Обозначая L = max | Xk \, имеем
к*к V /
ГЛАВНАЯ СПЕКТРАЛЬНАЯ ЗАДАЧА ДЛЯ КРАЕВЫХ ЗАДАЧ
195
ших чисел і, растущих, впрочем, при убывании ск со скоростью логарифма. Вычислители на это не очень надеются и стараются выбрать разумное приближение, возможно более близкое к ^Pa-. На этот счет есть достаточно надежные рецепты. Если А — разностная аппроксимация дифференциального оператора, то близка к
функции ик т = (— 1)*+,п.
Теперь добавим некоторые важные детали. Так как | Xk \ 1, то
при достаточно большом і величина A1U0 обращается либо в машинный нуль, либо в машинную бесконечность. Чтобы этого избежать, в процесс добавляется нормировка і-го приближения, после чего он выглядит так:
ui + l=Au‘, u, + l == м/ + 1/||м' + 1||.
При этом г-е приближение к собственному значению X^= (Aui, и1)/(и‘, и‘)
(если и1 = то формула дает = Xk). Критерием достигнутой точности служит соотношение
WAui - X^ui У < є,
где ? — заданная погрешность.
Полезно представить себе характерное значение скорости сходимости. Возьмем в качестве оператора А разностную аппроксимацию оператора Лапласа на сетке NxN. (Легко проверить, что суть дела не в шаге сетки h, а в числе узлов, приходящихся на линейный размер области; поэтому можно рассматривать задачу в квадрате 1x1.) Тогда
max
Xt\ =SN2 sin2 л ^SN2Il
Следующее за ним по модулю собственное значение есть 4N2 (^sin2 ?=± л + Sin2 J «
W 1--*L +1-4) =SA^fl-1-4
4N2 N2 J I 8 JV2
Итак,
~ і _ІіЕІ Xk 8 N2
Видно, что при больших N скорость сходимости степенного метода невелика.
Обратим внимание на то, что при исследовании спектра разностной аппроксимации дифференциального оператора вычисление
196
ПРИБЛИЖЕННЫЕ МЕТОДЫ ВЫЧИСЛИТЕЛЬНОЙ ФИЗИКИ
[Ч. II
max I Xt I не очень интересно. При N -*¦ оо эта величина стремится к бесконечности и особого смысла не имеет, хотя ее значение (или хотя бы ее оценка) нам понадобится. Интересной величиной является не max I Xk |, a max Xk. Спектры эллиптического дифференциального и разностного операторов устроены примерно так:
X1^ Xj ^ ... ^ Xk ... ,
• 00 при k -*¦ 00.
Нас интересует именно крайняя правая точка спектра. Итак, спектр А расположен на [—L, X1], причем JL-гж» | X11, где X1 может иметь любой знак.
Построим простой оператор с теми же собственными векторами, что и А, но с другим спектром: B=E + хА. Очевидно, его собственные векторы — те же Tpt, а собственные числа суть
