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

Вычислительные технологии

Том 5, № 4, 2000

МЕТОД КОЛЛОКАЦИЙ И НАИМЕНЬШИХ КВАДРАТОВ НА АДАПТИВНЫХ СЕТКАХ В ОБЛАСТИ С КРИВОЛИНЕЙНОЙ ГРАНИЦЕЙ *

В. В. Беляев, В. П. Шапеев Институт теоретической и прикладной механики СО РАН

Новосибирск, Россия e-mail: [email protected]

A version of the grid-projection method for solving boundary problems for elliptic equations in a domain with a curvilinear boundary on the rectangular adaptive grid is presented. The method is tested on the equations with a small parameter at highest derivatives. Convergence of the method and the grid adaptability to the peculiarity of the solution are shown.

Введение

Необходимость численного решения краевых задач математической физики с хорошей точностью стимулирует создание и использование новых методов и алгоритмов. В новых сеточных методах предлагается повышенная точность аппроксимации краевых условий. При ограниченных ресурсах ЭВМ точность численного решения в сеточных методах можно повысить за счет применения адаптивных сеток. Такие сетки сгущаются в подобластях особенностей решения задачи, где погрешность на равномерной сетке больше, чем в других частях расчетной области.

Здесь предлагается вариант метода коллокаций и наименьших квадратов (КНК) на прямоугольной адаптивной сетке в области с криволинейной границей. В данном способе реализации краевых условий сторона граничной ячейки совпадает с границей области.

Метод КНК численного решения краевых задач для одного уравнения второго порядка развит в работах [1,2], где показана успешность его применения к задачам с особенностями решения. Использование метода КНК для решения задач динамики вязкой жидкости [3,4] также демонстрирует его хорошие возможности. Существуют, кроме того, примеры использования проекционно-сеточного метода для решения обыкновенных дифференциальных уравнений [5].

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №99-01-00515.

© В. В. Беляев, В. П. Шапеев, 2000.

1. Метод коллокаций и наименьших квадратов

Рассматривается задача Дирихле для эллиптического уравнения в области с криволинейной границей П:

«11^X1X1 + «12^X1X2 + «22^X2X2 + Мхх + &2^Х2 + С^ = /, (жЬ Ж2) € П,

НдП = 0(Ж1,Ж2).

Здесь а11 = а11(ж1,ж2), а12 = а12(ж1,ж2), ..., / = /(ж1,ж2). Область покрывается сеткой с прямоугольными ячейками. При этом возникают регулярные (внутренние) и нерегулярные (граничные) ячейки. Рассмотрим г-ю ячейку сетки. Пусть (ж^, ж2^) — центр ячейки, к — полуширина, у1 = (ж1 — Жн)/к, у2 = (ж2 — )/к — локальные координаты, в которых ячейка имеет вид [—1,1] х [—1,1]. В г-й ячейке приближенное решение ищется в базисе полиномов второго порядка в виде

а(Ж1,Ж2) = и(у1,у2) = Ри + Р2гШ + Р3г^2 + Р4г(^12 — У22) + РбгШШ + РбгЫ + (1)

Коэффициенты гармонической части разложения (1) р^, , определяются с

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

Уравнения коллокации выписываются в четырех точках , ), ] = 1, ... , 4, и имеют вид

(аиМу1У1 + Й12«У1У2 + «22 иу2У2 + кб1«у1 + кб2«у2 + к2си)(£у ) = к2/(ху , Х2]),

где Хк] = Хк] + к£ку, к = 1, 2 и коэффициенты а11, а12, ... , с вычисляются в точках (ху, Х2]). Уравнения коллокации линейны относительно и, ... , & и могут быть представлены в виде

]>>]к№ = Р3, (2)

к=1

где коэффициенты Вук находятся дифференцированием левой части ^-го уравнения по ркг в точке ).

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

ди+ + ди- -

---+ пи+ = ^--+ Пи .

дп дп

Здесь п — внешняя нормаль к границе г-й ячейки, и+ и и- — пределы функции и при стремлении изнутри и извне к границе ячейки, п — параметр метода, зависящий от размеров ячейки и влияющий на обусловленность системы линейных уравнений, определяющих приближенное решение.

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

1/(п + 1).

Уравнения, возникающие из условий согласования, записываются в виде (2), где Вук находится дифференцированием выражения

1 /ди+ + + пи+

П + 1 \дп

6

по ркг в двух точках на каждой границе, пг — нормаль к границе г-й ячеики в ]-й точке, а Е] находится как значение в ]-й точке выражения

1

П + 1 \ дп

ди-

+ пи-

