CC BY
9
Вычислительные технологии
Том 23, № 3, 2018
Решение задачи Дирихле для уравнения Пуассона методом коллокации и наименьших квадратов в области с дискретно заданной границей
В. А. Беляев1'*, В. П. Шапеев1'2
1Институт теоретической и прикладной механики СО РАН, Новосибирск, Россия 2Новосибирский национальный исследовательский государственный университет, Россия *Контактный e-mail: [email protected]
Предложен и реализован новый вариант метода коллокации и наименьших квадратов повышенной точности для численного решения уравнения Пуассона. Реализованный алгоритм применяется в неканонических областях, границы которых заданы дискретно. Для приближенного и однозначного задания границы области по ее дискретным данным в прямоугольной системе координат строится параметрический двойной сплайн, в качестве компонент которого взяты два кубических сплайна. Используется идея присоединения вытянутых несамостоятельных нерегулярных граничных ячеек к соседним самостоятельным с целью уменьшения обусловленности глобальной системы линейных алгебраических уравнений.
Ключевые слова: метод коллокации и наименьших квадратов, краевая задача, уравнение Пуассона, повышенный порядок аппроксимации, дискретно заданная граница области, кубический сплайн.
Введение
Известно, что многие явления в природе происходят в областях с произвольной границей. Следовательно, необходимо иметь возможность моделировать эти явления в областях с формой, соответствующей действительности, чтобы на практике получать в расчетах моделируемых изделий, физических процессов и других объектов достоверные результаты. При этом зачастую для разных задач необходимо решать уравнение Пуассона в неканонических областях [1-10]. Однако многие численные методы, дающие высокую точность решения краевых задач для уравнений с частными производными (РВЕ) в прямоугольных областях, не имеют аналогов в случае областей иной формы.
Существуют численные методы, которые позволяют решать уравнение Пуассона в декартовых координатах [1, 2, 7-10]. К другой группе относятся численные методы, которые применяются для решения задач в ограниченном классе областей. Здесь уместно отметить, что задача решения с повышенной точностью РВЕ в области с криволинейной границей представляет собой трудную проблему для многих численных методов. Этот вопрос в литературе активно обсуждается для решения уравнения Пуассона в круге, кольце и цилиндре. Применение различных систем координат и сеток,
© ИВТ СО РАН, 2018
в частности полярной (цилиндрической) системы координат, не позволяет решить эту задачу с высокой точностью [3-6]. При переходе в полярную (цилиндрическую) систему координат в уравнении Пуассона перед производными появляются коэффициенты вида 1 /р и 1/р2, где р — радиальная компонента, поэтому при решении краевых задач для уравнения Пуассона при р ^ 0 возникает особенность, из-за которой снижается точность численного решения задачи. Некоторый обзор методов решения уравнения Пуассона, в частности в полярной системе координат, приведен в работе [5]. В данной работе, как и в [7-10], показано, что применение в методе коллокации и наименьших квадратов (КНК) регулярной сетки с прямоугольными ячейками в области с криволинейной границей и "одинарного" слоя нерегулярных ячеек (далее — н-ячеек) на ее границе позволяет эффективно решить эту проблему. Столь же эффективным для ее решения оказался несколько более сложный подход, изложенный в [11].
В методе КНК [5, 7-13] путем проектирования задачи для PDE в конечномерное линейное функциональное пространство ставится в соответствие приближенная задача, решение которой сводится к решению системы линейных алгебраических уравнений (СЛАУ). Решение последней определяет приближенное решение дифференциальной задачи. Поскольку в методе КНК приближенное решение кусочно-аналитическое (в данной работе кусочно-полиномиальное), в нем относительно просто реализуются варианты метода в областях различной формы [7-10] и на сетках с различной формой ячеек, в том числе варианты повышенной точности и на адаптивных сетках [11].
За последнее время метод КНК показал свою универсальность при решении краевых задач для PDE в неканонических областях [2, 5, 7-11]. В работах [7-10] продемонстрировано преимущество использования регулярной сетки с прямоугольными ячейками внутри области и одинарного слоя н-ячеек около границы области. В целях достижения аппроксимации повышенной точности дифференциальной задачи применялись "законтурные", не принадлежащие области решения дифференциальной задачи части граничных прямоугольных ячеек, отсеченных границей. Чтобы такой подход был применим, как и в случае других численных методов, предполагается, что решаемое уравнение справедливо не только в расчетной области, но и в малой окрестности ее границы. Кроме этого, "малые" по размеру (в сравнении с остальными) несамостоятельные н-ячейки, отсеченные границей области от прямоугольных ячеек начальной регулярной сетки, присоединялись к соседним самостоятельным ячейкам. В численных экспериментах установлено, что этот прием позволил существенно уменьшить обусловленность СЛАУ приближенной задачи и, как следствие, в сходящемся итерационном процессе сократить число итераций, необходимых для достижения заданной точности решения по сравнению со случаем, когда малые несамостоятельные н-ячейки наряду с другими ячейками использовались как самостоятельные для построения приближенного решения задачи. Важным преимуществом такого подхода является то, что он может быть применим к любой неканонической области с криволинейной границей, в отличие от других способов, которые в основном имеют более частный характер.
В данной работе уравнение Пуассона решается с повышенной точностью проекцион-но-сеточным методом КНК в областях, граница которых задана дискретным множеством точек, при этом некоторые из них могут быть точками ее излома. Излагаемый здесь подход может быть использован при решении краевых задач для других уравнений. В частном виде он применен в [14] при построении криволинейной образующей поверхности парового канала, отыскиваемой в процессе решения краевой задачи для трехмерного уравнения теплопроводности.
1. Постановка задачи и описание метода
Рассмотрим задачу Дирихле для уравнения Пуассона
Аи = /(хх,х2), (хх, х2) € О,
{
и
5о
= д(хх ,Х2)
(1)
в замкнутой области О, граница 6О которой задана дискретно конечной последовательностью точек, расположенных на ее границе (рис. 1). Здесь и(хх,х2) — искомая функция, /(х\,х2) и д(х\,х2) — заданные функции. Реализация рассматриваемого метода решения других краевых задач для уравнения Пуассона осуществляется аналогично алгоритму, описанному далее.
Пусть расстояние между заданными точками на границе порядка малой величины Н8. Если из этих данных или других соображений известно, что некоторые точки являются точками излома границы (например, кривая на рис. 3 — круг с "вырезом" имеет три точки излома), то разобьем ее на несколько участков так, чтобы концами этих участков были точки излома. Предположим, что между точками излома граница области имеет достаточную гладкость. Затем на каждом таком участке построим двойные сплайны (£), Х21 ^)). Здесь индекс г — номер участка в случае разбиения границы. Компоненты двойного сплайна (Ъ), Х2.— два кубических сплайна, являющиеся функциями параметра монотонно меняющегося при монотонном движении точки вдоль замкнутой границы. Пронумеруем все точки в возрастающей последовательности 0,..., Хр — 1 (рис. 1), где Хр — количество заданных точек. Во многих случаях в качестве параметра £ для наших целей удобно взять номер точки. Если известно, что входные данные соответствуют достаточно гладкой границе на всем ее протяжении,
' 0 Й ^ о > о о о 9 > : <А% Ф , , ❖ ❖ ❖ ❖ ^ о 14 ; □ □ : 1 15 16
<> о о ❖ ; ❖ 090 ❖ > , О О О ❖ : Т = "ГАИ] , ❖ ❖ ❖ о : □ <й <> о о : та: о о о Ъ : □ 12
' ❖ <> <> ❖ » ❖ ❖ ❖ ❖ > ; ❖ О5^ ❖ > , ❖ <> <> ❖ ' ❖ ❖ ❖ о , ❖ ❖ ❖ ❖ <> <> <> ❖ : о ^ л о 7 ^ ; ❖ О О О • □ □ - П° 8
<> ❖ ❖ ❖ : о о о ❖ > ; ❖ 010 ❖ ) , ❖ <> <> о ' ❖ ❖ ❖ о : ^ ° } ; ❖ 020 О > , ❖ ❖ ❖ О ООО О : о о ❖ ❖ > ; ❖ 030 ❖ > ❖ ❖ <□ ❖ : ^ л о!ъ 4 о • ; О О О Р о о о <в
Рис. 1. Область решения задачи: • — точки, дискретно задающие границу области (как и на рис. 2); Ь — номер точки. Номера точек возрастают в направлении против хода часовой стрелки
Рис. 2. Фрагмент расчетной области: о — точки записи уравнений коллокации; х — точки записи условий согласования; □ — точки записи краевых условий; о — центры граничных н-ячеек. Несамостоятельные н-ячейки 8 и 14 присоединяются к ячейкам 7, 10 соответственно
то для аппроксимации границы может быть достаточно построения единого сплайна с использованием всех заданных точек. В этом случае точки начала и конца сплайна совпадают. С целью получения "краевых условий" для построения сплайнов можно на концах каждого участка по координатам трех (или четырех) заданных точек построить интерполяционные полиномы Лагранжа (или Ньютона) для каждой компоненты двойного сплайна. Теперь для однозначного построения компонент двойных сплайнов достаточно положить, что их вторые производные на концах участков границы равны вторым производным от соответствующих интерполяционных полиномов. Известно, что в случае дискретных данных, соответствующих достаточно гладким кривым, построенные таким способом сплайны будут иметь четвертый порядок аппроксимации и сходимости при hs ^ 0 [15]. Если применяемый далее метод КНК решения краевой задачи для PDE в прямоугольной области будет четвертого порядка сходимости и имеет достаточную точность аппроксимации краевого условия задачи (1), то можно ожидать, что рассматриваемый здесь способ аппроксимации границы позволит получить в области с гладкой криволинейной границей также четвертый порядок сходимости. Чтобы достичь более высокого порядка сходимости, необходимо повысить степень базисных полиномов используемого пространства и повысить точность аппроксимации границы. Не представляет собой существенной трудности повысить на два порядка точность аппроксимации рассматриваемой здесь задачи.
Впишем расчетную область Q в прямоугольник, две стороны которого параллельны оси Х\, равны d\ и проходят через принадлежащие 5Q точки с минимальной и максимальной координатами по оси х2, а две другие стороны, параллельные оси х2, равны d2 и проходят через принадлежащие 5Q точки с минимальной и максимальной координатами по оси Х\. Полученную прямоугольную область размера d\ х d2 покроем регулярной сеткой с прямоугольными ячейками Qj размера 2h\ х 2h2, j = 1,..., N, N = N\ х N2 — количество ячеек. Как отмечено в [7], в случае наличия на границе области точек ее излома желательно изначально построить сетку так, чтобы они оказались на координатных линиях сетки. При решении задачи на последовательности сеток с кратным измельчением ячеек точки излома границы останутся на координатных линиях сетки. Иначе изначально и при дальнейшем измельчении ячеек будут присутствовать маленькие ячейки сетки. Из дальнейшего изложения будет ясно, что наличие таких ячеек зачастую уменьшает точность приближенного решения задачи, а также усложняет алгоритм построения сетки.
Прямоугольные ячейки, которые полностью расположены внутри области, назовем внутренними (рис. 2, ячейки 1-3, 5-7, 9, 10). Ячейки, пересеченные границей области, для краткости здесь назовем граничными (ячейки 4, 8, 11, 13, 14). Часть граничной ячейки, отсеченную границей и лежащую внутри области, назовем н-ячейкой. Н-ячейка в случае гладкой границы и достаточно мелких ячеек сетки может быть либо треугольного вида (далее — треугольная ячейка (ячейки 8, 11, 14)), либо четырехугольного (далее — четырехугольная ячейка (рис. 2, ячейки 4, 13)). При этом часть граничной ячейки, лежащую вне области, назовем законтурной. Прямоугольную ячейку сетки, от которой границей области отсечена н-ячейка, назовем материнской. Назовем условно внешней стороной н-ячейки часть границы 5Q области Q, которая оказалась внутри граничной ячейки. Стороны и отсеченные границей части сторон граничной ячейки, расположенные внутри области, назовем внутренними сторонами н-ячейки. Решение в прямоугольных ячейках, которые не содержат какую-либо часть области Q, естественно, не строим (ячейки 12, 15, 16).
Для удобства реализации метода в каждой ячейке области введем локальные координаты
Х\ — X\j Х2 — X2j
У1 = —z—-, У2 = —т—-, (2)
hi h2
где (xij, X2j) — центр j-й ячейки, hi = di/(2Ni), h2 = d2/(2N2), v(yi, У2) = u(xi,x2). Координаты центра ячейки в данном случае определяются как среднеарифметические величины от соответствующих координат вершин ее углов. Задача (1) после замены (2) в локальных переменных примет вид
1 d2v 1 d2v
h2w + h2 w = f (Xi(yi), Х2(У2)), (yi, У2) e Q c Q, (3)
V = g(Xl(Уl), X2(У2)), (У1,У2) e Ш П 6Qj, j = 1,...,N.
В каждой j-й ячейке сетки приближенное решение Uhj задачи (3) ищем в виде линейной комбинации с неопределенными коэффициентами базисных элементов пространства полиномов четвертой степени от двух переменных. Здесь в качестве базисных элементов фг (г = 1,..., 15) взяты мономы
15
Vhj (yi, У2) = ^2 brjфг = bij + b2jyi + b3jУ2 + b4jyi2 + b5jУ22 +
r=i
+bßj yrn + b7j yi2y2 + b8j yrn2 + bgj yi2y22 + bWj yi3 + bnj у23 + b^ yi3y2+
+bi3j ym3 + bUj yi4 + bi5j У24.
Глобальное приближенное решение будет состоять из "локальных" аналитических "кусков", связанных с центром соответствующей ячейки. Для определения неизвестных коэффициентов brj в каждой ячейке выписывается переопределенная "локальная" (частичная) система линейных алгебраических уравнений (СЛАУ). Она состоит из переопределенной системы уравнений, имеющей в каждой ячейке уравнения коллока-ции, условия согласования на общих сторонах, принадлежащие двум соседним ячейкам, и краевые условия на 5Q, если ячейка является граничной. Объединение всех частичных СЛАУ определяет глобальное решение. При реализации представленного в этой работе варианта метода КНК число точек записи уравнений коллокации, условий согласования и краевых условий в каждой ячейке в сумме было равно 32.
Если в расчетной сетке имеется или появилась после измельчения шагов сетки маленькая и/или вытянутая н-ячейка (рис. 2, ячейки 8, 14), то глобальная СЛАУ задачи, как правило, становится хуже обусловленной. Такие ячейки могут быть причиной понижения точности приближенного решения задачи на содержащей их сетке. Чтобы избежать этого, здесь используется идея присоединения таких н-ячеек к соседним ячейкам.
Назовем несамостоятельными н-ячейки, расположенные в граничных ячейках, в которых начало локальной системы координат находится вне расчетной области. Предлагается несамостоятельную н-ячейку присоединять к той соседней самостоятельной ячейке, с которой она имеет наибольшей длины сторону среди ее сторон, общих с другими соседними самостоятельными ячейками. Если несамостоятельная н-ячейка имеет одинаковую длину с двумя общими сторонами соседних с ней ячеек, то ее присоединяем к ячейке с наименьшим индексом j. Через эту сторону в присоединенную н-ячейку(ки) продолжается решение из ячейки, к которой ее (их) присоединили. При этом согласно
сформулированному правилу к ячейке могут быть присоединены несколько малых несамостоятельных н-ячеек. В этом случае внешняя сторона объединенной ячейки состоит из внешних сторон всех объединенных в ней ячеек. На этой внешней стороне расставляется столько точек для записи краевых условий, чтобы в сумме число точек записи уравнений коллокации, условий согласования и краевых условий в объединенной ячейке было равно 32. Эти краевые условия включаются в переопределенную систему уравнений объединенной ячейки. Для простоты реализации варианта полагаем, что начало локальной системы координат в объединенной ячейке совпадает с началом координат в исходной ячейке, к которой присоединяли ячейку(ки). Объединенную ячейку также считаем н-ячейкой.
В реализованном здесь варианте метода КНК несамостоятельные н-ячейки, имеющие в качестве соседних ячеек только несамостоятельные н-ячейки, не присоединяются к соседним несамостоятельным н-ячейкам и считаются условно самостоятельными. В них также по изложенному в данной работе алгоритму выписывается локальная СЛАУ из 32 уравнений и строится решение, согласованное с локальными решениями в соседних ячейках. Кроме того, при написании компьютерной программы устанавливалось следующее ограничение на присоединение ячеек. К самостоятельной ячейке, которая имеет Nq (q = 1,..., 4) прямолинейных сторон, максимум можно присоединить Nq — 1 несамостоятельных н-ячеек. Случай, когда возникает необходимость присоединить к самостоятельной ячейке (имеющей Nq прямолинейных сторон) Nq несамостоятельных н-ячеек, маловероятен на практике. При возникновении такой исключительной ситуации компьютерная программа выдает соответствующее сообщение пользователю. В проделанных авторами различных экспериментах такой случай не возникал.
Здесь затруднительно указать ограничения на вид неканонических областей и весь класс решаемых в них задач, для которых будет успешно применим предложенный здесь подход. Однако очевидно, что этот класс достаточно широк. Естественно, что если реализовать какой-то более сложный вариант метода КНК, то класс областей, для которых он будет применим, может расшириться.
Определим по следующему алгоритму расстановку точек коллокации, записи условий согласования в любой ячейке или краевых условий для граничной ячейки. В каждой ячейке для точек коллокации (рис. 2) возьмем следующие локальные координаты:
Если у ячейки есть соседняя ячейка, то их общую сторону (в случае граничных ячеек — общую сторону материнской ячейки) поделим на четыре равных отрезка. В середине каждого отрезка запишем условия согласования. Очевидно, что в некоторых граничных н-ячейках точки для записи уравнений коллокации или условий согласования могут быть вне области П (рис. 2, ячейки 4, 13 и 11). Чтобы такой подход был применим, как и в случае других численных методов, предполагается, что решаемое уравнение справедливо не только в расчетной области, но и в малой окрестности ее границы.
Пусть на внешней криволинейной стороне н-ячейки требуется расположить в точек (р\,р2, ...,Рз) для записи краевых условий. Зная вершины прямоугольной материнской ячейки в прямоугольной системе координат (х\, х2), находим соответствующие значения начала и конца внешней стороны ячейки Ь*, ¿** (¿** > ¿*) пересечением
сторон прямоугольной материнской ячейки (а в случае присоединения н-ячеек пересечением продолжения сторон материнской прямоугольной ячейки) с соответствующим участком сплайна. Эта задача аналогична задаче поиска значения одного действительного корня £ уравнения третьей степени на определенном отрезке. Для ее решения в данном варианте построения расчетной сетки использовался метод биссекции (метод деления отрезка пополам). Затем каждой точке р^ ставим в соответствие число
(2] — 1)(г** — г)
а^ = ¿* +--, где ] = 1,..., е. Для каждого числа aj находим, какому куску
2 $
двойного сплайна (I), Х21 (I)) соответствующая точкар^ принадлежит, где г —номер участка границы. Пусть она принадлежит к-му куску двойного сплайна: к — 1 < а^ < к (к = 0,..., Мр). Тогда точка р^ для записи краевых условий будет иметь следующие координаты в прямоугольной системе координат (х\, х2): (Х1к(а^), Х2к(а^)). Такой способ расстановки точек для записи краевых условий в случае прямолинейной внешней стороны н-ячейки эквивалентен делению этой стороны на в равных отрезков для записи в серединах этих отрезков краевых условий, как это было описано выше. Если самостоятельная н-ячейка имеет три, две или одну соседнюю ячейку, то на внешней ее стороне расставляем четыре, восемь или двенадцать точек соответственно для записи краевых условий по правилу, описанному выше.
Можно отметить, что существует много различных способов расстановки точек для записи уравнений коллокации, условий согласования и краевых условий. При этом во многих случаях полученные приближенные решения будут с высокой точностью совпадать друг с другом. Расстановка точек для записи перечисленных уравнений, в некоторой мере близкая к симметричной и равномерной, в методе КНК дает хорошие результаты.
На рис. 2 изображена расчетная область, покрытая регулярной сеткой с прямоугольными ячейками, в которых по описанному выше правилу расставлены точки записи уравнений коллокации, условий согласования и краевых условий.
Уравнения коллокации в каждой ]-й ячейке выписываются в 16 точках коллокации и имеют вид
1 9(у1^ У2с) 1 д2Ыгэ У2с) ( , , ^
ц—^— + ц—^— =1 Ыу1с),Х2(Шс)).
Здесь (у1с, у2с) — точки коллокации, где с = 1,..., 16.
В качестве условия согласования решения в четырех указанных выше точках согласования на каждой общей стороне между соседними ячейками требуем непрерывности линейной комбинации с весами функции V и ее производной по нормали:
к1 + к1 дь-
ц + = + ^.
Здесь п^ — внешняя нормаль к границе ]-й ячейки; = к1, если направление п^ совпадает с направлением оси у1, и = к2, если направление п^ совпадает с направлением оси у2; и V- — пределы функции при стремлении изнутри и извне к границе ]-й ячейки соответственно; к1, к2 — весовые параметры.
Если сторона ]-й ячейки совпадает с некоторым куском границы области, в указанных выше точках на ней выпишем краевые условия
= к3д(х1, Х2),
где к3 — весовой параметр краевого условия в переопределенной системе уравнений.
В каждой ячейке, объединяя уравнения коллокации, условия согласования и краевые условия (в случае граничных ячеек), относительно неизвестных коэффициентов brj получим переопределенную локальную СЛАУ вида
15
У^ В[г ь,.
Fh I = !,..., 32, j = 1,...,N,
r=1
определяющую приближенно локальное решение в окрестности начала системы координат в ячейке. Для того чтобы определить, что понимается под решением этой системы, рассмотрим функционал
ф = V ( > л в1гЪГт- Ц) . (5)
32/15 \ 2
£ £ Blr brj - fA .
1=1 v=i J
Он представляет собой сумму квадратов невязок всех уравнений в ячейке на приближенном решении задачи. Решение системы (4) на каждой итерации находится из условия минимума функционала Ф по коэффициентам brj, г = 1,..., 15 [13]. Система уравнений, полученная объединением уравнений во всех ячейках расчетной области (глобальная СЛАУ), решается в итерационном процессе Гаусса — Зейделя. В ней одна "глобальная итерация" состоит из последовательного решения локальных СЛАУ (4) во всех ячейках области [8].
Итерационный процесс продолжается до тех пор, пока не выполнится условие
тах 1ЬГуП+1 - Ьг™1 < е,
•Г]
где ЬГ^ — г-й (г = 1,..., 15) коэффициент полинома, аппроксимирующего решение в ячейке с номером ] на п-й итерации. Величина е — заданная константа, называемая псевдопогрешностью решения. Ее значение выбирается таким, чтобы погрешность решения глобальной СЛАУ была существенно меньше погрешности аппроксимации исходной задачи.
2. Численные эксперименты и результаты расчетов
В численных экспериментах на сходимость приближенного решения на последовательности сеток при измельчении шагов сетки вдвое рассматривалась задача (1) в областях, граница которых задана дискретно. В качестве тестовых решений брались аналитические функции, подстановкой которых в уравнения подбиралась их правая часть для того, чтобы эти функции удовлетворяли рассматриваемым дифференциальным уравнениям. Краевые условия брались из значений тестовых функций и их производных на 80. Во всех представленных ниже таблицах приведены значения погрешности приближенного решения в равномерной норме
II u - uh ||с= max max Iu(x1m, Х2т) - uhj(х1т, X2m)l,
j=1,...N m=1,...Lj
где N — количество расчетных ячеек в области; Lj — количество равномерно распределенных контрольных точек (х1т, х2т), взятых в j-й ячейке для подсчета в них погрешности; и — точное решение задачи (1); Uh — приближенное решение в области 0;
и^ — приближенное решение в ]-й ячейке. Порядок сходимости численного решения определяется следующим образом:
-, (6)
^N1 ,N2
где Е^1,и2 — значение ее абсолютной погрешности на сетке размера N = х Ы2; Еи1/2,и2/2 — значение ее абсолютной погрешности на сетке Ы1/2 х N2/2. На самом деле порядок сходимости численного решения вычислялся по формуле (6) формально, так как в случае нерегулярной сетки, строго говоря, определение порядка сходимости отсутствует. Тем не менее величина (6) показывает порядок малости погрешности решения относительно малой величины шага сетки. В численных экспериментах были взяты следующие значения параметров: Ь^ = 100, е = 10-14, к1 = 1, к2 = 5(| сов^)^ +1 вт(а)|Л,2), ((сов(а), вт(а)) — координаты внешней нормали к стороне ячейки (а = 0, п/2, ж, 3ж/2), к3 = 1.5. Во всех численных экспериментах здесь в начальном приближении решения взяты brj = 0.4, где г = 1,..., 15, ] = 1,..., N. Обозначим
через количество итераций. В расчетах использовано комбинированное применение операции продолжения на многосеточном комплексе в методе Федоренко и на промежуточных сетках комплекса метода ускорения сходимости итерационного процесса, основанного на подпространствах Крылова [16-19]. Такой прием ускорения итерационного процесса дал значительное сокращение времени решения задачи по сравнению со случаем без его применения.
Пример 1. Рассмотрим единичную окружность (рис. 2), граница которой задана 40 точками. В данном случае строим один двойной сплайн. Проведем также численные эксперименты и для случаев, когда задано 20, 30 и 60 точек границы. В табл. 1 представлены результаты численных экспериментов в этих случаях.
Из приведенных результатов видно, что для достижения лучшей точности необходимо во входных данных задачи задавать достаточное количество точек на границе расчетной области, т. е. максимально точно аппроксимировать границу области. Необходимое количество точек должно быть таким, чтобы погрешность интерполяции гра-
Таблица 1. Результаты численного эксперимента с решением и(ж1, х2) = 81п(ж1) + 8т(ж2) тестовой задачи в единичной окружности
мр х Ы2 Еи1,и2 NПег , ^N1/2^2/2
20 20х20 40х40 80х80 8.6 ■ 10-7 Расходимость Расходимость 92 —
30 20х20 40х40 80х80 1.0 ■ 10-9 4.6 ■ 10-11 5.4 ■ 10-12 97 204 472 4.4 3.1
40 20х20 40х40 80х80 1.0 ■ 10-9 3.2 ■ 10-11 1.6 ■ 10-12 96 204 452 5.0 4.3
60 20х20 40х40 80х80 1.0 ■ 10-9 3.6 ■ 10-11 1.1 ■ 10-12 96 202 463 4.9 5.0
Х2 1 +
0
1 x1
Х2
1 А Л А А
0
1 x1
Рис. 3. Области, границы которых заданы 50 точками (три точки излома) (а) и 200 точками (пять точек излома) (б): х — точки излома границы
а
Таблица 2. Результаты численного эксперимента с решением и(х 1,ж2) = 8т(ж1) + вт(ж2) тестовой задачи в области, граница которой на рис. 3, а задана 50 точками
Nx х N2 En^nI Niter , En1/2,n2/2 log2 ' ^n1,n2
10 х 10 1.7 ■ 10-8 81 —
20 х 20 1.0 ■ 10-9 96 4.1
40 х 40 3.6 ■ 10-11 205 4.8
80 х 80 1.9 ■ 10-12 452 4.2
ницы сплайнами не превышала погрешность аппроксимации дифференциальной задачи методом КНК.
Пример 2. Рассмотрим область (рис. 3, а), граница которой задана 50 точками. В данном случае строим три двойных сплайна, так как у нас три точки излома границы. Здесь точками границы являются точки, расположенные на трех окружностях: (хх — 0.5)2 + (х2 — 0.5)2 = 0.52, (хх — 0.15)2 + (х2 — 0.5)2 = 0.52, (xi — 0.75)2 + (х2 — 0.5)2 = 0.552. В табл. 2 приведены результаты численного эксперимента.
Пример 3. Рассмотрим область (рис. 3, б), граница которой задана 200 точками. В данном случае строим пять двойных сплайнов, так как у нас пять точек излома границы. Здесь граница области составлена из полуокружности (хх — 0.5)2 + (х2 — 0.5)2 = 0.52, х2 £ [0, 0.5] и следующих аналитических кривых, задаваемых уравнениями: ж2 = 0.5 + sin(4rai)/2, хх £ [0,0.25); ж2 = 0.5 — sin(4^i)/2, хх £ [0.25,0.5); ж2 = 0.5 + sin(4ra1)/2, хх £ [0.5, 0.75); ж2 = 0.5 — sin(4^xx)/2, хх £ [0.75,1].
Сначала в этой сложной области с изломами границы рассмотрим тестовую задачу с решением без особенностей. Результаты численного решения на последовательности сеток приведены в табл. 3. Близкие по точности решения результаты ранее были получены при решении этой задачи в четырехугольных областях [8-10]. Это говорит о достоинствах рассматриваемого здесь подхода. Некоторая потеря точности решения в области с криволинейной границей по сравнению со случаем прямоугольной области прежде всего связана с наличием н-ячеек. Стоит отметить, что здесь на грубой сетке возникают н-ячейки, в которых граница области разрезает их насквозь внутри
Таблица 3. Результаты численного эксперимента с решением и(х1,х2) = 8ш(ж1) + 8т(ж2) тестовой задачи в области, граница которой на рис. 3, б задана 200 точек
N1 х N2 En1,n2 Niter , ENi/2,N2/2 log2 ' ^N1, N2
40 х 40 3.6 ■ 10-11 266 —
80 х 80 1.9 ■ 10-12 612 4.2
Таблица 4. Результаты численного эксперимента в области, граница которой на рис. 3, б задана 200 точек с тестовым решением м(ж1,ж2) = е10х1 + е10ж2 + е5х1+5х2, имеющим большие градиенты. Здесь е = 0.5 ■ 10-10
N1 х N2 EN1,N2 Niter , EN1/2N2/2 log2 ' ^N1, N2
40 х 40 4.9 ■ 10-2 260 —
80 х 80 1.6 ■ 10-3 620 4.9
на несколько частей. Поэтому приходится начинать счет с более мелкой сетки размера 40 х 40, чтобы был применим предложенный здесь вариант метода КНК.
Теперь рассмотрим в этой области тестовую задачу с большими градиентами искомой функции в краевом условии и соответственно решения задачи. В табл. 4 приведены результаты численного эксперимента. Из нее видно, что решение этой задачи сходится с повышенным порядком, что обеспечено повышенной степенью базисных полиномов в рассматриваемом функциональном пространстве. Но в коэффициент остаточного члена аппроксимации задачи входят большие значения производных решения. Поэтому на рассматриваемой сетке погрешность решения существенно больше, чем в случае решения тестовой задачи без особенностей.
В численных экспериментах в некоторых других областях с криволинейными границами были получены аналогичные результаты.
Пример 4. Рассмотрим пример задачи Дирихле для уравнения Пуассона (1) с f (жьж2) = 1, g(x1, х2)\sn = 0 [20, 21] в правильном шестиугольнике (рис. 4). Пусть граница области задана 60 точками, расположенными на ней равномерно. В данном примере построим шесть двойных сплайнов, так как в этом случае имеется шесть точек излома границы.
Для оценки порядка сходимости численного решения в равномерной норме на последовательности сеток вычислялась величина
Е =| UNl,N2 - UNl/2,N2/2 I с= maxJ«Nb N2 (ХЦ, %2£,) - «^1/2,^2/2(^1^, Х2£ )\, (7)
где M — количество равномерно распределенных контрольных точек (ж^, х2^), взятых для простоты подсчета в прямоугольнике, включающем рассматриваемую область (в точках, лежащих вне шестиугольника, погрешность полагалась равной нулю); uNl,N2 — численное решение, полученное на сетке N1 х N2; uNl/2,N2/2 — численное решение на сетке ^1/2 х N2/2. Порядок сходимости решения формально определялся по формуле
Ер _ , ^ UNi,N2 — Uni/2,n2/2 ¡с
log2 -г = log
2
Ес 1 UN1/2N2/2 — UNi/4N2/4 ¡с
Рис. 4. Граница расчетной области в примере 4 с шестью точками излома (а), на (б) — Я — радиус круга, на границе которого расположены вершины шестиугольника, Яс — радиус окрестности с центром шестиугольника, Я8 — радиус окрестности с центром в точке излома
где Ер =|| и^,N2 — и^/2^2/2 Ус, Ес =|| и^/2^2/2 — ^¡4,^2/4 ||с и и^/4 — численное решение на сетке N1 /4 х И2/4 [8, 10]. В табл. 5 приведены полученные в численном эксперименте значения нормы (7) на последовательности сеток и соответствующие порядки ее сходимости при М = 250 000.
Замечание 1. В данном примере установлено, что при значительном увеличении точек для описания границы области значения норм (7) практически не изменяются.
Таблица 5. Результат численного эксперимента для примера 4
У ,N2 — ^N1 /2^2/2 Ус Е\яа=0 10^2 | Е\да = К/10 10^2 | Е\Да = Я/6 10®2 |
У ^40,40 — ^20,20 Ус 2.69 • 10-4 — 2.62 • 10-4 — 1.70 • 10-4 —
У ^80,80 — ^40,40 Не 9.21 • 10-5 1.54 4.75 • 10-5 2.46 2.87 • 10-5 2.56
У Щ60,160 — Щ0,80 У|с 2.78 • 10-5 1.72 6.39 • 10-6 2.89 3.17 • 10-6 3.17
У Щ20,320 — Щ60,160 У|с 1.36 • 10-5 1.36 4.02 • 10-7 3.99 1.05 • 10-7 4.91
У ,N2 — ^N1 /2^2/2 У|с Е\ис=К/2 ^ | Е \дс=я/6 ^ | Е\дс=к/10
У ^40,40 — ^20,20 У|с 1.04 • 10-5 — 9.48 • 10-6 — 9.45 • 10-6 —
У ^80,80 — ^40,40 Ус 1.63 • 10-6 2.67 1.40 • 10-6 2.75 1.39 • 10-6 2.76
У и1б0,160 — ^80,80 У|с 2.21 • 10-7 2.88 1.83 • 10-7 2.93 1.81 • 10-7 2.94
У и320,320 — ^160,160 У|с 2.83 • 10-8 2.96 2.31 • 10-8 2.98 2.27 • 10-8 2.99
Примечание: запись Е|даозначает величину нормы (7) в шестиугольнике, за исключением окрестностей радиуса Я8 = Л/10 точек излома его границы, Е=я/ю означает величину нормы (7) только в окрестности радиуса Яс = Л/10 центра шестиугольника (рис. 4).
Здесь 1о§2 ^ — порядок сходимости нормы (7), значение которой приведено в предыдущем
Ес
столбце.
Замечание 2. Результаты, полученные на последовательности сеток и приведенные в табл. 5, позволяют оценить значение погрешности приближенного решения.
Очевидно, в этом примере, как и в работах [20, 21], из-за наличия локальных особенностей в краевом условии в точках излома границы возникает особенность в решении задачи. Это приводит к увеличению погрешности численного решения задачи и снижению порядка его сходимости, значение которого существенно зависит от рассматриваемой подобласти, взятой в шестиугольнике. Из результатов, приведенных в табл. 5, видно следующее:
1) максимум погрешности решения находится в окрестности точек излома границы области (i?\_Rs=o);
2) при учете значений решения во всех контрольных точках наблюдается порядок его сходимости выше первого;
3) в подобластях, не содержащих окрестности вершин шестиугольника c Rs > 0, наблюдается повышенный порядок сходимости решения, при этом величина нормы (7) несколько меньше, чем в случае примеров с Rc > 0.
Аналогичный результат был получен в случае решения краевых задач для бигармо-нического уравнения с априори неизвестным решением, когда правая часть уравнения задавалась "независимо" от краевых условий [8, 10]. Там также наблюдался повышенный порядок сходимости численных решений, полученных с помощью аналогичных вариантов метода КНК.
Для того чтобы добиться более высокого порядка сходимости решения, можно увеличить порядок аппроксимации приближенной задачи (степень полиномов в базисе используемого пространства) или/и исключить главную часть значения асимптотики решения в окрестности вершин шестиугольника. Как известно, указанные приемы с этой целью используются при решении задач с аналогичными особенностями и другими методами, как это, например, сделано в случае прямоугольной области в работах [12, 22].
Можно также отметить, что эти результаты коррелируют с результатами работы [21], полученными при исследовании конечно-разностным методом поведения решения вышеупомянутой тестовой задачи. Эту задачу можно рассматривать в качестве эталонной при исследовании возможностей приближенных методов для решения краевых задач с особенностями для PDE.
Примененный здесь способ описания сплайнами дискретно заданной границы расчетной области можно использовать при численном решении краевых задач и другими методами. В работе [14] этот способ использовался при решении задачи конечно-разностным методом.
Заключение
Предложен и реализован новый вариант метода КНК с повышенным порядком точности для численного решения задачи Дирихле для уравнения Пуассона в областях, граница которых задана дискретно. Дифференциальная задача проектировалась в пространство с базисом из полиномов четвертой степени. Метод верифицирован на решении пяти тестовых задач. Из результатов, приведенных в табл. 1-5, видно, что численное решение имеет повышенный порядок сходимости. В случаях, когда в качестве правой части бралась функция, полученная подстановкой тестовой достаточно гладкой
функции от двух переменных, значения которой на границе расчетной области брались в качестве краевого условия, наблюдался четвертый порядок сходимости. При наличии особенности типа разрыва производных в точках излома границы области в численном решении задачи наблюдалось понижение порядка сходимости по отношению к предыдущему случаю.
Список литературы / References
[1] Johansen, H., Colella, P. A cartesian grid embedded boundary method for Poisson's equation on irregular domains // Comput. Phys. 1998. Vol. 147(1). P. 60-85.
[2] Ш^апеев А.В., ШШапеев В.П. Разностные схемы повышенной точности для решения эллиптических уравнений в области с криволинейной границей // Журн. вычисл. математики и матем. физики. 2000. Т. 40, №. 2. С. 223-232.
Shapeev, A.V., Shapeev, V.P. High-order accurate difference schemes for elliptic equations in a domain with a curvilinear boundary // Comput. Math. and Math. Phys. 2000. Vol. 40, No. 2. P. 213-221.
[3] Lai, M.C., Lin, W.W., Wang, W. A fast spectral/difference method without pole conditions for Poisson-type equations in cylindrical and spherical geometries // IMA J. Numer. Anal. 2002. Vol. 22, No. 4. P. 537-548.
[4] Lai, M.C. A simple compact fourth-order Poisson solver on polar geometry //J. Comput. Phys. 2002. Vol. 182, iss. 4. P. 337-345.
[5] Ворожцов Е.В., ШШапеев В.П. Численное решение уравнения Пуассона в полярных координатах методом коллокаций и наименьших невязок // Моделирование и анализ ин-форм. систем. 2015. Т. 22, № 5. С. 648-664.
Vorozhtsov, E.V., Shapeev, V.P. Numerical solution of the Poisson equation in polar coordinates by the method of collocations and least residuals // Modelling and Analysis of Inform. Systems. 2015. Vol. 22, No. 5. P. 648-664. (In Russ.)
[6] Shapeev, V.P., Vorozhtsov, E.V. Application of the method of collocations and least residuals to the solution of the Poisson equation in polar coordinates // J. of Multidisciplinary Eng. Sci. and Technology. 2015. Vol. 2, No. 9. P. 2553-2562.
[7] ШШапеев В.П., Беляев В^. Варианты метода коллокации и наименьших невязок повышенной точности в области с криволинейной границей // Вычисл. технологии. 2016. Т. 21, № 5. С. 95-110.
Shapeev, V.P., Belyaev, V.A. Versions of high order accuracy collocation and least residuals method in the domain with a curvilinear boundary // Comput. Technologies. 2016. Vol. 21, No. 5. P. 95-110. (In Russ.)
[8] Беляев В.А., ШШапеев В.П. Варианты метода коллокации и наименьших невязок для решения задач математической физики в трапециевидных областях // Вычисл. технологии. 2017. Т. 22, № 4. С. 22-42.
Belyaev, V.A., Shapeev, V.P. The versions of collocation and least residuals method for solving problems of mathematical physics in the trapezoidal domains // Comput. Technologies. 2017. Vol. 22, No. 4. P. 22-44. (In Russ.)
[9] Беляев В.А., ШШапеев В.П. Варианты метода коллокации и наименьших невязок для решения задач математической физики в выпуклых четырехугольных областях // Моделирование и анализ информ. систем. 2017. Т. 24, № 5. С. 629-648.
Belyaev, V.A., Shapeev, V.P. Versions of the collocation and least residuals method for solving problems of mathematical physics in the convex quadrangular domains // Modelling and Analysis of Inform. Systems. 2017. Vol. 24, No. 5. P. 629-648. (In Russ.)
[10] Belyaev, V.A., Shapeev, V.P. Versions of the collocation and least squares method for solving biharmonic equations in non-canonical domains // AIP Conf. Proceedings. 2017. Vol. 1893. P. 030102-1-030102-14. DOI: 10.1063/1.5007560.
[11] Беляев В.В., Ш^апеев В.П. Метод коллокаций и наименьших квадратов на адаптивных сетках в области с криволинейной границей // Вычисл. технологии. 2000. Т. 5, № 4. С. 13-21.
Belyaev, V.V., Shapeev, V.P. The collocation and least squares method on adaptive grids in a domain with a curvilinear boundary // Comput. Technologies. 2000. Vol. 5, No. 4. P. 13-21. (In Russ.)
[12] Исаев В.И., ШШапеев В.П. Варианты метода коллокаций и наименьших квадратов повышенной точности для численного решения уравнений Навье—Стокса // Журн. вычисл. математики и матем. физики. 2010. Т. 50, № 10. С. 1758-1770.
Isaev, V.I., Shapeev, V.P. High-accuracy versions of the collocations and least squares method for the numerical solution of the Navier—Stokes equations // Comput. Math. and Math. Phys. 2010. Vol. 50, No. 10. P. 1670-1681.
[13] Слепцов А.Г. Коллокационно-сеточное решение эллиптических краевых задач // Моделирование в механике: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ВЦ; ИТПМ. 1991. Т. 5(22), № 2. С. 101-126.
Sleptsov, A.G. Solution of elliptic boundary-value problems // Modelirovanie v Mekhanike. 1991. Vol. 5(22), No. 2. P. 101-126. (In Russ.)
[14] ШШапеев В.П., Черепанов А.Н. Конечно-разностный алгоритм для численного моделирования процессов лазерной сварки металлических пластин // Вычисл. технологии. 2006. Т. 11, № 4. С. 102-117.
Shapeev, V.P., Cherepanov, A.N. Finite-difference algorithm for numerical simulation of process of metal plates laser welding // Comput. Technologies. 2006. Vol. 11, No. 4. P. 102-117. (In Russ.)
[15] Барахнин В.Б., Шапеев В.П. Введение в численный анализ. СПб.; М.; Краснодар: Лань, 2005. 107 c.
Barakhnin, V.B., Shapeev, V.P. Introduction to numerical analysis. St. Petersburg; Moscow; Krasnodar: Lan, 2005. 107 p. (In Russ.)
[16] Saad, Y. Numerical methods for large eigenvalue problems. Manchester Univ. Press, 1991. 358 p.
[17] Слепцов А.Г. Об ускорении сходимости линейных итераций. II // Моделирование в механике: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ВЦ; ИТПМ. 1989. Т. 3, № 5. C. 118-125. Sleptsov, A.G. On acceleration for convergence of linear iterations. II // Modelirovanie v Mekhanike. 1989. Vol. 3, No. 5. P. 118-125. (In Russ.)
[18] Ворожцов Е.В., Шапеев В.П. Об ускорении итерационных процессов решения краевых задач комбинированием методов Крылова и Федоренко // Символ науки. 2015. № 10-2. С. 24-43.
Vorozhtsov, E.V., Shapeev, V.P. On acceleration of iterative processes for solving boundary value problems by combining Krylov and Fedorenko methods // Simvol Nauki. 2015. No. 10-2. P. 24-43. (In Russ.)
[19] ШШапеев В.П., Ворожцов Е.В. О комбинировании различных методов ускорения при итерационном решении уравнений с частными производными методом коллокаций и наименьших невязок // Моделирование и анализ информ. систем. 2017. Т. 24, № 1. С. 39-63. Shapeev, V.P., Vorozhtsov, E.V. On combining different acceleration techniques at the iterative solution of PDEs by the method of collocations and least residuals // Modelling and Analysis of Inform. Systems. 2017. Vol. 24, No. 1. P. 39-63. (In Russ.)
[20] Оганесян Л.А., Руховец Л.А. Вариационно-разностные методы решения эллиптических уравнений. Ереван: Изд. АН Арм. ССР, 1979. 236 с.
Oganesyan, L.A., Rukhovets, L.A. Variational-difference methods for solving elliptic equations. Yerevan: Ed. AN Arm. SSR, 1979. 236 p. (In Russ.)
[21] Ш^апеев В.П., ШШапеев А.В. Решение эллиптических задач с особенностями по схемам высокого порядка аппроксимации // Вычисл. технологии. 2006. Т. 11. Спецвыпуск, посвященный 85-летию со дня рождения Н.Н. Яненко. Ч. 2. С. 84-91.
Shapeev, V.P., Shapeev, A.V. Solutions of the elliptic problems with singularities using finite difference schemes with high order of approximation // Comput. Technologies. 2006. Vol. 11. Special issue, devoted to N.N. Yanenko's 85-th anniversary. Pt 2. P. 84-91. (In Russ.)
[22] Botella, O., Peyret, R. Benchmark spectral results on the lid-driven cavity flow // Computers & Fluids. 1998. Vol. 27, No. 4. P. 213-221.
Поступила в 'редакцию 18 октября 2017 г., с доработки —19 февраля 2018 г.
Solving the Dirichlet problem for the Poisson equation by the least squares collocation method with given discrete boundary domain
Belyaev, Vasiliy A.1'*, Shapeev, Vasiliy P.1,2
1Khristianovich Institute of Theoretical and Applied Mechanics SB RAS, Novosibirsk, 630090, Russia 2Novosibirsk State University, Novosibirsk, 630090, Russia
* Corresponding author: Belyaev, Vasiliy A., e-mail: [email protected]
The paper addresses a new version of the least squares collocation (LSC) method proposed and implemented for the numerical solution of boundary value problems for the Poisson's equation in case of given discrete domain boundary. The computer program builds a continuous double spline when boundary is smooth or a piecewise smooth double spline if boundary has salient points. In this version of the method we apply the idea of using parts of the cells of a regular grid (outside the domain). These parts of the cells are cut off by the boundary to construct the LSC method. It is assumed that the solution has no singularities at the boundary and in a certain small neighborhood of it. The differential equation for the problem is correct not only in the computational domain but also in a small neighborhood of its boundary. Then the idea of attaching "small" irregular cells to neighboring ones is used in the work with the aim of reducing the number of conditionality of the global system of linear algebraic equations. It is shown that the approximate solutions obtained by the LSC method converge with an increased order and coincide with the analytical solutions of the test problems with high accuracy.
Keywords: least squares collocation method, boundary value problem, Poisson's equation, high order approximation, discrete domain boundary, cubic spline.
Received 18 October 2017 Received in revised form 19 February 2018
© ICT SB RAS, 2018
