Fedorenko-RP-Vvedenie-v-vychislitelnuyu-fiziku (810773), страница 38
Текст из файла (страница 38)
О нем естественно судить по величине еоз ">', это следует из вырзжения р(с) = е" '(с,зр, + сзе~" "р + ...), Например, в некоторых реакторах 1 — Х, ж — (50 Ф 100), т.е. время 0.1 с уже является очень большим. Что же происходит с реактором? Все зависит от значения Х,: если Х, > О, реактор «взрывается», если 3 181 гл»аи»я спектг»льн»я з»д»чл для кг»звых з»д»ч 193 Х/ < О, реактор»тухнет». Рабочий режим реактора — это ситуация Х| = О. Разумеется, значение Х1 зависит от коэффициентов системы, т.е.
от физического состава реактора. Он поддается регулированию с помощью стержней, Цель этого регулирования — обеспечить значе- ние Х/ =О. Теперь понятно, почему в практике расчетов ядерных реакторов одной из главных вычислительных задач является вычисление край- ней точки спектра линейного дифференциального оператора (1).
На практике А — это не дифференциальный оператор, а конечно-раз- ностный, если мы решаем задачу (2) методом сеток, Для нас важно следующее обстоятельство: размерность конечномерного пространст- ва и очень велика (порядка 103 —: 103). Поэтому матрицы А мы обычно в явном виде не имеем. Что же реально мы имеем? Рассмотрим для простоты двумерный случай (функции зависят от х и у). Введем в зоне реактора сетку с шагом //, узлами х = /с//, у,„= т// (», л/ = О, 1, ..., л/) и сеточные функции Ф1, Ф2 .
Тогда вместо дифференциального уравне- ния можно рассматривать аппроксимирующее конечно-разностное уравнение: Ф1, — Ф1, Ф1», Ф/ к — 1. /~ Р1»ь1/2, а // Р1» п2, и» + Ф1, ч — Ф1, Ф1» — Ф!„1) +» ~Р1», +пг» Р1», -пз» вЂ” А11» Ф1» „, + А12» Ф2» — — — Ф1/,,„ », ~ ю, (второе уравнение запишется точно так же). Обычно в памяти ЭВМ мы имеем вектор (Ф1», Ф2» ) и коэффициенты Р1»~1/3,„, ..., А21»,„.
Иногда для них нет места в памяти и приходится использовать подпрограммы, которые по индексам //, л/ и какой-то относительно небольшой информации о структуре реактора позволяют получить Р1» пз и остальные коэффициенты схемы. Таким образом, имея точку и, можно вычислить точку (той же размерности) Аи. Это сравнительно «дешевая» операцию она требует О(/з/3) арифметических действий.
Степенной метод определения границ спектра матрицы. Ограничимся сравнительно простым, но важным в приложениях случаем, когда оператор А самосопряженный: А = А'. В этом случае все собственные значения Х» вешественны. Следующий алгоритм позволяет вычислять максимальное по модулю собственное значение и соответствующий собственный вектор.
Выбирается некоторый более 7 — 1833 пгивлиженные методы вычислительной еизикн [ч,п ! 94 или менее произвольный вектор ив (начзльное приближение). За- тем производятся итерации (! — номер итерации): и'+'=Аи', з=О, 1, 2, Нетрудно убедиться, что прн !- м вектор и' стремится к собственной функции, соответствующей шзх ! Х !. В самом деле, и' = Агин. Пусть 1р — собственные векторы матрицы А. Разложим ив в ряд по базису цк ив ~; с~~ы Тогда и' А' ,'~~ с„ц> = ~ с А'1р = ~ ск(Х~)' р .
Обознзчая к,=шах !Ль), имеем Очевидно, что компоненты суммы, у которых ! Хк! < Ь, стремятся к нулю и в конце концов остается только тот член, у которого ! Хь ! = к, (для простоты считаем, что такой член только один). Оценим скорость сходимости. Пусте К вЂ” номер максимального (по модулю) Х, К вЂ” 1 — номер следующего за ним (по модулю) собственного значения. Тогда Р,~ " = ()ьк)' с»Же+а са ~~ ~ 'рв ° а к Таким образом, и' состоит из слагаемого, пропорционального 1р», и убывающей при !-ь м погрешности. Ее можно оценить по норме Итак, погрешность убывает, как ! Хк/!ьк Разумеется, мы неявно предполагаем, что коэффициент ск не слишком мал, т.е.
что начальное приближение не слишком плохое. Плохое оно при с = О, однако и в этом случае метод дает правильный результат: за счет погрешностей округления е каком-то приближении обязательно появится ненулевой коэффициент ск. Но если'он слишком мзл, процесс' придется доводить до слишком боль- В !61 ГЛАВНАЯ СПЕКТРАЛЬНАЯ ЗАДАЧА ДЛЯ КРАЕВЫХ ЗА!!АЧ !95 ших чисел !, растущих, впрочем, при убывании с со скоростью логарифма. Вычислители на зто не очень надеются и стараются выбрать разумное приближение, возможно более близкое к тря. На этот счет есть достаточно надежные рецепты.
Если А — разностная аппроксимация дифференциального оператора, то 1р„близка к функции ик = ( — 1)" Р Теперь добавим некоторые важные детали. Так как ~ 2к~ ~ 1, то при достаточно большом ! величина А'ПВ обращается либо в машинный нуль, либо в машинную бесконечность. Чтобы этого избежать, в процесс добавляется нормировка !-го приближения, после чего он выглядит так: и2+! =Аи', и!+! = иР+!1!!иРР!'!!, При этом 2-е приближение к собственному значению 2!!!! = (Аи!, и!)/(и!, и!) (если и! = Ч',, то формула дает А<!! = Хк), Критерием достигнутой точности служит соотношение !!Аи! — Л<!!и!!! < е, где е — заданная погрешность, Полезно представить себе характерное значение скорости сходимости.
Возьмем в качестве оператора А разностную аппроксимацию оператора Лапласа на сетке 2!! х 2!!. (Легко проверить, что суть дела не в шаге сетки 2!, а в числе узлов, приходящихся на линейный размер области; поэтому можно рассматривать задачу в квадрате 1 х 1.) Тогда шах !А ! =8Д! Я1п — л 8Д! 1 — — . г ЗА! — ! г К = 2А! ~ З). Следующее за ним по модулю собственное значение есть 4!!! Гйп г/ гА' — ! г А! — 2 2!Ч л+зш — л 2А! 4дг 1 — — "+1 — — ", =8ЛЯ 1 — 5Я Итак, 2 К-! !Я ж! — —— В„г Видно, что при больших Ф скорость сходимостн степенного метода невелика.
Обратим внимание на то, что при исследовании спектра разностной аппроксимации дифференциального оператора вычисление Т* ПРИБЛИЖЕННЫЕ МЕГОДЫ ВЫЧИСЛИТЕЛЬНОЙ ФИЗИКИ (ч. и шах ! Х ! не очень интересно. При ЛР- эта величина стремится к бесконечности и особого смысла не имеет, хотя ее значение (или хотя бы ее оценка) нам понадобится. Интересной величиной является не шах ~ Х 1, а пзах Х», Спектры эллиптического дифференциального и разностного операторов устроены примерно так: л,>Х » ...
Х ..., Х» — пРИ /Р Нас интересует именно крайняя правая точка спектра. Итак, спектр А расположен на 1 — Е, Х,), причем /.м ) Х, ~, где Х, может иметь любой знак. Построим простой оператор с теми же собственными векторами, что и А, но с другим спектром: В = В+ ЕА. Очевидно, его собственные векторы — те же 1р, а собственные числа суть р» = 1 + ТХ», Вф» = ф~ + ТА1Р» = тр» + Т1.»~Р». Подберем т таким, чтобы шах ! р»~ = ~ ф, ~. Очевидно, это достига- ется при т = 1/В. В самом деле, р» > 1 — ~ ЗЛ» ~ > 1 — »,//. = О, шах ! р ! = шах р = шах (1 + ЗХ ) = 1 + т шах Л = 1 + ЕХН Спектр В устроен так: ~, может быть больше илн меньше единицы в зависимости от знака Х„ ~» е [О, р,).
Теперь можно применить степенной метод к оператору В: и'+' = Ви~, и'+' = и'+'Ци'+'~~, ру+П = (и~Ф', и'), Как было выЯснено, и'- С~РН Р1Й- Дг Можно оЦенить н скоРость сходимостн: погрешность убывает, как р ~+,~ ~1 т("ч»з) ! так как обычно ~ Х, ~ «2„~ Х ~ « /.. Для того чтобы составить себе представление о скорости сходи- мости, обратимся опять-таки к разносгному оператору Лапласа на сетке М х л/. В атом,случае, как показывают простые вычисления, (Х, — Хз)// - пз/2Ф. Скорость сходямости невысока.
Легко понять, что ее можно повысить примерно вдвое, взяв т 2//.. Метод обратной итерации. Запишем уравнение (1) в форме Ви = ри (где В = А ', р = » .') н применим степенной метод к оператору В, определяя его максимальное по модулю собственное зна- з 161 ГЛАВНАЯ СПЕКТРАЛЬНАЯ ЗАДАЧА ДЛЯ КРАЕВЫХ ЗАДАЧ 197 чепце. Пусть, как это часто бывает, [А,[ =ппп [)ь [, тогда В1ие выделяет именно то, что нас интересует. Используя для иллюстрации оператор Лапласа, находим, что спектр В = /ь ' расположен на интервале [1/Д„01, А1 < О, а скорость сходимости степенного метода для В определяется величиной р /~г Для лапласиана это дает (в терминах 5 14) Д1,/А, ЕЛО.4, т.е. скорость сходимости практически на зависит от числа узлов сетки /Ь/. Итак, метод обратной итерации очень эффективен (погрешность убывает, как (0.4) '), но сама стандартная итерация весьма громоздка.
Поскольку явною выражения для оператора В мы не имеем, метод реализуется в форме Аи' Ь ' = и' и каждая итерация требует приближенного решения системы линейных уравнений с матрицей А. Если ппп [ 1ьь[ ~ [ 1ь1 [, то предварительно следует сдвинуть спектр, т.е. перейти к оператору А' = А — аЕ, подобрав соответствующее значение а.
Легко понять, что итерации будут тем эффективнее, чем ближе а к Хг Однако при этом заметно осложняется решение уравнения (А — аЕ)и'+' = и', так как оператор А — аЕ приближается к вырожденному оператору А — 1ь1Е. Вычисление второй собственной функции. В некоторых задачах представляет интерес вторая собственная функция (соответствующая точке спектра 1ь ), а иногда и последующие.
Степенной метод позволяет вычислить ее, правда, после того, как уже вычислена первая собственная функция. Используется тот же алгоритм, но в подпространстве, ортогональном найденной собственной функции рг В выражении и'=В'из ~х' с„ф )' р слагаемое с (р )'1р будет играть главную роль в случае, если начальное приближение ие выбрано так, что с, = О. Другими словами, начальное приближенис должно быть ортогонально первой собственной функции 1р, (счнтаем, что она уже найдена), так как с, = (ие, 1р,).
П.ооцесс итераций организуется так. Возьмем какую-то функцию из. Сгроецнруем ее на надпространство, ортогональное функции р,: ие:= и~ — (ие, ~р,)~р1. Выполняем стандартные итерации степенного метода: и'+' = Ви'. В результате погрешностей округления (н неточности определения зр,) функция и' содЕржит пусть маленькую, но все же ненулевую компоненту с, ~р,, поэтому нужно проводить систематическую ортогоналнзацию и нормировку: и' 1=и1 — (и, зр~)зр1, и~~~ = и~ ~/[[и1+'[[.