В случае, когда сторона ячейки совпадает с границей области, в двух ее точках выписываются краевые условия. Если ) — локальные координаты такой точки, то уравнение, выписываемое в ней, имеет вид (2), где Б]1 = 1, Б]2 = , Бу3 = ^, Б]4 = — ^, = 6], Б]6 = + , Е] = §(¿1, ¿2), ¿1 = К + ХН, ¿2 = К + ^ Всего для г-й ячейки выписываются четыре уравнения вида (2) из уравнений колло-кации и восемь — из условий согласования и краевых условий. При этом получается переопределенная система из 12 линейных уравнений относительно р1г, ..., р6г, из которой методом наименьших квадратов (МНК) находятся коэффициенты полинома. Это вызвано тем, что, как правило, МНК (требование минимума функционала невязки) позволяет уменьшить абсолютное значение погрешности в задачах аппроксимации. Составим из сумм квадратов невязок этих уравнений два функционала:

4

ф1(р) = £

]=1 12

6

Б]кРкг —

Ф2(Р) = £

]=5

к=1 6

У^ Б]кРкг —

к=1

Минимизируя Ф1() по 6г при фиксированных ц, а Ф2() по ц при фиксированном 6г, I = 1, ... , 5, получим систему линейных уравнений для определения кг, к = 1, ... , 6:

Е

к=1

12

7 у Бз1Бзк ]=5

Ркг

12

]=5

I = 1, ..., 5,

6

Е

к=1

У", Б]бБ]к

.3 = 1

Ркг = £ Бз6Fз. (3) ]=1

Решение отыскивается итерациями по области; при нахождении приближенного решения в г-й ячейке считается, что в остальных ячейках оно уже найдено (коэффициенты Ркт, к =1,..., 6, т = г берутся с текущей итерации, если на ней решение в т-й ячейке уже вычислялось, или с предыдущей итерации).

2

2

6

4

4

2. Адаптация сетки

После нахождения приближенного решения v(x1,x2) производится апостериорная оценка погрешности. В каждой ячейке с помощью метода коллокаций и наименьших квадратов вычисляется "уточненное" решение в виде полинома третьей степени /ш(х1,х2). В качестве оценки погрешности рассматривается функция

е(х1,х2) = |^(хьх2) — а(хьх2)|.

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

3. Регулярные и нерегулярные ячейки

В квадратной ячейке (рис. 3) в качестве точек коллокации берутся точки (\/2/2, ±\/2/2), (—\/2/2, ±\/2/2), а условия на границах ячейки ставятся в восьми точках (по две на каждой границе). Так, на нижней границе — в точках (—С, —1) и (£, —1), аналогично для других границ.

Назовем эти точки соответственно:

— обычными точками коллокации (на рисунке они обозначены символами "о");

— обычными точками согласования — на границе между ячейками (на рисунке — символы "|"), £ = л/3/3;

— обычными точками краевых условий — граница ячейки является границей области (на рисунке — символы " х"), £ = л/2/2 .

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

1. Если ячейка К содержит менее двух обычных точек согласования или не содержит обычных точек коллокации (рис. 2), то решение для нее не определяется (она не рассматривается как самостоятельная), а за счет ее части, лежащей внутри области, дополняется соседняя ячейка Ь. Для нахождения решения в ячейке Ь вместо условий согласования между ячейками К и Ь используются краевые условия на границе области. Из обычных точек краевых условий на границе между К и Ь перпендикулярно ей проводятся прямые, и в точках их пересечения с границей области ставятся краевые условия. Таким образом, удается избежать существования сильно вытянутых ячеек, что обычно является одной из причин плохой обусловленности системы (3).

Рис. 1. Рис. 2.

2. Для ячейки, центр которой лежит внутри области достаточно далеко (на расстоянии большем Н/2) от ее границы, точки постановки краевых условий определяются следующим способом (рис. 3, а):

— обычные точки краевых условий, лежащие вне области, соединяются отрезками прямых с центром ячейки;

— краевые условия ставятся в точках пересечений этих отрезков с границей области.

Подобным образом находятся и точки коллокации (рис. 3, б). Центр ячейки соединяется

отрезками с углами ячейки, лежащими вне области. На каждом отрезке выбирается точка, делящая часть отрезка между центром ячейки и границей области в том же отношении, в котором обычная точка коллокации делит отрезок, соединяющий центр и вершину угла ячейки.

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

