Федоренко Р.П. Введение в вычислительную физику (1185915), страница 39
Текст из файла (страница 39)
Таким способом получается система уравнений ~~; Ха «и(«е) = ~~~„ЬЕ „Ха — Х ,'Г се а ха аа О, й = О, 1, ..., М, (11) а О а О а О 1 с „3 К(«Ф, Р,) «м'(В) е«Е. о Теперь для определения (х„) н Х мы имеем л«+ 1 уравнений (11), которым можно придать стандартную форму спектральной задачи Ах = ХСх. Опыт показал, что сравнительно небольшие значения л«, приводящие к не очень трудоемкой алгебраической проблеме собственных значений, при таком подходе позволяют получить с хорошей точноспио сравнительно большое число точек спектра дифференциального оператора (примерно Ф/2 и даже больше), Например, для уравнения Ламе такой алгоритм уже при Л«= 21 дает дня первых восьми собственных значений величины, совпадающие с 1 1В1 ГЛАВНАЯ СПЕКТРАЛЬИАЯ ЗАДАЧА ДЛЯ КРАЕВЫХ ЗАДАЧ 191 «каноннческими табличными» в девяти десятичных знаках, следующие четыре собственных значения совпадают с табличными в восьми десятичных знаках, н т.д.
При оценке трудоемкости и точности алгоритма нужно иметь в виду, что она существенно зависит от того, как вычисляются интегралы, определяющие коэффициенты матриц А и С. Наиболее эффективные результаты получаются в том случае, когда интегралы Ь „, с» „«берутся» в конечном виде и вычисления проводятся по каким-то не очень сложным формулам с очень высокой точностью.
Если же таких формул нет н приходится вычислять интегралы (а их не так уж мало, примерно 21тз) по какнм-то квадратурным формулам, трудоемкость алгоритма заметно возрастает. й 16. Главная спектральная задача для краевых задач математической фкзмкк В приложениях часто возникает задача, которую формально можно записать в простой форме.
Пусть А — линейный дифференциальный оператор, соответствующий некоторой краевой задаче для 'уравнений с частными производными. Нужно найти главное собственное число н соответствующую ему собственную функцию: Аи = 1!и. Главным собственным числом называют обычно крайнюю точку спектра, например с наибольшим значением Ве 1. Поясним суть дела, рассмотрев важный в приложениях пример — математическую модель ядерного реактора. Разумеется, мы ограничимся сравнительно простой моделью. Ядерный реактор будем представл1пь себе в виде некоторого прямоугольного тела (например, в виде трехмерноГо куба). Распределение нейтронов в реакторе описывается системой двух- групповых уравнений диффузии: дФ! — — — б»т Р! кгаб Ф, — АЫФ, + Аыфз, ! (2) — — = б»т Рз кгаб Фз — Аз»Ф» + АмФ1, дФЗ с краевыми условиями на границе Г куба Ф, =Фз = О н какими-то начальными данными. Здесь Ф1(1, х, у, з) (1 ° 1, 2) — 'функции, описывающие распределение быстрых (Ф1) и медленных (Фз)'нейтронов.
Уравнения (2) описывают нх эволюцию во времени с учетом следующих процессов: 192 НРИБЛНЖЕННЫЕ МЕТОДЫ ВЫЧИСЛИТЕЛЬНОЙ ФИЗИКИ 1Ч. и 1) диффузия (члены <ИТЮ, ягаб Фз); 2) поглощение нейтронов (члены — Аыф, и — АЗВФВ); 3) рождение быстрых нейтронов при поглощении медленных (член А,ВФ ) н наоборот (член АВ,Ф,). Коэффициенты системы 1!9, А суть некоторые функции х, у, з, определяемые физическими консгантамн материалов, из которых составлен реактор. Зздача линейна относительно Фг Ее можно записать в компактной форме: <р,=Ар, (3) где ~р = (Ф,,ФВ), А — линейный дифференциальный оиерзтор эллиптического типа, главная дифференциальная часть которого 6190 хгаб при постоянном Ю есть просго ЮЬ (свойства оператора А во многом похожи на свойства оператора Лапласа). Чтобы представить себе характер процессов, протекающих в реакторе, воспользуемся методом Фурье для решения уравнения (3).
Так как А не зависит от 1, частные решения (3) можно искать в виде ~р(1) е"'и. Подставляя в (3), имеем 1.сын .4еыи еиАН или Аи = Хи, т.е. Х должно быть собственным числом, и — соответствующей собственной функцией оператора А. Из теории эллиптических уравнений мы знаем, что имеется дискретное множество собственных значений ХВ н соответствующих собственных функций зЦ„, образующих полную систему.
Занумеруем собственные значения в порядке убывани Ке ХВ: ке 1ье- — В при е- Ф . Начальную функцию р„разложим в ряд по ф,: ре = ~~, сьзрь. Тогда сразу имеем решение: (4) ~р(г) =~ с е~'зрю Нетрудно убедиться, что при достаточно большом времени г в решении (4) выделяется главный член с наибольшим значением Ке Х: р(1) ж с~с" 'зрг Кстати, нужно понимать, какое время является большим для процесса, описываемого системой (2). О нем естественно судить по величине еоз ">', это следует из вырзжения р(с) = е" '(с,зр, + сзе~" "р + ...), Например, в некоторых реакторах 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.»~Р». Подберем т таким, чтобы шах ! р»~ = ~ ф, ~.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
