Формализация вычислений при численном решении краевых задач Текст научной статьи по специальности «Математика»

УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА Том 150, кн. 1 Физико-математические пауки 2008

УДК 519.632.4

ФОРМАЛИЗАЦИЯ ВЫЧИСЛЕНИЙ ПРИ ЧИСЛЕННОМ РЕШЕНИИ КРАЕВЫХ ЗАДАЧ

С. И. Мартыненко

Аннотация

Рассматривается один из подходов к решению проблемы формализации вычислений при численном решении широкого класса краевых задач со скоростью сходимости, близкой к оптимальной. Дано описание основных компонент универсальной мпогосеточпой технологии и результатов вычислительных экспериментов.

Ключевые слова: краевая задача, мпогосеточпые методы, метод конечных элементов. вычислительный эксперимент.

Введение

Начиная с середины 1980-х годов габариты и стоимость вычислительной техники стали стремительно уменьшаться, а производительность стала стремительно нарастать. К компьютерам получили массовый доступ инженеры, физики, химики и другие специалисты, которые не имели достаточной подготовки в области вычислительной математики, но у которых были свои задачи, зачастую связанные с необходимостью численного решения дифференциальных уравнений в частных производных. Сразу возникла проблема эффективного использования компьютерной техники пользователями, которые не обладают достаточными навыками программирования н познаниями в области численных методов.

Выход был найден в разработке автономных1 программ, используя которые инженер только формулирует задачу, а детали вычислительного алгоритма ему могут быть даже неизвестны. Применительно к техническим приложениям, работу автономных программ упрощенно можно представить следующим образом: конструктор проектирует некоторый узел при помощи графической программы, например AutoCAD. Затем геометрия узла переносится в вычислительный модуль, конструктор задает граничные и начальные условия, после чего проводит тепловой, прочностной, гидродинамический или другой расчет и анализирует результаты. Как правило, после анализа полученных результатов нужно внести изменения в конструкцию н повторить расчет. Обычно конструктор выполняет несколько подобных «итераций», чтобы получить оптимальную, с его точки зрения, конструкцию. Еще большую практическую ценность будет представлять возможность расчета машины в целом, а не только отдельных ее узлов, поскольку поэлементное моделирование часто связано с погрешностями постановки граничных условий.

Однако широкому внедрению методологии математического моделирования в машиностроение препятствует несовершенство численных методов. Применение традиционных численных методов подразумевает возможность контроля математиком хода вычислительного процесса и внесения изменений в написанную им же

англоязычной литературе, помимо термина autonomous software, иногда используется термин black box software, то есть программы, устроенные по принципу «черного ящика» [1].

программу. Ранее такой подход был оправдан, но с появлением мощных персональных компьютеров возникла острая потребность в новых вычислительных методах, обладающих достаточной точностью, эффективностью, высоким уровнем формализации и параллелизма без контроля со стороны пользователя. Уже сейчас такие программные продукты, как ANSYS, STAR-CD, FLUENT, CFX, PHOENICS и др., получили широкое применение в НИИ и ОКБ для решения разнообразных прикладных задач. Конечно, подобные программные средства еще далеки от совершенства и часто вызывают обоснованные претензии пользователей, что обусловлено использованием в них традиционных численных методов. Однако на автономные продукты существует устойчивый спрос, а значит, будет и предложение. Конкуренция на рынке программного обеспечения неизбежно приведет к повышению их потребительских качеств.

Сформулируем задачу о построении универсального алгоритма для автономных программных продуктов. Пусть дана исходная дифференциальная задача

£(u) = f (ж),

где ж G Q, Q - область N-мерного пространства, f (x) - заданная функция, L -(не)линейный дифференциальный оператор. Полагается, что граничные условия учтены оператором L и правой частыо f(ж). В области Q строится сетка и/или триангуляция, на которой осуществляется контрольно-объемная и/или конечноэлементная аппроксимация дифференциальной задачи. Обозначим результирующую СЛАУ, полученную в результате аппроксимации, как

A x = b.

Предположим, что имеется множество итерационных методов M, при помощи которых можно решить полученную СЛАУ. Тогда задача о построении универсального алгоритма состоит в отыскании такой стратегии применения итерационных методов из множества M к решению СЛАУ A x = b, которая позволит достичь скорости сходимости, близкой к оптимальной, то есть для численного решения краевых задач нужно выполнить O(N lgk N) арифметических операций, где N есть число неизвестных, a k = 1, 2.