3. Центр ячейки лежит за пределами области, но часть ячейки, лежащая внутри области, содержит хотя бы одну точку коллокации и не менее двух обычных точек согласования. В ячейке, которая содержит две обычные точки коллокации, лежащие внутри области (рис. 4, а), выбираются еще две — с координатами (±£/2,£), где ( = —\/2/2. Аналогично для случая, показанного на рис. 4, б.

Иначе (рис. 5, а), уравнения коллокации выписываются в точках (С, С), (С, С), (С, С),

(С, С), где С = —^2/2, С = (С — 1)/2.

Если ячейка содержит четыре обычных точки согласования (рис. 5, б), то условия согласования выписываются только в них. В случае, когда в ячейке три обычных точки согласования (рис. 6, а), кроме них берется еще одна — (—1, —(1 + \/3)/3). Если ячейка содержит только две обычных точки согласования (рис. 6, б), то на нижней и левой границах ставится по два условия согласования — в обычных точках согласования и, дополнительно , в точках (—1, — (1+л/3)/3) и (—(1+л/3)/3, — 1). В четырех точках, равномерно

Рис.

4.

Рис. 6.

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

4. Тестирование

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

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

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

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

Рис. 7.

В качестве примера использовалась модельная задача типа диффузии-конвекции, рассмотренная в [2]. Уравнение

eAu + (xi — ai)ux1 + (x2 — a2)ux2 = 0

(4)

имеет семейство точных решений. Если для произвольного угла р ввести обозначения ц = cosр, V = sinр, p = — (ца1 + va2), q = vai — ца2,

yi = ц/Х\ + vx2 + p, y2 = — vxi + цх2 + q, G(y) = 0.5+—= e-í2dt,

Vn J

0

то функция u(xi, x2) = G(yi)G(y2) будет являться решением уравнения (4). Уравнение (4) решалось в четверти единичного круга с краевыми условиями

и(х 1 ,Х2 )|ап = С(у1)С(у2). Если параметр £ принимает малые значения, то функция и(х1,х2) имеет внутренние

Таблица 1

N 400 1600 6400 25600

A 3.5 ■ 10-5 4.4 ■ 10-6 6.6 ■ 10-7 3.2 ■ 10-7

B 3.7 ■ 10-5 4.6 ■ 10-6 7.3 ■ 10-7 4.0 ■ 10-7

слои, что существенно влияет на точность расчетов. В табл. 2 представлена зависимость погрешности от числа ячеек сетки для случая е = 10-4, ^ = п/6, а1 = 0.6, = 0.3. Для сравнения в табл. 3 приведена зависимость погрешности от числа ячеек при е = 1. В этом случае решение не имеет больших градиентов и заданная точность достигается на значительно меньшем числе ячеек. На рис. 8, а изображена адаптивная сетка, построенная программой при отыскании численного решения задачи, на рис. 8, б — график точного решения уравнения (4), имеющего ярко выраженный внутренний слой вдоль кривой в расчетной области. Данная сетка содержит ячейки восьми разных размеров, то есть длины сторон наибольшей и наименьшей ячеек отличаются в 256 раз. Видно, что алгоритм, реализованный в программе, измельчает расчетную сетку на участках с большими градиентами решения. В другом примере решалось уравнение Ди = 19000х18 + 6у с точным решением и = 50х20 + у3 — у, из которого брались краевые условия. Решение в заданной области принимает значения от —1 до 50, а его градиенты изменяются от 0 до 1000.

Т а б л и ц а 2

N 16 31 367 2083 6724 8446 17968 29680 69718

Error 0.63 0.59 0.46 0.22 4.9 ■ 10-2 1.9 ■ 10-2 8.3 ■ 10-3 3.5 ■ 10-3 5.0 ■ 10-4

Т а б л и ц а 3

N 16 61 232 862 3118 6241 16171

Error 5.5 ■ 10-4 2.0 ■ 10-5 3.7 ■ 10-6 1.9 ■ 10-6 7.9 ■ 10-7 3.4 ■ 10-7 2.5 ■ 10-8

В табл. 4 представлена зависимость величины погрешности решения от числа ячеек N для адаптивной и равномерной сеток. Анализ таблицы показывает, что при одинаковой погрешности адаптивная сетка содержит на порядок меньшее число ячеек, чем равномерная. Использование адаптивной сетки позволяет существенно снизить количество вычислений, необходимых для получения решения с заданной точностью, что дает возможность решать сложные задачи при ограниченных мощностях вычислительной техники.

Т а б л и ц а 4

Адаптивная сетка Равномерная сетка

N Error N Error

28 9.77 61 9.77

67 2.15 232 2.15

127 7.0 ■ 10-1 880 5.55 ■ 10-1

379 1.85 ■ 10-1 3400 1.48 ■ 10-1

862 8.83 ■ 10-2 13288 3.84 ■ 10-2

1495 4.75 ■ 10-2 52366 1.02 ■ 10-2

4753 1.29 ■ 10-2

5617 1.17 ■ 10-2

15823 4.12 ■ 10-3

Рис. 8.

Список литературы

[1] Слепцов А. Г. Коллокационно-сеточное решение эллиптических краевых задач. Моделирование в механике, 5(22), №2, Новосибирск, 1991, 101-126.

[2] Слепцов А. Г., Шокин Ю. И. Адаптивный проекционно-сеточный метод для эллиптических задач. Журн. вычисл. матем. и матем. физики, 37, №5, 1997, 572-586.

[3] Семин Л. Г., Слепцов А. Г., Шапеев В. П. Метод коллокаций — наименьших квадратов для уравнений Стокса. Вычисл. технологии, 1, №2, 1996, 90-98.

[4] Семин Л. Г., Шапеев В. П. Метод коллокаций и наименьших квадратов для уравнений Навье —Стокса. Там же, 3, №3, 1998, 72-84.

[5] Ascher U.M., Mattheij R. M.M., Russell R. D. Numerical solution of boundary value problems for ordinary differential equations. Prentice Hall Series in Comput. Math., 1988.

Поступила в редакцию 6 января 2000 г.