Fedorenko-RP-Vvedenie-v-vychislitelnuyu-fiziku (810773), страница 33
Текст из файла (страница 33)
если Яю= — хзг, то ( — 51«, «) = = Ци, 1«), т.е. х е (У„Уз). Сведя исследование итерационного процесса с регуляризатором В к исследованию простой итерации с оператором 5, можно воспользоваться уже знакомой нам теорией. В частности, если границы у, и у близки друг к другу (т.е. оператор В 'Р «хорошо обусловлен»), то в качестве оптимального итерационного параметра можно взять т = 2/(у, + уз), что приведет к сходимости с множителем (уз — т,)/( у, + уз) за одну итерацию. Правда, не следует забывать, что «цена» такой итерации зависит от затрат вычислительной работы на решение уравнения Ве = г.
Если оператор В 'Р «плохо обусловлен (у /у,з~!), то можно использовать чебышеаское ускорение и получить процесс, в котором средняя за итерацию эффективность соответствует множителю (1 — 2 ~7;/Я). Применение общей схемы. Конкретные итерационные алгоритмы получаются при конкретном выборе регуляризатора В.
Приведем некоторые примеры. Метод простой итерации. Он получается при В = Е. Метод переменных направлений (точнее, некоторые его обобщения). Он получается при выборе в качестве регуляризатора «факто- ризованной» конструкции В»в Š— а — Š— оз— [ч. 1 основы вычнслнтхльной млтемлтнкн 1то (Напомним, что мы договорились производные понимать в смысле их простейших разностных аппроксимаций.) Оператор В легко «обращается». Решение уравнения Во = з сводится к решению последовательности задач: а) Š— а — н'=г, б) Š— аз — о=о'.
! хг) ' ~ г) Каждая из этих задач расщепляется на серии «одномерных», легко решаемых прогонкой систем уравнений. Здесь мы сталкиваемся с характерной ситуацией: конструкция оператора В содержит некоторые параметры (а„а в данном случае). Вместе с параметром т они должны быть найдены таким образом, чтобы получить возможно более высокую скорость сходимости. Точный анализ сходнмости в общем случае провести не удается. Поэтому задача «оптимизации параметров» обычно решается раздельно: сначала за счет выбора параметров регулярнзатора стремятся уменьшить обусловленность 5 — величину уз/у,, т.е. сблизить оценки в неравенствах энергетической эквивалентности операторов В и — Р. Теоретической предпосылкой для существования хороших оценок такого типа является известный факт из теории эллиптических операторов; два любых эллиптических дифференциальных оператора одного порядка энергетически эквивалентны друг другу.
Следствием этого является и энергетическая эквивалентность их разностных аппроксимаций с постоянными ун у, не зависящими по существу от шага сетки Л. Рассматриваемый здесь факторизованный оператор В является, как легко заметить, аппроксимацией дифференциального оператора четвертого порядка (правда, вырожденного). Дифференциальные операторы разных порядков не могут быть энергетически эквивалентными.
Это приводит к существенной зависимости констант у, от Л: у /у, = О(Л ') при «оптимальном» выборе ан а. Попеременно-треугольный метод. В качестве регуляризатора используется факторизованный оператор В=Я1йз Е+ а — „+ — д Š— а — „+— д д1З д д 2 Очевидно, В= Š— аз~ — + — ~. Оператор ~ — + — ) имеет второй дх ду~ ' дх ду порядок, но, к сожалея)по, не является строй эллиптическим. Сла- гаемое Е при подходящем выборе о. придает разностной аппроксимации В «эллиптический» характер.
Оператор В легко обратим: решение уравнения Ви= г требует числа операций, пропорционального числу неизвестных, т.е. числу з! 41 гашение эллиптнчвскнх э»д«ч методом свток узлов сетки. Чтобы убедиться в этом, выпишем подробно разностную аппроксимацию Я, и Рз: «! На рис. 16 показана сетка и расположение в ней шаблонов операторов Я,, Яз; такая аппроксимация обеспечивает важное свойство Я =Я; Из этого рисунка очевиден алгоритм «обращеннй» Вп Я при нзвестнык значениях о на границе. Решение, например, уравнения Я,о = з осуществляется «маршевым» алгоритмом вычисления слева-направо и снизу-вверх, начиная с левою нижнего угла области.
На каждом шаге такого алгоритма в выражении (М|о)» из трех значений о два (о«, и о„ ,) и уже известны н можно вычислить 1 Г о 1+2«ул ~ кч л ( «-к ««,т-1) Устойчивость этого «марша» легко устанавливается (прн о > О). Обращение йз осуществляется аналогично, но в обратном направлении (начиная с правого верхнего угла области). Оптимизация оценок эквивалентности за счет о приводит к уз/у, = О(л ').
Выбирая далее последовательность итерационных параметров т,. в соответствии с теорией чебышевского ускорения, получаем алгоритм, в котором число итераций, необходимых для уменьшения погрешности начального приближения в е ' раз, есть (для сетки ЖхуУ) У(е) = О(~/У!и е 1). Ограничимся здесь этими общими сведениями, отправляя интересующихся деталями (как оценивать у„уз, как выбирать и и т.п.) к специальной литературе. «Двухступенчатые» итерационные методы.
Почти очевидно, в каком направлении следует искать операторы В, наилучшие с точки зрения оценок энергетической эквивалентности: оператор В должен быть возможно более похожим на — Уэ (идеальный случай: В = — Ю, у, = уз »» 1„ достаточно одной итерации; к сожалению, она просто эквивалентна решению исходной задачи). Однако по мере сближения В с — Р возрастают трудности решения уравнения Во = г.
основы вычислительной млтемлтики !ч.1 Удачным компромиссом является, например, выбор В этом случае у~у, ие зависит от л, а для решения разиостиого уравнения — Ьн= з можно использовать построенные в последние годы эффективные алгоритмы для решения в прямоугольной области уравнения с постоянными коэффициентами. Ограничимся здесь только названиями методов, тем более что многие из иих уже оформлеиы как стаидартиые быстро работающие программы математического обеспечения современных ЭВМ: это методы циклической редукции, быстрого преобразоваиия Фурье (см. 5 24), маршевый и иекоторые другие.
Кстати, включение некоторых из иих в арсенал средств практических вычислений было связаио с анализом вычислительной неустойчивости и разработкой вычислительно устойчивых модификаций (как это было с чебышевским ускорением). Простейшие формы алгоритмов неустойчивы. Можно сближать В с — ет, решая уравнение Ви = х подходящим итерационным методом, например методом Переменных иаправлений. При этом, естественно, иет необходимости в очень точном решении уравнения, достаточно ограничиться каким-то числом «виутренних» итераций. Так мы приходим к семейству «двухступенчатых» итерационных алгоритмов.
Их оптимизация, в частности выбор наилучшего (с точки зрения эффективности процесса в целом) числа внутренних итераций, связана с достаточно сложиым анализом. Такие итерационные процессы д теория их оптимизации построены Е. Г. Дьяконовым. Миогосеточиый метод. Опишем конструкцию своеобразного итерационного алгоритма, получившего в последние годы широкое применение по причинам, которые удобно обьяснить несколько позже. Метод достаточно сложен алгоритмически, сложен ои и для теоретического анализа даже в самом простом случае. Поэтому мы ограиичимся самым общим описанием и качественным обьясиеиием механизма, обеспечивающего высокую скорость сходимости итераций. Исходной идеей этой конструкции является следующее замечаиие.
Все собственные функции оператора о (разиостиого) зГ « = з1п — зш — условно разделим иа две части: гладкие крв . Етв , » ж л ( р < ЛЧ2 и 4 < йЧ2) и негладкие (р з 'ЛЧ2 или д ~ лЧ2) функции, т.е. разделим низкие и высокие гармоники и частоты Х некоторой условной границей. Легко построить итерациоиный метод, эффективно гасящий иегладкую компоненту погрешности (и иевязки).
В самом деле, метод простой итерации и' = и'+ т(Ьи' — ~), как было показано, гасит (р, д)-компоиеиту погрешноСти, умножая ее за один шаг на 1тз 8 141 гашение эллиптических э»д»ч методом снох 1 — тХ . Высокие частоты расположены, очевидно, между Х, ли —— = (4/Лз) зшз(пН/4Я) 2/Л» и Х„,, „,, 8/Лз. Выбирая т оптимальным для этой части спектра, т.е. т = 2Л»/(2 + 8) = 0.2Л», получаем убывание негладких компонент погрешности (невязки) в процессе итераций со скоростью (О.б) ', т.е.
достаточно быстрое. На остальной части спектра сходимость, конечно, очень медленная. Так, компонента (1, 1) погрешности убывает с показателем 1 — 2пэт 1 — 0.4(п/Ф)» за шаг. В целом итерационный процесс оказывается неэффективным. Однако неболыпое число таких итераций «сглаживает» невязку: высокие гармоники в ней гасятся, основную роль играет гладкая компонента. Таким образом, после 1 итераций имеем приближение и', удовлетворяющее уравнениям (А и')» — /» — — г,', (внутри) и' — р = 0 (на границе) (это просто определение невязки г). Если бы мы могли решить уравнения (Ьн) = г' (внутрн) и»,„= 0 (на границе) то функции ю была бы поправкой в том смысле, что точное решение и» = и» вЂ” е» .
На первый взгляд, найти поправку ю — задача такой же степени трудности, как и исходная. Но после /итераций с т = 0.2Л» ситуация изменилась: невязка г' стала гладкой функцией и уравнение для поправки можно решать на другой, более грубой сетке. Предположим для простоты изложения, что М= 2~, и наряду с основной сеткой с шагом Л = 1/М введем вспомогательную сетку с шагом Н= 2Л.
Узлы этой сетки совпадают с четными узлами основной сетки (т.е. с теми узлами, (Л, »л), для которых четны оба индекса Л н лч). На этих сетках мы будем рассматривать близкие по смыслу функции — в этих случаях будем использовать одинаковые буквы для Н-сеткн и Л-сетки (большие и малые соответственно). Итак, имеем 1'-е приближение и' и его невязку г' . Возьмем ограничение невязки на Н-сетку, т е., проще говоря, А» —— г „ (Л, т = О, 1, ..., Н/2), и решим на Н-сетке уравнение (ВЖ)» = Р 1ч.1 174 ОСНОВЫ ВЫЧИСЛИТЕЛЬНОИ МАТЕМАТИКИ Теперь вопрос в следующем: из того, что (17Вг) „= г ь 7, следует ли (хотя бы приближенно), что (йв)4 и г' и? Если да, то проведенная коррекция явно целесообразна. Вычисления (простые, но довольно о о О о ! о о о О о о о о о о о о о О о г 6 о ъ о Г Асс'-н') Ь Ь Ьг В Рис, 17 громоздкие) показывают, что ответ неоднозначен.
Точнее, соотношение Ь м ж гс выполняется «в среднем», в слабом смысле слова. Поясним это на одномерном примере. На рис. 17 показаны этапы процесса: а — гладкая на 7ь-сетке невязка г', б — ее ограничение (с нулевыми значениями на границе). здесь 4у — аппроксимация оператора Лапласа на Н-сетке.
Эта задача заметно проще исходной хотя бы потому, что в ней в четыре раза меньше неизвестных и число обусловленности для системы меньше (тоже в четыре раза). Тем не менее решение такой системы не настолько проще решения исходной задачи на й-сетке, чтобы молсно было пренебречь проблемой решения вспомогательной зщ1ачи.
Пока будем считать, что вспомогательная задача так или иначе решена (приближенно, вообще говоря), Интерполируя (линейно по х и у, например) функцию 1т на узлы основной й-сетки, получаем функцию ье . Вычитая ее из и', приходим к новому приближению й = и' — и. Что можно сказать о невязке этой функции? Вычислим ее во внутренних узлах: г = (сь(и' — ь»)) — / = г, — (с»и) 1 те з 14) гашение эллиптических злллч методом сеток иа Н-сетку Я; в — функция г = Ь(и' — и) . Обратим внимание на характер функции г: она в среднем близка к нулю и состоит, в основном, из негладких собственных функций.