При этом на универсальный алгоритм накладываются следующие ограничения:

1. Минимальное использование ресурсов компьютера. Темпы развития современных математических моделей опережают совершенствование аппаратной части вычислительной техники. Поэтому в отдельных случаях ресурсы компьютера могут быть задействованы преимущественно для хранения разностной задачи, а не промежуточных вычислений. Минимальный объем хранимых данных составит ^ ~ 2N вещественных чисел, то есть в оперативной памяти могут размещаться толь-

xb

A

после каждой итерации по нелинейности. Рациональное использование ресурсов компьютера позволит применять достаточно мелкие вычислительные сетки для численного решения краевых задач в автономных программах.

2. Возможность изметния способа и/или порядка аппроксимации в процессе решения задачи.

3. Отсутствие проблемно-зависимых компонент. В универсальном алгоритме не допускается применение компонент, которые зависят от решаемой задачи и обеспечивают высокую скорость сходимости. Однако возможно использование проблемно-зависимых параметров (например, верхней и/или нижней релаксации) при условии, что существует эффективная методика определения их оптимальных значений без контроля со стороны пользователя.

4. Возможность эффективного распараллеливания вычислений вне зависимости от используемых итерационных методов из множества М.

Повидимому. отсутствие проблемно-зависимых компонент в универсальном алгоритме не позволит достичь оптимальной скорости сходимости (к = 0), хотя этот факт не доказан.

Очевидна невозможность применения только одного итерационного метода в автономных программах, поэтому нужно научиться комбинировать различные алгоритмы для получения универсальной и высокоэффективной вычислительной технологии. Для этого понадобилась принципиально новая идея, которую предложил выдающийся отечественный математик Радий Петрович Федоренко. В 1961 и 1964 гг. им были опубликованы две работы, в которых использовалась высокая сходимость некоторых итерационных методов для высокочастотных гармоник [2. 3]. В 1966 г. Н.С. Бахвалов доказал оптимальность метода по числу арифметических операций для достижения точности, согласованной с порядком аппроксимации [4]. Миосеточные методы для МКЭ были изложены в монографиях В.В. Шайдурова [5] и М. А. Ольшанского [6].

Перспектива решать краевые задачи с оптимальными вычислительными усилиями выглядела очень заманчиво, но развитие многосеточных методов пошло по традиционному для 1980-х годов пут адаптации отдельных компонент алгоритма к решаемой задаче. Достаточно быстро удалось разработать высокоэффективные многосеточные алгоритмы для решения уравнения Пуассона на равномерных сетках, однако усложнение решаемых краевых задач (нелинейность, анизотропия, разрывность коэффициентов и т. д.) быстро превратили классические многосеточные методы (КуМ) в труднообозримое семейство алгоритмов, практически не поддающееся формализации. Из-за отсутствия четких критериев оптимальности адаптации компонент К\]М к конкретной задаче трудно судить о достижимости оптимальной скорости сходимости. Сразу были предприняты многочисленные попытки разработать универсальный вариант КуМ для решения широкого класса прикладных задач путем последовательного улучшения отдельных компонент [7. 81.

Хотя до сих пор не прекращаются попытки разработать универсальный вариант К\]М, стала очевидной необходимость поиска иной формы алгоритмизации основополагающей идеи Р.П. Федоренко для последующей реализации в автономных программах. Универсальная Многосеточная Технология (УуТ) определяется как наиболее рациональная организация вычислений при решении краевых задач без контроля хода вычислительного процесса пользователем автономной программы.

1. Универсальная многосеточная технология

Для построения универсального алгоритма представляется перспективной взаимная адаптация исходной дифференциальной задачи, вычислительной сетки, способа аппроксимации, итерационных методов решения СЛАУ и архитектуры компьютера. Поэтому УуТ2 состоит из двух частей: аналитической (адаптация краевой задачи к численным методам) и вычислительной (аппроксимация адаптированной краевой задачи и решение результирующей СЛАУ при помощи оригинального многосеточного метода).

2Термин «технология» используется в традиционном смысле, а именно как «совокупность приемов обработки», указывая па наличие аналитической и вычислительной частей.

1.1. Аналитическая часть технологии. Для наглядности рассмотрим краевую задачу вида

для иллюстрации основных компонент У\,]Т. Аналитическая часть У\,]Т для краевых задач типа (1) состоит в представлении искомого решения м(ж) в виде суммы двух функций е(ж) и й(ж), то есть

В последующих многосеточных итерациях3 сеточный аналог функции и(ж) будет служить приближением к решению разностной краевой задачи, а сеточный аналог функции с(ж) - поправкой, вычисляемой па грубых сетках. Представление (2), называемое £-модификацией реш ения и (ж), является одной из форм адаптации решаемых краевых задач к У\]Т.

Подстановка представления (2) в (1) приводит к следующей Е-модифициро-ваниой форме исходной краевой задачи

В (3) члены с поправкой c(x) переносятся в левую часть, а остальные - в правую.

Е-модификация решения, используемая в УуТ, похожа на представление решения, которое применяется в КмМ, но между между ними существует два принципиальных отличия:

• Е-модификация в УуТ осуществляется перед дискретизацией исходной краевой задачи для более точной формулировки разностных задач на грубых сетках и возможности гибкого изменения типа п/пли порядка аппроксимации;

•Е

дач к УуТ. Представление искомого решения в виде произведения двух функций (так называемая П-модификация) может оказаться более предпочтительной для ряда нелинейных задач [9].

ЕП

кативиого выделения особенностей, соответственно.

1.2. Вычислительная часть технологии. Вычислительная часть У^]Т связана с компьютерным счетом и состоит из построения мелкой и грубых сеток, аппроксимации модифицированной краевой задачи иитегро-иитерполяциоииым методом и решения полученных сеточных уравнения при помощи унифицированного многосеточного метода.

3Данные названия функций и (approximation to the solution) и c (coarse grid correction) пришлось позаимствовать из K\]\l для сохранения устоявшейся терминологии.

(1)

О < x < 1, м(О) = uq, м(1) = мі, k(x) ^ а > О, q(x) ^ О

u(x) = c(x) + u(x).

(2)

^О) = mq — м(О), c(1) = м1 — м(1), (За)

где правая часть r(x) имеет вид

(ЗЬ)

£ £ £ £ £ £ £ £ £

X Т '7* 7* /|* /г» <7» гр О™

^ иЬ 2 «я/д а/^ о/д 0/0 О/у «£/д

«.V „.V V -V V „V „V „V «.V „.V

А -^1 х2 х3 х4 х5 х6 х7 х8 х9 х10

СГг I « I • I I • I • I I • I • I I • I

С1

Самая мелкая сетка

у Первая сетка первого уровня

Ч---------------

с;

С\

I • I«I

I« I ■ I

I»I ■ I

Самая мелкая сетка

Вторая сетка первого уровня

01

0\

-ЬН-

I • I • I

I • I • I

Самая мелкая сетка

I • I • I

Третья сетка первого уровня

Рис. 1. Построение грубых сеток в универсальной мпогосеточпой технологии

1.2.1. Построение самой мелкой сетки. Первый этап вычислительной части УуТ состоит в построении самой мелкой сетки в области И = [0,1] для последующей аппроксимации модифицированной краевой задачи интегро-интерполяционным методом. Сетка С0 состоит го двух множеств точек (0;1) и С (0;1), которые задаются следующими соотношениями:

(0;1) = {хГ1 хГ = н(*-1)> * =1 2,

С(0;1) = {х1 | х1 = 0.5 (х! + ж!+1), *

н0

• , НХ

1,2,.

+ 1, н = (НХ)-1},

• •, нХ}.

В общем случае до получения разностного аналога £-модифицированной краевой задачи будем называть х! и х^ точками сетки. В КуМ построение грубых сеток зависит от способа аппроксимации и/или расположения контрольных объемов на самой мелкой сетке [1]. В У\,]Т такой зависимости нет. поэтому конфигурация контрольных объемов может быть задана после построения грубых сеток, но перед аппроксимацией модифицированной краевой задачи. Другими словами, искомая сеточная функция, которая является решением разностной краевой задачи, может быть определена как в точках х!, так и в х^, но па построение грубых сеток это ие влияет.

1.2.2. Построение грубых сеток. Построение грубых сеток в УуТ осуществляется посредством удаления двух точек из множеств (0;1) и С(0;1) самой мелкой сетки , как показано та рис. 1. Выделим первую точку х! из множества (0;1), которая будет первой точкой х! более грубой сетки С1. Далее пропустим две точки х! и х! из множества (0;1). Четвертая точка х| будет второй точкой х! более грубой с етки С1, седьмая точка х! - третьей точкой х! более грубой сетки С]; и т. д. То чка х| из множес тва С (0;1), расположенная посередине между точками х! и х| го множества С!(0;1), будет первой точкой х( более грубой с етки С];. Удаляя то две точки из мн ожества С (0;1), нетрудно получить остальные точки х( грубой сетки С1.

Далее аналогично построим еще одну грубую сетку С2,, но построение начнем с точек х! и хЗ. самой мелкой с етки , то есть со сдвигом па одну точку (рис. 1). Наконец, третья грубая сетка СЗ, строится тем же самым образом, но начиная с

Самая мелкая сетка С? Нулевой уровень

Рис. 2. Мпогосеточпая структура

точек Ж3 и Х4 самой мелкой сетки С°. Непосредственно из рис. 1 следуют основные свойства грубых сеток в УуТ:

1) грубые сетки С}, ^2 и С3 не имеют общих точек, то есть

С п ст = 0, п = т;

2) мелкая сетка С? представима в виде объединения грубых сеток С}, С} и

, то есть

с? = 0 С;

к=1

3) все сетки геометрически подобны, однако шаг грубых сеток С}, С^ и С^ в три раза больше, чем шаг сетки С°;

4) вне зависимости от способа определения сеточных функций на самой мелкой сетке каждый контрольный объем на сетках С}, С} и С3 является объединением трех контрольных объемов на сетке С}.

Самая мелкая сетка С} образует нулевой сеточный уровень, а три грубые сетки С}, С2 и С3 - первый сеточный уровень. Далее построение еще более грубых сеток осуществляется рекуррентным образом: каждая сетка С^, г = 1,..., 3Ь уровня Ь рассматривается как самая мелкая сетка для трех грубых сеток С^+}, у = = 1,3Ь+1 следующего уровня Ь + 1. Девять еще более грубых сеток, полученных из трех сеток первого уровня, образуют второй сеточный уровень и так далее, как показано на рис. 2. Построение грубых сеток завершается, когда на грубых сетках останется всего несколько точек хч и х(. В дальнейшем совокупность самой мелкой и всех грубых сеток будет называться миогосеточиой структурой.

Номер сеточного уровня с самыми грубыми сетками, обозначаемый как Ь+. вычисляется перед построением грубых сеток. Положим, что большинство самых грубых сеток имеют три точки х^ ш х(. Тогда число точек самой мелкой сетки С? теть N ИЛИ « 3Ь++1. Отсюда Ь+ определяется как

ЇІЗ-1

где квадратные скобки означают целую часть.

Сеточный уровень Ь состоит из сеток О^ к = 1,..., 3L с шагом Л.3Ь, где Н есть шаг самой мелкой сетки О}. Сетки О^ состоят из двух множеств точек О4(Ь;к) и О((Ь;к), которые могут быть как узлами, так и гранями контрольных объемов. Для удобства работы с многосеточной структурой обозначим отображение индексов точек данной сетки О^ та индексы точек самой мелкой сетки О} как Ж|4| и Ж{4}, где {г} и г есть индексы точек самой мелкой О0 и данной сетки о£ соответственно. Например, отображение точек из множества О¥(1;1) сетки О}, показанной на рис. 1, есть ж{}} = ж^, Ж|2} = ж|, Ж|3} = ж^, ... Отображение индексов позволяет формулировать разностные краевые задачи на многосеточной структуре практически так же, как и в односеточных алгоритмах.

Поскольку Ж-мериая сетка представима в виде произведения N одномерных сеток, то построение грубых сеток осуществляется независимо в каждом пространственном направлении. В этом случае сеточный уровень Ь состоит из 3NL сеток.

Перечислим основные свойства грубых сеток У\,]Т, которые будут использованы для численного решения краевых задач.

Свойство 1: каждый контрольный объем на сетках ОL представим в виде объединения 3L контрольных объемов на самой мелкой сетке О0. В сочетании со свойством аддитивности определенного интеграла относительно подобластей это позволит существенно расширить класс краевых задач, решаемых унифицированным образом.

Свойство 2:каждая сетка О^' (Ь = Ь+) представима в виде объединения трех соответствующих ей грубых сеток, что позволит исключить интерполяцию ИЗ УуТ.

О0

одного уровня:

3Ь

О0 =и О^ Ь = 0,...,Ь+

к=1

Свойство 3: сетки одного уровня не имеют общих точек, то есть

О^: п От = 0, п = т, Ь = 1,...,Ь+,

что позволит эффективно распараллеливать вычисления и экономно использовать память компьютера.

1.2.3. Аппроксимация модифицированных краевых задач на многосеточной структуре. Аппроксимация на многосеточной структуре пнтегро-ннтерполяцнонным методом мало отличается от «односеточной» аппроксимации [Ю, И].

Поскольку па границах области П = [0,1] заданы условия Дирихле, то точки сетки Ж|4} будут узлами, ж{} - гранями контрольных объемов, а сам объем определяется как

^ = {ж 1 Ж{ч-1} < ж < жм}.

Интегрирование уравнения (За) по контрольному объему V приводит к следующему соотношению

£ £

Х&} Ж{г>

ад(ж^о) — и|(ж/. 1г) 1 с 1 г

~ ~ПуГ, ч(х) с(ж) йж = —Т г(х)<Ь,

££

Х{*-1> Х{^-1>

где т = — Использование простейших интерполяций приводит к следующей разностной схеме

/ 7 \ С{*+1} ~ С{ > } /7,\ С{*} С{*—1} / \ _ / X

( ){*+1} /?2з2Ь \ }{*} /?2з2Ь ('3,){'} с{'} ){*}’

(«)и = У

т£

есть среднее значение функции д(ж) та контрольном объеме V и

/ v S -1

I x{i}

1 f dx

■ («){*}

1 dx

h3L J k(x)

■"{i-1} \ ■"{i-l}

Правая часть разностной E-модифицированной задачи

_.f {i}

dx

rf x{i-1}

является среднеарифметическим значением невязок, вычисленных в тех узлах самой мелкой сетки, которые расположены внутри данного контрольного объема.

Первое отличие аппроксимации модифицированных краевых задач на многосеточной структуре от односеточной аппроксимации состоит в вычислении коэффициентов (#){*}, (&){*} и правой части (г){4} . Алгоритм быстрого вычисления интегралов на многосеточной структуре изложен в [12]. Второе отличие состоит в аппроксимации граничных условий на грубых сетках, чьи граничные узлы не совпадают с границами области. Аппроксимация условий Дирихле и Неймана подробно изложена в [9. 12].

1

1.2.4. Многосеточные итерации. Многосеточная итерация У^]Т схематично показана на рис. 3. В У\,]Т используется только простейший (пилообразный) цикл (sawtooth cycle), который является частным случаем V-цикла [1]. Вычисления начинаются с уровня L+, состоящего из самых грубых сеток. Сначала на каждой сетке уровня L+ вычисляют коэффициенты (q}{j|, (k}{j| и правую часть (r}{j|, а затем выполняют сглаживающие итерации по мере достижения критерия останова. Поскольку сетки уровня L+ содержат лишь несколько узлов, то, как правило, требуется несколько сглаживающих итераций для вычисления поправки па данном уровне. Далее осуществляется переход к следующему уровню L+ — 1 с более мелкими сетками. Данный переход представляет собой изменение отображения индексов н не вносит никаких погрешностей в поправку, вычисленную на сетках уровня L+. На уровне L+ — 1 вычисления проводятся аналогичным образом: сначала определяются коэффициенты (q}{j|, (k}{j| и правая часть (r}{j}, а затем выполняются сглаживающие итерации. После выполнения сглаживающих итераций на самой мелкой сетке осуществляется пересчет приближения к решению (й := й + с) и обнуление поправки (с := 0). На этом многосеточная итерация заканчивается, проверяется критерий останова на самой мелкой сетке и в случае необходимости осуществляется следующая многосеточная итерация, начиная с уровня L+ (рис. 3).

Рис. 3. Пилообразный цикл и мпогосеточпая итерация

В многосеточных итерациях У^]Т отсутствует предварительное сглаживание, чтобы избежать дополнительных (по отношению к односеточным алгоритмам) проблем при решении нелинейных задач.

2. Вычислительные эксперименты

Результаты вычислительных экспериментов, связанных с решением различных двух- и трехмерных краевых задач, приводятся в [9]. В данной статье будет показана лишь возможность гибкого изменения способа п/пли порядка аппроксимации.

2.1. Изменения способа и/или порядка аппроксимации краевых задач. В качестве модельной рассмотрим задачу Дирихле для уравнения Пуассона

d2u d2u

+ dtf = U ап = ’ ( 5

где область Q является единичным квадратом.

Для наглядности рассмотрим наименее эффективный вариант У\,]Т . в котором множество используемых итерационных методов M состоит только из одного элемента, а именно метода Зейделя с точечным упорядочиванием неизвестных. Краевая задача (4) записывается в следующей Е-модифицированной форме:

d2c d2c

= г{х,у), с = -и, оа)

dx2 dy2 an

d2u d2u

Ф, у) = -f(x, у)- — -—. (5Ь)

Интегрирование модифицированной задачи (5) по контрольному объему приводит к следующей разностной схеме:

c{j-m ~ 2cwу + c{i+m , чи-1} ~ 2с{ія + c{ij+1} _ /,д (r>

^232 L + h232L v/{u}’ vDa/

f f

Ж{° 2 ~ 2 ~

«ш-jispr / / (ew

f f x{i-i} y{j-i}

В первом тесте (Тест 1а) правая часть Е-модифицированной задачи (6Ь) на самой мелкой сетке аппроксимируется на пятиточечном шаблоне:

/,,\ _ f U'i-ij ~ 2ui:j + iii+lj Uij-i — 2iiij + Wjj+i | гл/и2\

V Hj — Jij h2 /?2 + >■

Во втором тесте (Тест 16) правая часть Е-модифицированной задачи (6Ь) на самой мелкой сетке аппроксимируется на девятиточечном шаблоне:

h2 / d2/ d2A 1

= -fij + (^5^2 + Q^IJ ~ ^2 (“i-y+1 + 4“iJ + l + {l"i+lj+1 +

+ — 20Uij + + Ui-ij-i + 4Mjj_i + Ui+ij-1) + O(h4).

Предположим, что точное решение (ue) рассматриваемой краевой задачи (4) имеет вид

ue(x,y)= Q(x)Q(y), Q(e) = 20( ее + (1 — e)g — ^. (7)

Подстановка точного решения (7) в исходную задачу (4) определяет правую часть / (х,У)-

Скорость сходимости У^]Т и точность численного решения будем оценивать при помощи следующих параметров:

а) относительной евклидовой нормы невязки

1ИЙ№ — Ь||.

11Ы1 ’

б) погрешности численного решения

E(q) = max |«e(xv, yj) — й(|)|, ij

где q — номер многосеточной итерации.

Во всех вычислительных экспериментах итерационное решение разностных краевых задач будем начинать с нулевого начального приближения, то есть

U(0) =0 ^ R(0) = 1 и E(0) = max |ме(ж^, yj|. ij ij

Численное решение данной модельной краевой задачи осуществлялось на шестиуровневой многосеточной структуре (L+ = 5) с самой мелкой сеткой 501 х 501 (h = 1 /500). Сглаживающая итерация состояла из четырех итераций метода Зейде-ля с различными точечными упорядочиваниями неизвестных. Две сглаживающие итерации выполнены на каждом сеточном уровне. Результаты вычислительного эксперимента сведены в табл. 1.

Вычислительный эксперимент продемонстрировал важное преимущество У^]Т по сравнению с КуМ, которое возникает при Е-модификации краевой задачи: для получения высокоточного решения достаточно повысить порядок аппроксимации Е

Е

дачи использовать различные методы и различные сетки. Таким образом, У\,]Т можно обобщить для решения краевых задач на неструктурированных сетках.

Проиллюстрируем основные положения данного обобщения на примере решения задачи Дирихле в единичном квадрате. Построим триангуляцию для апЕ

Е

модифицированной краевой задачи (рис. 4).

Табл. 1

Сходимость технологій! в Тесте 1

Тест 1а Тест 16

Я К £ Я К £

0 1 17.95 0 1 17.95

1 8.653 • 1СГ3 1.753 • 1СГ1 1 8.656 • 10~3 1.753 • 10-1

2 5.690 • 1СГ6 1.126 • 1СГ3 2 5.696 • 10-6 1.122 • 10~3

3 4.571 • 10“' 1.414 • 1(ГБ 3 4.576 • 10“' 8.181 • 10~6

4 5.671 • 1СГ9 6.081 • 10~6 4 5.677 • 10-9 9.866 • 10~8

5 8.342 • 1СГ11 5.985 • 10~6 5 1.049 • Ю-10 1.388 • Ю-10

6 3.288 • 1СГ12 5.985 • 10~6 6 1.901 • 10~п 2.419 • 10~п

7 2.468 • 1СГ12 5.985 • 10~6 7 5.798 • 10~12 3.778 • 10~п

Рис. 5. Интерполяция невязки с триангуляции па ВСС (слева) и интерполяция поправки с ВСС па триангуляцию (справа)

Тогда многосеточная итерация У\,]Т представима в виде следующей последовательности действий:

1. Аппроксимация правой части Е-модифицированной краевой задачи при помощи МКЭ на исходной триангуляции.

2. Вычисление невязки.

3. Интерполяция невязки с исходной триангуляции на ВСС (рис. 5).

4. Аппроксимация левой части Е-модифицированной краевой задачи при помощи иитегро-иитерполяциоииого метода на ВСС.

5. Выполнение многосеточной итерации (рис. 3).

Табл. 2

Сходимость технологии в Тесте 1в

Я к £

0 1 17.95

1 1.536 • 1СГ2 2.045 • 10-1

2 8.566 • 1СГБ 1.956 • 10~3

3 6.123 • 1СГ' 1.841 • 10-6

4 7.252 • 1СГ9 8.401 • 10~6

5 9.663 • 1СГ11 6.005 • 10~6

6 4.446 • 1СГ12 5.985 • 10~6

7 3.474 • 1СГ12 5.985 • 10~6

6. Интерполяция поправки с ВСС на исходную триангуляцию (рис. 5).

7. Пересчет приближения к решению, обнуление поправки.

8. Проверка критерия останова, возврат к п. 1 (если необходимо).

Третий тест (Тест 1в) предназначен для иллюстрации сходимости У\]Т для МКЭ и по своим основным условиям совпадает с первым тестом (Тест 1а). Результаты данного теста, сведенные в табл. 2. показывают что использование двух сеток для аппроксимации левой и правой части £-модифицированной краевой задачи не оказало заметного влияния на скорость сходимости У\,]Т. В отличие от алгебраических многосеточных методов, рассмотренный вариант У\,]Т может быть

использован для решения нелинейных задач. По-видимому, при использовании

неструктурированных сеток невозможно исключить интерполяцию из миогосеточ-

ного алгоритма. Наличие такого проблемно-зависимого компонента, как интерполяция, неизбежно приведет к снижению уровня формализации при решении задач с разрывными коэффициентами. Отметим, что в УуТ интерполяция используется только на самой мелкой сетке.

Вычислительная область

Равномерная сетка в вычислительной области

Рис. 7. Управляющая функция Т : Рис. 8. Сетки в вычислительной и

Хг ^ Хг физической облаСТЯХ

Рассмотрим один из возможных вариантов построения ВСС. Предположим, что область вписана в единичный N-мерный куб. Пример двухмерной области показан па рис. 6. Разобьем область па треугольные элементы с вершинами ),

к = 1, 2,..., К. Для области, показанной па рис. 6, К = 409. Составим упорядоченное по возрастанию множество абсцисс вершин, из которого удалены совпадающие элементы

Тх = {ж* | 0 = Ж1 < Х2 < ... < ж* < ж*+1 < ... < жк = 1}, К ^ К .

В рассматриваемом примере К = 368. Далее построим па единичном отрезке равномерную сетку

Т' ! .Г ' 1 *=1,2,..., К}.

К — 1

Функция Т : X* ^ Хг, обычно называемая управляющей, отображает равномерную сетку в вычислительной области на упорядоченное множество абсцисс вершин (рис. 7). Однако в данном случае число узлов ВСС будет приблизительно равным К2. Потребуем, чтобы число узлов ВВС было приблизительно равным числу вершин исходной триангуляции. Для этого построим еще одну равномерную сетку

Т' ! Ж: ж: (' 'г *=1,2,..., К},

где

К

+ 1;

причем квадратные скобки означают целую часть. В рассматриваемом примере К = 21. Для отыскания положения узлов ВВС в физической области воспользуемся сплайн-интерполянтом управляющей функции Т : X* ^ ж*. Па рис. 8 показана полученная равномерная сетка в вычислительной области, соответствующая ей сетка в физическом области и силайн-интериолянт управляющей функции Т. Аналогичным образом строится сетка по каждому пространственному направлению. В данном примере число узлов ВСС составит К2 = 441, что несколько больше числа К = 409

Рис. 9. Исходная триангуляция и ВСС

Уравнение вида

в вычислительной области записывается как

А (,\' •''' + А (\vjv_ + А = _Я£1М1

дх V У'у% дх] ду V х'хК ду) дг V х'ху’у дг) Ку'уК ’

то есть с точностью до вида коэффициентов и правой части совпадает с исходным уравнением. Производные х'х, у'у и могут быть вычислены путем дифференцирования силайн-интерполяитов соответствующих управляющих функций. В [9] показано, что в данном случае сглаживающие процедуры на основе метода Зейделя могут оказаться неэффективными. Поэтому на уровнях с грубыми сетками следует применять более сильные сглажнвателн типа предобусловлеииых методов сопряженных градиентов. Пример построения комбинированной сглаживающей процедуры приведен в [9].

Summary

S.I. Martynenko. Formalization of Computations at Numerical Solution of Boundary Value Problems.

The paper represents a new robust multigrid technique for solving boundary value problems in black box manner. To overcome problem of robustness, the technique was based on incorporating the adaptation of boundary value problems to numerical methods, control volume discretisation and a new multigrid solver into a united computational algorithm. Special multiple coarse grid correction strategy makes it possible to obtain problem-independent, transfer operators. As a result, the most of the modes are approximated on coarse grids in order to ensure the problem of smoothing on the finest grid to be the least demanding. Detailed description of the robust multigrid technique and examples of application for solving benchmark problems are given in the paper.

Key words: boundary value problem, multigrid methods, finite element method, numerical experiment.

Литература

1. Wesseling P. An Introduction to Multigrid Methods. Chichester: John Wiley & Sons,

1991. 284 p.

2. Федоренко Р.П. Релаксационный метод решения разностных эллиптических уравнений // Журп. вычисл. матем. и матем. физ. 1961. Т. 1. .V 5. С. 922 927.

3. Федоренко Р.П, Скорость сходимости одного итерационного метода // Журп. вычисл. матем. и матем. физ. 1964. Т. 4. .V 3. С. 227 235.

4. Бахвалов Н. С. О сходимости одного релаксационного метода при естественных ограничениях па эллиптический оператор // Журп. вычисл. матем. и матем. физ. 1966. Т. 6, Л» 5. С. 861 883.

5. Шайдуров В.В. Мпогосеточпые методы конечных элементов. М.: Наука, 1989. 288 с.

6. Ольшанский М.А. Лекции и упражнения по мпогосеточпым методам. М.: ФИЗМА-

ЛИТ, 2005. 168 с.

7. Dendy Jr. J.E. Black box multigrid // J. Comput. Pliys. 1982. V. 48. P. 366 386.

8. Dendy Jr. J.E. Black box multigrid for systems // Appl. Math. Comp. 1986. V. 19.

P. 57 74.

9. Martynenko S.I. Robust Multigrid Technique for black box software // Comp. Met.li. in Appl. Math. 2006. V. 6, No 4. P. 413 435.

10. Калиткии П.П. Численные методы. М.: Наука, 1978. 512 с.

11. Самарский А.А. Теория разностных схем. М.: Наука, 1983. 616 с.

12. Мартыненко С.И. Универсальная мпогосеточпая технология для численного решения краевых задач па структурированных сетках // Вычислительные методы и программирование. 2000. Т. 1, Раздел 1. С. 85 104.

Поступила в редакцию 25.06.07

Мартыненко Сергей Иванович кандидат физико-математических паук, старший научный сотрудник отдела «Спеццвигатели и химмотология» ФГУП Центральный

Институт Авиационного Моторостроения им. П.И. Баранова, г. Москва.

E-mail: martyn_ sQmail.ru