1. На главную
  2. Электронные версии
  3. Журнал вычислительной математики и математической физики
  4. T. 62, № 4, 2022

Журнал вычислительной математики и математической физики, 2022, T. 62, № 4, стр. 531-552

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

В. А. Беляев 1, 2, *, Л. С. Брындин 1, 2, **, С. К. Голушко 2, 3, ***, Б. В. Семисалов 2, 4, ****, В. П. Шапеев 1, 2, *****

1 Институт теоретической и прикладной механики им. С.А. Христиановича СО РАН
630090 Новосибирск, ул. Институтская, 4/1, Россия

2 Новосибирский государственный университет
630090 Новосибирск, ул. Пирогова, 2, Россия

3 Федеральный исследовательский центр информационных и вычислительных технологий
630090 Новосибирск, пр-т Акад. Лаврентьева, 6, Россия

4 Институт математики им. С.Л. Соболева СО РАН
630090 Новосибирск, пр-т Акад. Коптюга, 4, Россия

* E-mail: belyaevasily@mail.ru
** E-mail: bryndin-1996@mail.ru
*** E-mail: s.k.golushko@gmail.com
**** E-mail: vibis@ngs.ru
***** E-mail: shapeev.vasily@mail.ru

Поступила в редакцию 10.02.2020
После доработки 05.03.2021
Принята к публикации 16.11.2021

Полный текст (PDF)

Аннотация

Предложены и реализованы новые h-, p- и hp-варианты метода коллокации и наименьших квадратов, позволяющие находить приближенные решения краевых задач для неоднородного бигармонического уравнения в нерегулярных и многосвязных областях. Приведены формулы для операции продолжения при переходе с грубой сетки на более мелкую на многосеточном комплексе в случае применения различных пространств полиномов. Экспериментально показано, что численные решения краевых задач, полученные разработанными вариантами метода, сходятся с повышенным порядком к аналитическим решениям в случаях, когда последние не имеют особенностей. Приведено сравнение полученных результатов с результатами других авторов, использовавших конечно-разностный, конечно-элементный и другие методы, основанные на применении полиномов Чебышёва. Рассмотрены примеры задач с особенностями. Разработанные варианты метода использованы для моделирования изгиба упругой изотропной пластины нерегулярной формы, находящейся под действием поперечной нагрузки. Библ. 30. Фиг. 5. Табл. 12.

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

ВВЕДЕНИЕ

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

Данная работа посвящена численному решению с повышенной точностью различных краевых задач для бигармонических уравнений, к решению которых сводится анализ поведения многих задач механики сплошной среды (МСС): гидродинамики при малых числах Рейнольдса, линейной теории упругости, теории тонких пластин и др.

При численном решении краевых задач для бигармонического уравнения во многих случаях применяется метод конечных разностей (МКР), главными преимуществами которого являются простота построения расчетной сетки и скорость решения соответствующих задач линейной алгебры. При этом существует ограниченное количество публикаций, посвященных применению МКР для решения бигармонического уравнения в нерегулярных областях (см. [1]–[4]). Другим широко используемым методом для решения таких задач является метод конечных элементов (МКЭ) с использованием неструктурированных сеток (см. [5]–[7]). Как отмечается в [8], несмотря на то, что техника решения УЧП четвертого порядка с помощью МКЭ и МКР хорошо развита, существует не так много результатов, полученных при решении краевых задач в нерегулярных областях со сложными граничными условиями. Кроме того, работа с несколькими тысячами элементов в расчетной области, имеющих разную форму и ориентацию, – достаточно сложная задача (см. [8], [9]).

Условно численные методы решения бигармонического уравнения можно разбить на две основные группы: методы, использующие непосредственную дискретизацию заданного уравнения и краевых условий (см. [8], [10], [11]), и смешанные (см. [3], [4]). Например, в [3], [4] отмечается сложность численного решения бигармонического уравнения из-за наличия производных четвертого порядка в дифференциальном уравнении, поскольку, как известно, переход от системы дифференциальных уравнений к одному уравнению более высокого порядка сопровождается существенным увеличением обусловленности системы линейных алгебраических уравнений (СЛАУ) приближенной задачи. Тем не менее, как сказано в [12], “…построение эффективных методов численного решения бигармонического уравнения может уменьшить влияние этого обстоятельства”. Чтобы избежать проявлений этого недостатка, в [3] вводится дополнительная переменная и решение исходной задачи Дирихле для бигармонического уравнения сводится к последовательному решению двух уравнений Пуассона.

В последнее время при решении задач МСС все более популярными становятся спектральные методы, методы спектральных элементов и методы, основанные на применении полиномов Чебышева достаточно высоких степеней, главным преимуществом которых является экспоненциальный характер уменьшения погрешности в случае достаточно гладких решений (см. [8], [10], [11]). Однако, когда решения не являются достаточно гладкими, такие методы становятся менее эффективными и менее гибкими, например, по сравнению с МКЭ. Существуют и другие методы решения бигармонического уравнения: метод быстрых мультиполей (см. [13]), метод сплайн-коллокации (см. [14]), миметический метод (см. [15]), смешанные методы (см. [9]) и др. Отметим, что большинство работ посвящено решению краевой задачи Дирихле для бигармонического уравнения и существенно меньше задачам, в краевые условия которых входят производные более высоких порядков. Их наличие приводит к дополнительным трудностям и уменьшает возможности получения решений высокой точности.

Данная работа посвящена разработке и реализации h-, p- и hp-вариантов метода коллокации и наименьших квадратов (h-, p- и hp-МКНК) (см. [16]–[24]) для решения бигармонических уравнений в нерегулярных областях, в том числе в многосвязных областях. В проекционно-сеточном МКНК путем проектирования исходной задачи для УЧП в конечномерное линейное функциональное пространство ставится в соответствие приближенная задача, сводящаяся к решению переопределенной СЛАУ, которое определяет приближенное решение дифференциальной задачи. В h-МКНК уточнение решения происходит за счет измельчения шагов сетки при фиксированной степени базисных полиномов используемого полиномиального пространства. В p-МКНК уточнение решения происходит за счет увеличения степени базисных полиномов используемого полиномиального пространства при фиксированной величине шага сетки. В hp-МКНК уточнение решения происходит за счет измельчения шагов сетки и увеличения степени базисных полиномов используемого полиномиального пространства. В данной работе исследованы возможности p-МКНК1 при решении задачи в пространстве полиномов достаточно высокой степени в одной ячейке, совпадающей со всей расчетной областью, и p-МКНКh с фиксированным количеством ячеек в количестве больше одной.

1. ПОСТАНОВКА ЗАДАЧИ И СЛУЧАЙ ДИСКРЕТНОГО ЗАДАНИЯ ВНЕШНЕЙ ГРАНИЦЫ ОБЛАСТИ

Рассмотрим краевую задачу Дирихле для неоднородного бигармонического уравнения

(1.1)
${{\Delta }^{2}}u = f({{x}_{1}},{{x}_{2}}),\quad ({{x}_{1}},{{x}_{2}}) \in \Omega ,$
(1.2)
${{\left. u \right|}_{{{{\gamma }_{i}}}}} = {{g}_{{1i}}}({{x}_{1}},{{x}_{2}}),\quad i = \overline {0,{{N}_{\gamma }}} ,$
(1.3)
${{\left. {{{u}_{n}}} \right|}_{{{{\gamma }_{0}}}}}\mathop = \limits^{{\text{п}}{\text{.в}}{\text{.}}} {{g}_{{20}}}({{x}_{1}},{{x}_{2}}),$
(1.4)
${{\left. {{{u}_{n}}} \right|}_{{{{\gamma }_{i}}}}} = {{g}_{{2i}}}({{x}_{1}},{{x}_{2}}),\quad i = \overline {1,{{N}_{\gamma }}} ,$
в нерегулярной многосвязной области $\Omega \subset {{\mathbb{R}}^{2}}$ с внешней замкнутой границей ${{\gamma }_{0}}$ и внутренними границами ${{\gamma }_{i}}$ (окружности), $i = \overline {1,{{N}_{\gamma }}} $ (фиг. 1), где ${{N}_{\gamma }}$ – количество круглых отверстий внутри области $\Omega $. Здесь $u({{x}_{1}},{{x}_{2}})$ – искомая функция, $f({{x}_{1}},{{x}_{2}})$, ${{g}_{{1i}}}({{x}_{1}},{{x}_{2}})$ и ${{g}_{{2i}}}({{x}_{1}},{{x}_{2}})$ – заданные функции, ${{\Delta }^{2}} = \frac{{{{\partial }^{4}}}}{{\partial x_{1}^{4}}} + 2\frac{{{{\partial }^{4}}}}{{\partial x_{1}^{2}\partial x_{2}^{2}}} + \frac{{{{\partial }^{4}}}}{{\partial x_{2}^{4}}}$, ${{u}_{n}} = \frac{{\partial u}}{{\partial n}}$, ${\mathbf{n}}$ – внешняя единичная нормаль к соответствующей границе ${{\gamma }_{i}}$, $i = \overline {0,{{N}_{\gamma }}} $. Решение задачи (1.1)–(1.4) понимается в обобщенном смысле. Предполагается, что $u \in {{C}^{4}}(\Omega )$, ${{\left. u \right|}_{{{{\gamma }_{i}}}}} \in {{C}^{1}}({{\gamma }_{i}})$, $i = \overline {1,{{N}_{\gamma }}} $, и $u$ – непрерывная кусочно-гладкая функция, дифференцируемая как минимум один раз на каждой из частей ${{\gamma }_{0}}$, отделенных друг от друга точками излома. Условие (1.2) ставится всюду на ${{\gamma }_{i}}$, $i = \overline {0,{{N}_{\gamma }}} $, и (1.4) всюду на ${{\gamma }_{i}}$, $i = \overline {1,{{N}_{\gamma }}} $. При построении обобщенного решения краевое условие (1.3) требуется выполнить почти всюду (п.в.) на ${{\gamma }_{0}}$ (см. [26]): всюду на ${{\gamma }_{0}}$ за исключением точек излома в случае их наличия, поскольку направление нормали в них не определено.

Фиг. 1.

Область решения задачи, где $t$ – номера точек “$ \bullet $”, дискретно задающих ее границу.

На практике зачастую встречаются задачи, граница которых полностью или частично задана дискретным набором точек. Пусть внешняя граница ${{\gamma }_{0}}$ области $\Omega $ задана дискретно конечной последовательностью точек, расположенных приблизительно равномерно относительно друг друга (фиг. 1). Восполним в этом случае границу области по входным дискретным данным векторной сплайн-функцией, в качестве компонент которой взяты два кубических сплайна, являющиеся функциями параметра $t$, монотонно меняющегося при монотонном движении точки вдоль замкнутой границы. Во многих случаях в качестве параметра $t$, $t = \overline {0,{{N}_{p}} - 1} $, удобно взять номер точки (фиг. 1). Более подробное описание техники восполнения границы области можно найти в [17]. Отметим, что если из входных данных или других каких-то соображений известно, что некоторые точки являются точками излома границы, то следует строить отдельные векторные сплайн-функции на участках, концами которых являются точки излома.

2. ОПИСАНИЕ МКНК

2.1. Построение сеток и выбор базисных элементов

Расчетную область $\Omega $ впишем в прямоугольник размера ${{d}_{1}} \times {{d}_{2}}$ (см. [17]). В нем построим сетку ${{N}_{1}} \times {{N}_{2}}$ с прямоугольными ячейками ${{\Omega }_{j}}$ размера $2{{h}_{1}} \times 2{{h}_{2}}$, $j = \overline {1,{{N}_{{{\text{cells}}}}}} $, где ${{N}_{{{\text{cells}}}}}$ – количество расчетных ячеек.

Ячейки сетки, полностью расположенные внутри области, назовем внутренними (фиг. 2, ячейки 9 и 11). Прямоугольные ячейки, пересеченные границами области ${{\gamma }_{i}}$, $i = \overline {0,{{N}_{\gamma }}} $, назовем граничными. Часть граничной ячейки, отсеченную границами ${{\gamma }_{i}}$, $i = \overline {0,{{N}_{\gamma }}} $, и лежащую внутри области, назовем н-ячейкой (фиг. 2, ячейки 1–8, 10, 12–16). Считаем, что на фиг. 2 номера граничных ячеек совпадают с номерами их внутренних н-ячеек. Часть граничной ячейки, расположенную вне области, назовем законтурной. Прямоугольную ячейку сетки, от которой границами области ${{\gamma }_{i}}$, $i = \overline {0,{{N}_{\gamma }}} $, отсечена н-ячейка, назовем материнской. Внутренние прямоугольные ячейки полностью совпадают со своими материнскими ячейками. Стороны любых ячеек, расположенные внутри области, назовем внутренними сторонами ячеек. Части границ ${{\gamma }_{i}}$, $i = \overline {0,{{N}_{\gamma }}} $, которые оказались внутри граничной ячейки, назовем ее внешними сторонами.

Фиг. 2.

Фрагмент расчетной области (a) и его различные части (б)–(г) при $K = 4$; значок “$ \bullet $” обозначает точки, дискретно задающие границу области, “$\diamondsuit $” – точки коллокации, “$ \times $” – точки согласования, “$\square $” – точки записи краевых условий, “$ \circ $” – центры некоторых н-ячеек. Темным цветом выделены материнские ячейки, содержащие внутри себя н-ячейки 1–6, которые присоединяются к ячейкам 7–11 соответственно. У граничных ячеек 2, 3, 5 тонкими стрелочками указаны их нерегулярные части, т.е. сами н-ячейки.

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

(2.1)
${{y}_{1}} = ({{x}_{1}} - {{x}_{{1j}}}){\text{/}}{{h}_{1}},\quad {{y}_{2}} = ({{x}_{2}} - {{x}_{{2j}}}){\text{/}}{{h}_{2}},$
где $({{x}_{{1j}}},{{x}_{{2j}}})$ – центр $j$-й материнской ячейки, ${{h}_{1}} = {{d}_{1}}{\text{/}}(2{{N}_{1}})$, ${{h}_{2}} = {{d}_{2}}{\text{/}}(2{{N}_{2}})$. Положим далее $v({{y}_{1}},\;{{y}_{2}})$ = $u({{x}_{1}}({{y}_{1}}),{{x}_{2}}({{y}_{2}}))$.

Задача (1.1)–(1.4) после замены (2.1) в локальных переменных примет следующий вид:

$\frac{1}{{h_{1}^{4}}}\frac{{{{\partial }^{4}}v}}{{\partial {{y}_{1}}^{4}}} + \frac{2}{{h_{1}^{2}h_{2}^{2}}}\frac{{{{\partial }^{4}}v}}{{\partial {{y}_{1}}^{2}\partial {{y}_{2}}^{2}}} + \frac{1}{{h_{2}^{4}}}\frac{{{{\partial }^{4}}v}}{{\partial {{y}_{2}}^{4}}} = f({{x}_{1}}({{y}_{1}}),{{x}_{2}}({{y}_{2}})),\quad ({{y}_{1}},{{y}_{2}}) \in \Omega ,$
(2.2)
$\begin{gathered} {{\left. v \right|}_{{{{\gamma }_{i}}}}} = {{g}_{{1i}}}({{x}_{1}}({{y}_{1}}),{{x}_{2}}({{y}_{2}})),\quad ({{y}_{1}},{{y}_{2}}) \in {{\gamma }_{i}} \cap \delta {{\Omega }_{j}},\quad i = \overline {0,{{N}_{\gamma }}} , \\ {{\left. {\left( {\frac{{{{n}_{1}}}}{{{{h}_{1}}}}\frac{{\partial v}}{{\partial {{y}_{1}}}} + \frac{{{{n}_{2}}}}{{{{h}_{2}}}}\frac{{\partial v}}{{\partial {{y}_{2}}}}} \right)} \right|}_{{{{\gamma }_{0}}}}}\mathop = \limits^{{\text{п}}{\text{.в}}{\text{.}}} {{g}_{{20}}}({{x}_{1}}({{y}_{1}}),{{x}_{2}}({{y}_{2}})),\quad ({{y}_{1}},{{y}_{2}}) \in {{\gamma }_{0}} \cap \delta {{\Omega }_{j}}, \\ \end{gathered} $
${{\left. {\left( {\frac{{{{n}_{1}}}}{{{{h}_{1}}}}\frac{{\partial v}}{{\partial {{y}_{1}}}} + \frac{{{{n}_{2}}}}{{{{h}_{2}}}}\frac{{\partial v}}{{\partial {{y}_{2}}}}} \right)} \right|}_{{{{\gamma }_{i}}}}} = {{g}_{{2i}}}({{x}_{1}}({{y}_{1}}),{{x}_{2}}({{y}_{2}})),\quad ({{y}_{1}},{{y}_{2}}) \in {{\gamma }_{i}} \cap \delta {{\Omega }_{j}},\quad i = \overline {1,{{N}_{\gamma }}} ,$
где $({{n}_{1}},{{n}_{2}})$ – компоненты внешней единичной нормали ${\mathbf{n}}$ к соответствующей границе ${{\gamma }_{i}}$ в точке $({{x}_{1}},{{x}_{2}})$, $i = \overline {0,{{N}_{\gamma }}} $, $\delta {{\Omega }_{j}}$ – внешняя сторона $j$-й ячейки, $j = \overline {1,{{N}_{{{\text{cells}}}}}} $.

Приближенное решение ${{v}_{{hj}}}$ задачи (2.2) в каждой $j$-й ячейке ищется в виде линейной комбинации с неизвестными коэффициентами базисных элементов некоторого функционального пространства. В случае применения полиномов Чебышёва представление приближенного решения, по аналогии с [19], имеет вид

(2.3)
$\begin{gathered} {{v}_{{hj}}}({{y}_{1}},{{y}_{2}}) = \sum\limits_{{{i}_{1}} = 0}^{{{K}_{1}} - 1} \,\sum\limits_{{{i}_{2}} = 0}^{{{K}_{2}} - 1} \,{{c}_{{{{i}_{1}}{{i}_{2}},j}}}{{\phi }_{{{{i}_{1}}}}}({{y}_{1}}){{\phi }_{{{{i}_{2}}}}}({{y}_{2}}),\quad ({{y}_{1}},{{y}_{2}}) \in [ - 1,1] \times [ - 1,1], \\ {{\phi }_{{{{i}_{1}}}}}({{y}_{1}}) = \cos ({{i}_{1}}\arccos ({{y}_{1}})),\quad {{\phi }_{{{{i}_{2}}}}}({{y}_{2}}) = \cos ({{i}_{2}}\arccos ({{y}_{2}})). \\ \end{gathered} $
Если в качестве базисных элементов взяты мономы, то приближенное решение имеет вид
(2.4)
${{v}_{{hj}}}({{y}_{1}},{{y}_{2}}) = \sum\limits_{{{i}_{1}} = 0}^K \,\sum\limits_{{{i}_{2}} = 0}^{K - {{i}_{1}}} \,{{c}_{{{{i}_{1}}{{i}_{2}},j}}}y_{1}^{{{{i}_{1}}}}y_{2}^{{{{i}_{2}}}}.$
В МКНК для определения неизвестных коэффициентов ${{c}_{{{{i}_{1}}{{i}_{2}},j}}}$ в каждой ячейке выписывается переопределенная “локальная” СЛАУ, состоящая из уравнений коллокации, условий согласования на общих сторонах, принадлежащих двум соседним ячейкам, и краевых условий на ${{\gamma }_{i}}$, $i = \overline {0,{{N}_{\gamma }}} $, если ячейка является граничной. Совокупность всех “локальных” СЛАУ – глобальная СЛАУ, которая определяет глобальное решение задачи.

Если в расчетной сетке имеется или появилась после измельчения шагов сетки маленькая и/или вытянутая н-ячейка (фиг. 2, ячейки 1–6), то глобальная СЛАУ задачи, как правило, становится хуже обусловленной (см. [16], [17]). Такие ячейки могут быть причиной понижения точности приближенного решения задачи, увеличения количества итераций и времени счета, а также расходимости итерационного процесса (см. [16]). Чтобы справиться с этой проблемой, здесь используется идея присоединения таких н-ячеек к соседним ячейкам. Подробную информацию о технике присоединения можно найти, например, в [17]. Здесь кратко обозначим основную идею. Назовем несамостоятельными н-ячейки, расположенные в таких граничных ячейках, в которых начало локальной системы координат находится вне расчетной области (фиг. 2, ячейки 1–6). Все остальные ячейки, включая внутренние и некоторые исключительные ячейки, речь о которых пойдет ниже, по определению считаются самостоятельными. Предлагается несамостоятельную н-ячейку присоединять к той соседней самостоятельной ячейке, с которой она имеет наибольшей длины внутреннюю сторону среди ее всех внутренних сторон, общих с другими соседними самостоятельными ячейками (фиг. 2, ячейки 1–6 присоединяются к ячейкам 7–11 соответственно). Через такие стороны в присоединенные н-ячейки продолжается решение из ячейки, к которой их присоединили. Согласно сформулированному выше правилу, к самостоятельной ячейке могут быть присоединены несколько малых/вытянутых несамостоятельных н-ячеек. В этом случае внешняя сторона объединенной ячейки состоит из внешних сторон всех объединенных в ней ячеек. Объединенную ячейку также считаем граничной н-ячейкой. При этом, если круглое отверстие в случае многосвязной области полностью оказалось внутри прямоугольной ячейки, даже если начало локальной системы координат находится вне расчетной области $\Omega $, то такую ячейку считаем самостоятельной (фиг. 2, ячейка 16). Прямоугольная ячейка, которая изначально содержит внутри себя несколько непересекающихся участков границ ${{\gamma }_{i}}$, $i = \overline {0,{{N}_{\gamma }}} $, также считается самостоятельной.

В p-МКНК$_{1}$ исходная область также включается в прямоугольник (фиг. 3) и в расчетной области строится приближение решения единым полиномом высокой степени. Аналогично по формуле (2.1) вводятся локальные координаты и ${{N}_{{{\text{cells}}}}} = {{N}_{1}} = {{N}_{2}} = j = 1$, $2{{h}_{1}} = {{d}_{1}}$, $2{{h}_{2}} = {{d}_{2}}$.

Фиг. 3.

Двусвязная область. Схема расположения точек записи краевых условий и точек коллокации в трех случаях для p-МКНК1: (a) – в корнях полиномов Чебышёва в прямоугольнике при ${{K}_{1}} = {{K}_{2}} = 6$, (б) – равномерно в прямоугольнике при $K = 6$, (в) – равномерно внутри исходной расчетной области при $K = 6$.

2.2. Уравнения приближенной задачи и их решение

Будем понимать под степенью переопределенности “локальных” СЛАУ величину $\eta $, равную отношению количества уравнений в ней к количеству неизвестных. В работах [16], [17], [22], посвященных решению бигармонического уравнения h-МКНК было установлено, что хорошая обусловленность приближенной задачи достигается при $\eta \approx 3$.

Уравнения коллокации, умноженные на $h_{1}^{2}h_{2}^{2}$ в каждой $j$-й ячейке, $j = \overline {1,{{N}_{{{\text{cells}}}}}} $, выписываются в ${{N}_{{{\text{col}}}}}$ точках коллокации и имеют вид

(2.5)
${{k}_{c}}\left( {\frac{{h_{2}^{2}}}{{h_{1}^{2}}}\frac{{{{\partial }^{4}}{{v}_{{hj}}}}}{{\partial y_{1}^{4}}} + 2\frac{{{{\partial }^{4}}{{v}_{{hj}}}}}{{\partial y_{1}^{2}\partial y_{2}^{2}}} + \frac{{h_{1}^{2}}}{{h_{2}^{2}}}\frac{{{{\partial }^{4}}{{v}_{{hj}}}}}{{\partial y_{2}^{4}}}} \right) = {{k}_{c}}h_{1}^{2}h_{2}^{2}f({{x}_{1}}({{y}_{{1c}}}),{{x}_{2}}({{y}_{{1c}}})),$
где $({{y}_{{1c}}},{{y}_{{2c}}})$ – точки коллокации, $c = \overline {1,{{N}_{{{\text{col}}}}}} $, ${{k}_{c}}$ – положительный весовой множитель уравнения коллокации. В случае применения базисных полиномов Чебышёва количество ${{N}_{{{\text{nu}}}}}$ неизвестных коэффициентов ${{c}_{{{{i}_{1}}{{i}_{2}},j}}}$ в (2.3) будет ${{N}_{{{\text{nu}}}}} = {{K}_{1}}{{K}_{2}}$. Аналогично [19], полагаем ${{N}_{{{\text{col}}}}} = {{N}_{{{\text{nu}}}}}$ и выписываем уравнения коллокации в точках с координатами $(\alpha _{{{{i}_{1}}}}^{1},\alpha _{{{{i}_{2}}}}^{2})$, ${{i}_{1}} = \overline {1,{{K}_{1}}} $, ${{i}_{2}} = \overline {1,{{K}_{2}}} $, $\alpha _{{{{i}_{1}}}}^{1}$ и $\alpha _{{{{i}_{2}}}}^{2}$ – корни многочленов Чебышева первого рода степеней ${{K}_{1}}$ и ${{K}_{2}}$ соответственно. На фиг. 3а изображена расстановка точек коллокации в случае применения полиномов Чебышева при ${{K}_{1}} = {{K}_{2}} = 6$ для p-МКНК1.

В случае применения мономов степени $K$ в качествe базисных элементов количество ${{N}_{{{\text{nu}}}}}$ неизвестных коэффициентов ${{c}_{{{{i}_{1}}{{i}_{2}},j}}}$ в (2.4) будет ${{N}_{{{\text{nu}}}}} = (K + 1)(K + 2){\text{/}}2$. Полагаем ${{N}_{{{\text{col}}}}} = ([\sqrt {{{N}_{{{\text{nu}}}}}} ] + {{1)}^{2}}$, где “$[{\kern 1pt} \cdot {\kern 1pt} ]$” – целая часть числа. Затем располагаем точки коллокации в каждой ячейке равномерно. На фиг. 3б изображена расстановка точек коллокации в случае применения мономов при $K = 6$ для p-МКНК1, а на фиг. 2 изображены фрагмент расчетной области и его различные части с сеткой при $K = 4$. При этом во всех вариантах МКНК (кроме p-МКНК1) в оказавшихся внутри отверстия точках коллокации н-ячеек, к которым не присоединялись н-ячейки, уравнение коллокации (2.5) не выписывалось (см. фиг. 2, ячейки 13 и 16). Авторы этой работы дополнительно рассмотрели p-МКНК1 без использования законтурной части. На фиг. 3в изображена равномерная расстановка точек коллокации в случае применения мономов при $K = 6$, которая была осуществлена следующим образом. Сначала полагалось ${{N}_{{{\text{col}}}}} = ([\sqrt {{{N}_{{{\text{nu}}}}}} ] + {{1)}^{2}}$ и точки равномерно располагались (фиг. 3б) в прямоугольнике, который включает нерегулярную область. Далее компьютерная программа вычисляла количество точек, которые расположены внутри исходной нерегулярной области. Если количество таких точек оказалось меньше, чем значение ${{N}_{{{\text{nu}}}}}$, то последовательно полагалось ${{N}_{{{\text{col}}}}} = ([\sqrt {{{N}_{{{\text{nu}}}}}} ] + {{2)}^{2}},([\sqrt {{{N}_{{{\text{nu}}}}}} ] + {{3)}^{2}},...$ до тех пор, пока очередное ${{N}_{{{\text{col}}}}}$ не оказалось больше, чем ${{N}_{{{\text{nu}}}}}$.

Во всех вариантах (кроме p-МКНК1) в каждой $j$-й ячейке, $j = \overline {1,{{N}_{{{\text{cells}}}}}} $, в качестве условия согласования решения в $[{{N}_{{{\text{col}}}}}{\text{/}}4]$ точках согласования (фиг. 2) на каждой общей стороне между соседними ячейками требовалась непрерывность линейной комбинации с весами значений искомого решения и его производных по нормали:

(2.6)
${{k}_{{{{m}_{0}}}}}{{v}_{{hj}}} + \frac{{{{k}_{{{{m}_{1}}}}}}}{{{{h}_{{nj}}}}}\frac{{\partial {{v}_{{hj}}}}}{{\partial {{n}_{j}}}} = {{k}_{{{{m}_{0}}}}}{{\hat {v}}_{h}} + \frac{{{{k}_{{{{m}_{1}}}}}}}{{{{h}_{{nj}}}}}\frac{{\partial {{{\hat {v}}}_{h}}}}{{\partial {{n}_{j}}}},$
(2.7)
${{h}_{1}}{{h}_{2}}\left( {\frac{{{{k}_{{{{m}_{2}}}}}}}{{h_{{nj}}^{2}}}\frac{{{{\partial }^{2}}{{v}_{{hj}}}}}{{\partial n_{j}^{2}}} + \frac{{{{k}_{{{{m}_{3}}}}}}}{{h_{{nj}}^{3}}}\frac{{{{\partial }^{3}}{{v}_{{hj}}}}}{{\partial n_{j}^{3}}}} \right) = {{h}_{1}}{{h}_{2}}\left( {\frac{{{{k}_{{{{m}_{2}}}}}}}{{h_{{nj}}^{2}}}\frac{{{{\partial }^{2}}{{{\hat {v}}}_{h}}}}{{\partial n_{j}^{2}}} + \frac{{{{k}_{{{{m}_{3}}}}}}}{{h_{{nj}}^{3}}}\frac{{{{\partial }^{3}}{{{\hat {v}}}_{h}}}}{{\partial n_{j}^{3}}}} \right),$
где ${{{\mathbf{n}}}_{{\mathbf{j}}}}$ – внешняя единичная нормаль к границе $j$-й ячейки, ${{v}_{{hj}}}$ и ${{\hat {v}}_{h}}$ – пределы приближенного решения задачи при стремлении изнутри и извне к границе $j$-й ячейки соответственно, ${{h}_{{nj}}} = {{h}_{1}}$, если направление ${{{\mathbf{n}}}_{{\mathbf{j}}}}$ совпадает с направлением оси ${{y}_{1}}$, и ${{h}_{{nj}}} = {{h}_{2}}$, если направление ${{{\mathbf{n}}}_{{\mathbf{j}}}}$ совпадает с направлением оси ${{y}_{2}}$, ${{k}_{{{{m}_{0}}}}},{{k}_{{{{m}_{1}}}}},{{k}_{{{{m}_{2}}}}},{{k}_{{{{m}_{3}}}}}$ – положительные весовые множители условий согласования. Здесь $\frac{{\partial {{v}_{{hj}}}}}{{\partial {{n}_{j}}}} = {{n}_{{1j}}}\frac{{\partial {{v}_{{hj}}}}}{{\partial {{y}_{1}}}} + {{n}_{{2j}}}\frac{{\partial {{v}_{{hj}}}}}{{\partial {{y}_{2}}}}$, где $({{n}_{{1j}}},{{n}_{{2j}}})$ – вектор единичной внешней нормали ${{{\mathbf{n}}}_{{\mathbf{j}}}}$ к границе $j$-й ячейки. При отсутствии множителя ${{h}_{1}}{{h}_{2}}$ в (2.7) ${{k}_{{{{m}_{2}}}}},\;{{k}_{{{{m}_{3}}}}}$ нужно согласовывать со значениями шагов сетки (см. п. 2.3). При этом в оказавшихся внутри отверстия точках согласования н-ячеек, к которым не присоединялись н-ячейки, уравнения согласования (2.6), (2.7) не выписывались (фиг. 2, ячейки 12–15).

В p-МКНК1$(j = 1)$ на внешней границе ${{\gamma }_{0}}$ в ${{N}_{{{\text{nu}}}}}$ точках, приблизительно равномерно расположенных на ней (фиг. 3), выписываются краевые условия

(2.8)
${{k}_{{{{b}_{0}}}}}{{v}_{{hj}}} = {{k}_{{{{b}_{0}}}}}{{g}_{1}}({{x}_{1}},{{x}_{2}}),$
(2.9)
${{k}_{{{{b}_{1}}}}}\left( {\frac{{{{n}_{1}}}}{{{{h}_{1}}}}\frac{{\partial {{v}_{{hj}}}}}{{\partial {{y}_{1}}}} + \frac{{{{n}_{2}}}}{{{{h}_{2}}}}\frac{{\partial {{v}_{{hj}}}}}{{\partial {{y}_{2}}}}} \right) = {{k}_{{{{b}_{1}}}}}{{g}_{2}}({{x}_{1}},{{x}_{2}}),$
где $({{n}_{1}},{{n}_{2}})$ – вектор единичной внешней нормали ${\mathbf{n}}$ к границе ${{\gamma }_{0}}$, ${{k}_{{{{b}_{0}}}}},{{k}_{{{{b}_{1}}}}}$ – положительные весовые множители краевых условий Дирихле. Информацию о технике расстановки точек записи краевых условий на дискретно заданной границе области можно найти в [17]. Если нерегулярная область является многосвязной с внутренними контурами – окружностями радиусов ${{r}_{i}}$, $i = \overline {1,{{N}_{\gamma }}} $, то на каждом внутреннем контуре дополнительно равномерно располагается $[2{{r}_{i}}{{N}_{{{\text{nu}}}}}{\text{/}}\sqrt {{{d}_{1}}{{d}_{2}}} ]$ точек для записи краевых условий (фиг. 3).

Во всех вариантах МКНК (кроме p-МКНК$_{1}$) на внешней стороне граничной $j$-й ячейки, в том числе в объединенной граничной н-ячейке, $j = \overline {1,{{N}_{{{\text{cells}}}}}} $, которая может совпадать с некоторой частью одной или частями нескольких границ ${{\gamma }_{i}}$, $i = \overline {1,{{N}_{\gamma }}} $, области $\Omega $, в $(4 - {{N}_{s}})[{{N}_{{{\text{col}}}}}{\text{/}}4]$ точках на ней (фиг. 2), выписываются краевые условия (2.8), (2.9), где ${{N}_{s}} = \overline {1,3} $ – количество ячеек, соседних с рассматриваемой. При этом, если к граничной н-ячейке не присоединялась какая-либо н-ячейка, то вместо $(4 - {{N}_{s}})[{{N}_{{{\text{col}}}}}{\text{/}}4]$ точек записи краевых условий на ее границе равномерно располагается столько точек записи краевых условий, сколько точек коллокации и согласования располагается внутри отверстия (фиг. 2, ячейки 12–15). Если отверстие оказалось полностью внутри некоторой материнской ячейки (фиг. 2, ячейка 16), то на границе данного отверстия располагается ${{N}_{{{\text{col}}}}} + [{{N}_{{{\text{col}}}}}{\text{/}}4] - {{N}_{{{\text{cic}}}}}$ точек записи краевых условий, где ${{N}_{{{\text{cic}}}}}$ – количество точек коллокации, которые лежат в материнской ячейке, но вне данного отверстия.

Аналогично работам [16], [17], здесь в целях достижения аппроксимации повышенной точности дифференциальной задачи применялись “законтурные” точки коллокации и согласования в граничных н-ячейках около внешней дискретно заданной границы ${{\gamma }_{0}}$ (фиг. 2). У н-ячеек, которые являются самостоятельными, имеют четыре соседа и которые пересечены границами ${{\gamma }_{i}}$, $i = \overline {1,{{N}_{\gamma }}} $, круглых отверстий, “законтурные” части в произвольном случае зачастую малы по сравнению с “законтурными” частями граничных н-ячеек, которые пересечены внешней границей ${{\gamma }_{0}}$. Здесь предложено не использовать у таких ячеек их законтурные части для записи уравнений приближенной задачи. Однако зачастую длина границ ${{\gamma }_{i}}$, $i = \overline {1,{{N}_{\gamma }}} $, существенное меньше длины внешней границы ${{\gamma }_{0}}$. В таких случаях требование на количество краевых условий становится жестче: необходимо по возможности выписывать краевые условия в каждой ячейке, которая содержит границу ${{\gamma }_{i}}$, $i = \overline {1,{{N}_{\gamma }}} $, даже если ячейка имеет четыре соседа, т.е. ${{N}_{s}} = 4$.

Объединяя уравнения приближенной задачи (2.5)–(2.9) в соответствии с конкретным вариантом МКНК и типом ячеек, относительно искомых ${{c}_{{{{i}_{1}}{{i}_{2}},j}}}$ имеем переопределенную локальную СЛАУ $Ax = b$, где $A$ – прямоугольная вещественная матрица, $x$, $b$ – векторы неизвестных и правых частей уравнений соответственно. Глобальная СЛАУ здесь решается с помощью метода итераций по подобластям (см. [16]). При построении решения в каждой ячейке делается $QR$-декомпозиция матрицы $A$. Итерационный процесс продолжается до тех пор, пока не выполнится условие

(2.10)
$\mathop {\max }\limits_{{{i}_{1}}{{i}_{2}}j} \left| {c_{{{{i}_{1}}{{i}_{2}},j}}^{{n + 1}} - c_{{{{i}_{1}}{{i}_{2}},j}}^{n}} \right| < \epsilon ,$
где $c_{{{{i}_{1}}{{i}_{2}},j}}^{n}$ – коэффициент полинома в $j$-й ячейке на $n$-й итерации, $j = \overline {1,{{N}_{{{\text{cells}}}}}} $, $\epsilon $ – наперед заданная малая константа, называемая псевдопогрешностью решения.

2.3. Ускорение итерационного процесса

Для значительного уменьшения количества итераций и времени, необходимого для проведения расчета, здесь будем использовать: комбинированное применение (см. [24]) операции продолжения на многосеточном комплексе в методе Федоренко (см. [26]); диагонального предобуславливателя (см. [27]); на всех сетках комплекса метода ускорения сходимости итерационного процесса, основанного на подпространствах Крылова (см. [28]); свойства матриц локальных СЛАУ в линейной задаче с постоянными коэффициентами в МКНК (см. [22]). Значения весовых множителей в МКНК влияют на точность приближенного решения задачи, скорость сходимости итерационного процесса, обусловленность локальных СЛАУ (см. [21]). Поэтому в нем одним из эффективных способов уменьшения времени решения задачи является выбор оптимальных значений этих множителей. Авторы представленной работы воспользовались опытом и результатами предыдущих работ, посвященных решению бигармонического уравнения (см. [17], [22]). Однако в некоторых примерах можно заметить различающиеся значения весовых множителей. Их поиск был осуществлен по алгоритмам из [21], [22].

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

Пусть у нас имеется решение, построенное на сетке размера ${{N}_{1}}{\text{/}}2 \times {{N}_{2}}{\text{/}}2$ с ячейками размера $4{{h}_{1}} \times 4{{h}_{2}}$, которое в $j$-й ячейке имеет вид

${{v}_{{hj}}}({{Y}_{1}},{{Y}_{2}}) = \sum\limits_{{{i}_{1}} = 0}^{{{K}_{1}}} \,\sum\limits_{{{i}_{2}} = 0}^{{{K}_{2}}} \,{{c}_{{{{i}_{1}}{{i}_{2}},j}}}{{\phi }_{{{{i}_{1}}}}}({{Y}_{1}}){{\phi }_{{{{i}_{2}}}}}({{Y}_{2}}),$
где ${{\phi }_{{{{i}_{1}}}}}({{Y}_{1}}){{\phi }_{{{{i}_{2}}}}}({{Y}_{2}})$ – базисный элемент пространства, в котором ищется решение.

Каждая материнская ячейка сетки размера ${{N}_{1}}{\text{/}}2 \times {{N}_{2}}{\text{/}}2$ (грубая сетка) содержит в себе четыре материнские ячейки сетки размера ${{N}_{1}} \times {{N}_{2}}$ (мелкая сетка). Пусть для любой точки на плоскости: $({{Y}_{1}},{{Y}_{2}})$ – ее координаты в локальной системе координат, связанной с $j$-й ячейкой грубой сетки, $({{y}_{1}},{{y}_{2}})$ – в локальной системе координат, связанной с одной из четырех ячеек мелкой сетки, содержащихся в $j$-й ячейке грубой сетки, $({{x}_{1}},{{x}_{2}})$ – в глобальной системе координат. Обозначим $({{X}_{{1j}}},{{X}_{{2j}}})$ – центр $j$-й ячейки грубой сетки, $({{x}_{{1j}}},{{x}_{{2j}}})$ – центр одной из четырех ячеек, содержащихся в $j$-й ячейке грубой сетки. Тогда из формул (2.1) следует:

${{x}_{1}} = 2{{h}_{1}}{{Y}_{1}} + {{X}_{{1j}}},\quad {{x}_{2}} = 2{{h}_{2}}{{Y}_{2}} + {{X}_{{2j}}},$
${{y}_{1}} = (2{{h}_{1}}{{Y}_{1}} + {{X}_{{1j}}} - {{x}_{{1j}}}){\text{/}}{{h}_{1}},\quad {{y}_{2}} = (2{{h}_{2}}{{Y}_{2}} + {{X}_{{2j}}} - {{x}_{{2j}}}){\text{/}}{{h}_{2}}.$
Отсюда получаем
${{Y}_{1}}({{y}_{1}}) = 0.5{{y}_{1}} + d{{y}_{1}},\quad {{Y}_{2}}({{y}_{2}}) = 0.5{{y}_{2}} + d{{y}_{2}},$
где $d{{y}_{1}} = ({{x}_{{1j}}} - {{X}_{{1j}}}){\text{/}}2{{h}_{1}}$, $d{{y}_{2}} = ({{x}_{{1j}}} - {{X}_{{1j}}}){\text{/}}2{{h}_{2}}$.

Запишем решение в ячейке мелкой сетки в виде

(2.11)
${{v}_{{hj}}}({{y}_{1}},{{y}_{2}}) = \sum\limits_{{{i}_{1}} = 0}^{{{K}_{1}}} \,\sum\limits_{{{i}_{2}} = 0}^{{{K}_{2}}} \,{{\tilde {c}}_{{{{i}_{1}}{{i}_{2}}}}}{{\phi }_{{{{i}_{1}}}}}({{y}_{1}}){{\phi }_{{{{i}_{2}}}}}({{y}_{2}}) = \sum\limits_{{{i}_{1}} = 0}^{{{K}_{1}}} \,\sum\limits_{{{i}_{2}} = 0}^{{{K}_{2}}} \,{{c}_{{{{i}_{1}}{{i}_{2}},j}}}{{\phi }_{{{{i}_{1}}}}}(0.5{{y}_{1}} + d{{y}_{1}}){{\phi }_{{{{i}_{2}}}}}(0.5{{y}_{2}} + d{{y}_{2}}).$

Базис из полиномов Чебышёва. Найдем коэффициенты ${{\tilde {c}}_{{{{i}_{1}}{{i}_{2}}}}}$ в случае, когда ${{\phi }_{{{{i}_{1}}}}},{{\phi }_{{{{i}_{2}}}}}$ – полиномы Чебышёва степени ${{i}_{1}}$ и ${{i}_{2}}$ соответственно.

Для удобства введем матрицы ${{B}_{1}} = (b_{{i,l}}^{1})_{{i = 0,l = 0}}^{{{{K}_{1}},{{K}_{1}}}}$ и ${{B}_{2}} = (b_{{i,l}}^{2})_{{i = 0,l = 0}}^{{{{K}_{2}},{{K}_{2}}}}$ такие, что

(2.12)
${{\phi }_{{{{i}_{1}}}}}(a{{y}_{1}} + b) = \sum\limits_{l = 0}^{{{i}_{1}}} \,b_{{i,l}}^{1}{{\phi }_{l}}({{y}_{1}}),\quad {{\phi }_{{{{i}_{2}}}}}(c{{y}_{2}} + d) = \sum\limits_{l = 0}^{{{i}_{2}}} \,b_{{i,l}}^{2}{{\phi }_{l}}({{y}_{2}}),$
где $a,b,c,d$ – произвольные константы. Выпишем матрицу ${{B}_{1}}$ следующим образом. Для ${{i}_{1}} = 0$ и ${{i}_{1}} = 1$ имеем
${{\phi }_{0}}(a{{y}_{1}} + b) = 1 = {{\phi }_{0}}({{y}_{1}}),$
${{\phi }_{1}}(a{{y}_{1}} + b) = a{{y}_{1}} + b = a{{\phi }_{1}}({{y}_{1}}) + b{{\phi }_{0}}({{y}_{1}}).$
Для ${{i}_{1}} = 2$
${{\phi }_{2}}(a{{y}_{1}} + b) = 2(a{{y}_{1}} + b{{)}^{2}} - 1 = {{a}^{2}}{{\phi }_{2}}({{y}_{1}}) + 4ab{{\phi }_{1}}({{y}_{1}}) + (2{{b}^{2}} - 1 + {{a}^{2}}){{\phi }_{0}}({{y}_{1}}).$
Следовательно, матрица ${{B}_{1}}$ имеет вид
${{B}_{1}} = \left( {\begin{array}{*{20}{c}} 1&0&0&0& \ldots \\ b&a&0&0& \ldots \\ {2{{b}^{2}} - 1 + {{a}^{2}}}&{4ab}&{{{a}^{2}}}&0& \ldots \\ \vdots & \vdots & \vdots & \vdots & \vdots \end{array}} \right),$
а матрица ${{B}_{2}}$ получается аналогичным способом. Заполняя $i$-ю строку матрицы ${{B}_{1}}$ ($i \geqslant 2$) и, используя рекуррентное соотношение для полиномов Чебышёва, будем иметь
${{\phi }_{{i + 1}}}(a{{y}_{1}} + b) = 2(a{{y}_{1}} + b){{\phi }_{i}}(a{{y}_{1}} + b) - {{\phi }_{{i - 1}}}(a{{y}_{1}} + b).$
Подставим в правую часть полученного соотношения выражения для ${{\phi }_{i}}$ и ${{\phi }_{{i - 1}}}$ из (2.12) и проведем тождественное преобразование

$\begin{gathered} {{\phi }_{{i + 1}}}(a{{y}_{1}} + b) = 2(a{{y}_{1}} + b)\sum\limits_{l = 0}^i \,b_{{i,l}}^{1}{{\phi }_{l}}({{y}_{1}}) - \sum\limits_{l = 0}^{i - 1} \,b_{{i - 1,l}}^{1}{{\phi }_{l}}({{y}_{1}}) = a\sum\limits_{l = 1}^i \,b_{{i,l}}^{1}(2{{y}_{1}}{{\phi }_{l}}({{y}_{1}}) - {{\phi }_{{l - 1}}}({{y}_{1}})) + \\ \, + a\sum\limits_{l = 1}^i \,b_{{i,l}}^{1}{{\phi }_{{l - 1}}}({{y}_{1}}) + 2a{{y}_{1}}b_{{i,0}}^{1}{{\phi }_{0}}({{y}_{1}}) + 2b\sum\limits_{l = 0}^i \,b_{{i,l}}^{1}{{\phi }_{l}}({{y}_{1}}) - \sum\limits_{l = 0}^{i - 1} \,b_{{i - 1,l}}^{1}{{\phi }_{l}}({{y}_{1}}). \\ \end{gathered} $

Используя рекуррентное соотношение для полиномов Чебышёва и тождество ${{y}_{1}}{{\phi }_{0}}({{y}_{1}}) = {{\phi }_{1}}({{y}_{1}})$, перепишем последнее соотношение в виде

${{\phi }_{{i + 1}}}(a{{y}_{1}} + b) = a\sum\limits_{l = 1}^i \,b_{{i,l}}^{1}{{\phi }_{{l + 1}}}({{y}_{1}}) + a\sum\limits_{l = 1}^i \,b_{{i,l}}^{1}{{\phi }_{{l - 1}}}({{y}_{1}}) + 2ab_{{i,0}}^{1}{{\phi }_{1}}({{y}_{1}}) + 2b\sum\limits_{l = 0}^i \,b_{{i,l}}^{1}{{\phi }_{l}}({{y}_{1}}) - \sum\limits_{l = 0}^{i - 1} \,b_{{i - 1,l}}^{1}{{\phi }_{l}}({{y}_{1}}).$
Поменяем пределы суммирования у первых двух слагаемых:
${{\phi }_{{i + 1}}}(a{{y}_{1}} + b) = a\sum\limits_{l = 2}^{i + 1} \,b_{{i,l - 1}}^{1}{{\phi }_{l}}({{y}_{1}}) + a\sum\limits_{l = 0}^{i - 1} \,b_{{i,l + 1}}^{1}{{\phi }_{l}}({{y}_{1}}) + 2ab_{{i,0}}^{1}{{\phi }_{1}}({{y}_{1}}) + 2b\sum\limits_{l = 0}^i \,b_{{i,l}}^{1}{{\phi }_{l}}({{y}_{1}}) - \sum\limits_{l = 0}^{i - 1} \,b_{{i - 1,l}}^{1}{{\phi }_{l}}({{y}_{1}}).$
Сгруппировав слагаемые в полученном соотношении, для элементов $b_{{i + 1,l}}^{1}$ матрицы ${{B}_{1}}$ получим рекуррентные соотношения

$\begin{gathered} b_{{i + 1,0}}^{1} = ab_{{i,1}}^{1} + 2bb_{{i,0}}^{1} - b_{{i - 1,0}}^{1}, \\ b_{{i + 1,1}}^{1} = ab_{{i,2}}^{1} + 2ab_{{i,0}}^{1} + 2bb_{{i,1}}^{1} - b_{{i - 1,1}}^{1}, \\ b_{{i + 1,l}}^{1} = ab_{{i,l - 1}}^{1} + ab_{{i,l + 1}}^{1} + 2bb_{{i,l}}^{1} - b_{{i - 1,l}}^{1},\quad l = \overline {2,i - 1} , \\ b_{{i + 1,i}}^{1} = ab_{{i,i - 1}}^{1} + 2bb_{{i,i}}^{1}, \\ b_{{i + 1,i + 1}}^{1} = ab_{{i,i}}^{1}. \\ \end{gathered} $

Воспользуемся равенством (2.11) в представлении решения (2.3), формально заменив ${{K}_{1}}$ на ${{K}_{1}} - 1$, ${{K}_{2}}$ на ${{K}_{2}} - 1$. В формуле (2.12) положим $a = c = 0.5,$ $b = d{{y}_{1}},$ $d = d{{y}_{2}}$ и сконструируем матрицы ${{B}_{1}},{{B}_{2}}$ по полученным формулам. В результате будем иметь

$\begin{gathered} {{v}_{{hj}}}({{y}_{1}},{{y}_{2}}) = \sum\limits_{{{i}_{1}} = 0}^{{{K}_{1}} - 1} \,\sum\limits_{{{i}_{2}} = 0}^{{{K}_{2}} - 1} \,{{c}_{{{{i}_{1}}{{i}_{2}},j}}}{{\phi }_{{{{i}_{1}}}}}(0.5{{y}_{1}} + d{{y}_{1}}){{\phi }_{{{{i}_{2}}}}}(0.5{{y}_{2}} + d{{y}_{2}}) = \\ = \sum\limits_{{{i}_{1}} = 0}^{{{K}_{1}} - 1} \,\sum\limits_{{{i}_{2}} = 0}^{{{K}_{2}} - 1} \,{{c}_{{{{i}_{1}}{{i}_{2}},j}}}\left( {\sum\limits_{{{k}_{1}} = 0}^{{{i}_{1}}} \,b_{{{{i}_{1}},{{k}_{1}}}}^{1}{{\phi }_{{{{k}_{1}}}}}({{y}_{1}})} \right)\left( {\sum\limits_{{{k}_{2}} = 0}^{{{i}_{2}}} b_{{{{i}_{2}},{{k}_{2}}}}^{2}{{\phi }_{{{{k}_{2}}}}}({{y}_{2}})} \right) = \sum\limits_{{{i}_{1}} = 0}^{{{K}_{1}} - 1} \,\sum\limits_{{{i}_{2}} = 0}^{{{K}_{2}} - 1} \left[ {\sum\limits_{{{k}_{1}} = {{i}_{1}}}^{{{K}_{1}} - 1} \,\sum\limits_{{{k}_{2}} = {{i}_{2}}}^{{{K}_{2}} - 1} \,{{c}_{{{{k}_{1}}{{k}_{2}},j}}}b_{{{{k}_{1}},{{i}_{1}}}}^{1}b_{{{{k}_{2}},{{i}_{2}}}}^{2}} \right]{{\phi }_{{{{i}_{1}}}}}({{y}_{1}}){{\phi }_{{{{i}_{2}}}}}({{y}_{2}}). \\ \end{gathered} $
Следовательно, коэффициенты ${{\tilde {c}}_{{{{i}_{1}}{{i}_{2}}}}} = \sum\nolimits_{{{k}_{1}} = {{i}_{1}}}^{{{K}_{1}} - 1} {\sum\nolimits_{{{k}_{2}} = {{i}_{2}}}^{{{K}_{2}} - 1} {{{c}_{{{{k}_{1}}{{k}_{2}},j}}}b_{{{{k}_{1}},{{i}_{1}}}}^{1}b_{{{{k}_{2}},{{i}_{2}}}}^{2}} } $.

Если ${{K}_{1}} - 1 = {{K}_{2}} - 1 = K - 1$, то можно показать, что для осуществления операции продолжения необходимо $O({{N}_{{{\text{cells}}}}}{{K}^{4}})$ арифметических операций, в то время как для выполнения одной итерации требуется $O({{N}_{{{\text{cells}}}}}{{K}^{6}})$ арифметических операций.

Базис из мономов. Воспользуемся соотношением (2.11) в представлении решения (2.4), положив ${{K}_{1}} = K,$ ${{K}_{2}} = K - {{i}_{1}}$. Применяя формулу бинома Ньютона, получаем

$\begin{gathered} {{v}_{{hj}}}({{y}_{1}},{{y}_{2}}) = \sum\limits_{{{i}_{1}} = 0}^K \,\sum\limits_{{{i}_{2}} = 0}^{K - {{i}_{1}}} \,{{c}_{{{{i}_{1}}{{i}_{2}},j}}}{{\phi }_{{{{i}_{1}}}}}(0.5{{y}_{1}} + d{{y}_{1}}){{\phi }_{{{{i}_{2}}}}}(0.5{{y}_{2}} + d{{y}_{2}}) = \\ = \sum\limits_{{{i}_{1}} = 0}^K \,\sum\limits_{{{i}_{2}} = 0}^{K - {{i}_{1}}} {{c}_{{{{i}_{1}}{{i}_{2}},j}}}\left( {\sum\limits_{{{k}_{1}} = 0}^{{{i}_{1}}} \,C_{{{{i}_{1}}}}^{{{{k}_{1}}}}d{{y}_{1}}^{{{{i}_{1}} - {{k}_{1}}}}{{{0.5}}^{{{{k}_{1}}}}}{{\phi }_{{{{k}_{1}}}}}({{y}_{1}})} \right)\left( {\sum\limits_{{{k}_{2}} = 0}^{{{i}_{2}}} \,C_{{{{i}_{2}}}}^{{{{k}_{2}}}}d{{y}_{2}}^{{{{i}_{2}} - {{k}_{2}}}}{{{0.5}}^{{{{k}_{2}}}}}{{\phi }_{{{{k}_{2}}}}}({{y}_{2}})} \right) = \\ = \sum\limits_{{{i}_{1}} = 0}^K \,\sum\limits_{{{i}_{2}} = 0}^{K - {{i}_{1}}} \left[ {\sum\limits_{{{k}_{1}} = {{i}_{1}}}^K \,\sum\limits_{{{k}_{2}} = {{i}_{2}}}^{K - {{k}_{1}}} \,{{c}_{{{{k}_{1}}{{k}_{2}},j}}}C_{{{{k}_{1}}}}^{{{{i}_{1}}}}d{{y}_{1}}^{{{{k}_{1}} - {{i}_{1}}}}{{{0.5}}^{{{{i}_{1}}}}}C_{{{{k}_{2}}}}^{{{{i}_{2}}}}d{{y}_{2}}^{{{{k}_{2}} - {{i}_{2}}}}{{{0.5}}^{{{{i}_{2}}}}}} \right]{{\phi }_{{{{i}_{1}}}}}({{y}_{1}}){{\phi }_{{{{i}_{2}}}}}({{y}_{2}}), \\ \end{gathered} $
где $C_{n}^{k} = \frac{{n!}}{{k!(n - k)!}}$. Следовательно, ${{\tilde {c}}_{{{{i}_{1}}{{i}_{2}}}}} = \sum\nolimits_{{{k}_{1}} = {{i}_{1}}}^K {\sum\nolimits_{{{k}_{2}} = {{i}_{2}}}^{K - {{k}_{1}}} {{{c}_{{{{k}_{1}}{{k}_{2}},j}}}C_{{{{k}_{1}}}}^{{{{i}_{1}}}}d{{y}_{1}}^{{{{k}_{1}} - {{i}_{1}}}}{{{0.5}}^{{{{i}_{1}}}}}C_{{{{k}_{2}}}}^{{{{i}_{2}}}}d{{y}_{2}}^{{{{k}_{2}} - {{i}_{2}}}}{{{0.5}}^{{{{i}_{2}}}}}} } $.

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

3. РЕЗУЛЬТАТЫ ЧИСЛЕННЫХ ЭКСПЕРИМЕНТОВ

Рассмотрим задачу (1.1)–(1.4) и другие краевые задачи для бигармонического уравнения. В представленных ниже таблицах приведены значения относительной погрешности приближенного решения в дискретном аналоге нормы ${{L}_{2}}$:

(3.1)
${{\left\| {{{E}_{r}}} \right\|}_{2}} = \sqrt {\frac{{\sum\limits_{j = 1}^{{{N}_{{{\text{cells}}}}}} \sum\limits_{m = 1}^{{{Q}_{j}}} {{{(u({{x}_{{1m}}},{{x}_{{2m}}}) - {{u}_{{hj}}}({{x}_{{1m}}},{{x}_{{2m}}}))}}^{2}}}}{{\sum\limits_{j = 1}^{{{N}_{{{\text{cells}}}}}} \sum\limits_{m = 1}^{{{Q}_{j}}} {{{(u({{x}_{{1m}}},{{x}_{{2m}}}))}}^{2}}}}} ,$
относительной погрешности приближенного решения в бесконечной норме
(3.2)
${{\left\| {{{E}_{r}}} \right\|}_{\infty }} = \frac{{\mathop {\max }\limits_{j = \overline {1,{{N}_{{{\text{cells}}}}}} } (\mathop {\max }\limits_{m = \overline {1,{{Q}_{j}}} } \left| {u({{x}_{{1m}}},{{x}_{{2m}}}) - {{u}_{{hj}}}({{x}_{{1m}}},{{x}_{{2m}}})} \right|)}}{{\mathop {\max }\limits_{j = \overline {1,{{N}_{{{\text{cells}}}}}} } (\mathop {\max }\limits_{m = \overline {1,{{Q}_{j}}} } \left| {u({{x}_{{1m}}},{{x}_{{2m}}})} \right|)}},$
и абсолютной погрешности приближенного решения в бесконечной норме
(3.3)
${{\left\| {{{E}_{a}}} \right\|}_{\infty }} = \mathop {\max }\limits_{j = \overline {1,{{N}_{{{\text{cells}}}}}} } (\mathop {\max }\limits_{m = \overline {1,{{Q}_{j}}} } \left| {u({{x}_{{1m}}},{{x}_{{2m}}}) - {{u}_{{hj}}}({{x}_{{1m}}},{{x}_{{2m}}})} \right|),$
где ${{Q}_{j}}$ – количество равномерно распределенных контрольных точек $({{x}_{{1m}}},{{x}_{{2m}}})$, взятых в $j$-й материнской ячейке для подсчета в них погрешности, $u$ – точное решение задачи, ${{u}_{{hj}}}$ – приближенное решение в $j$-й ячейке. Здесь ${{Q}_{j}} = 10\,000$ в p-МКНК1, а в остальных вариантах МКНК ${{Q}_{j}} = 100$. Для сравнения с МКЭ (см. далее (3.8)) в табл. 11 и 12 приведены значения абсолютной погрешности в норме ${{L}_{2}}$, где интеграл считался с помощью прямого произведения одномерных квадратурных формул Гаусса c 20 узлами по каждому из направлений

${{\left\| {{{E}_{a}}} \right\|}_{{{{L}_{2}}}}} = \sqrt {\int\limits_\Omega {{{(u - {{u}_{h}})}}^{2}}d\Omega } .$

Порядок сходимости погрешности численного решения в случае h-МКНК определим по формуле

(3.4)
$R = \mathop {\log }\nolimits_2 \frac{{{{E}_{p}}}}{{{{E}_{c}}}},$
где ${{E}_{c}}$ – значение погрешности (${{\left\| {{{E}_{r}}} \right\|}_{2}}$, ${{\left\| {{{E}_{r}}} \right\|}_{\infty }}$ или ${{\left\| {{{E}_{a}}} \right\|}_{\infty }}$) на сетке размера ${{N}_{1}}$ $ \times $ ${{N}_{2}}$, ${{E}_{p}}$ – значение погрешности на сетке размера ${{N}_{1}}{\text{/}}2$ $ \times $ ${{N}_{2}}{\text{/}}2$. Обозначим через ${{N}_{{{\text{iter}}}}}$ – количество итераций, при котором останавливался итерационный процесс (2.10).

Во всех численных экспериментах во всех вариантах МКНК (кроме p-МКНК1) в итерационном процессе на самой грубой сетке здесь в начальном приближении решения взяты ${{c}_{{{{i}_{1}}{{i}_{2}},j}}} = 0.4$.

В данной работе во всех рассмотренных примерах известно точное решение задач. Однако в работах по МКНК показана возможность построения сходящегося решения с хорошей точностью и порядком сходимости зачастую не хуже второго при применении полиномов четвертой степени также и в случаях, когда решение различных краевых задач для бигармонического уравнения неизвестно (см. [16], [17]).

3.1. Задача Дирихле в односвязной звездной области

Рассмотрим задачу Дирихле (1.1)–(1.3) для бигармонического уравнения с тестовым решением $u({{x}_{1}},{{x}_{2}}) = x_{1}^{2} + x_{2}^{2} + {{e}^{{{{x}_{1}}}}}\cos ({{x}_{2}})$ в невыпуклой области (фиг. 3, без отверстия), внешняя граница ${{\gamma }_{0}}$ которой задана соотношениями

$\begin{array}{*{20}{l}} {{{x}_{1}} = (0.6 + 0.25\sin (5t))\cos (t),\quad t \in [0,2\pi ),} \\ {{{x}_{2}} = (0.6 + 0.25\sin (5t))\sin (t).} \end{array}$
В табл. 1 приведены результаты численных экспериментов в случае применения h-МКНК, где в качестве базисных элементов взяты мономы при $K = 4$. В качестве весовых множителей здесь и в p-МКНКh взяты ${{k}_{c}} = {{k}_{{{{m}_{0}}}}} = {{k}_{{{{m}_{2}}}}} = {{k}_{{{{b}_{0}}}}} = {{k}_{{{{b}_{1}}}}} = 1$, ${{k}_{{{{m}_{1}}}}} = 0.01,$ ${{k}_{{{{m}_{3}}}}} = 0.015$ из [17], $\epsilon {{ = 10}^{{ - 10}}}$ в h-МКНК, $\epsilon {{ = 10}^{{ - 14}}}$ в p-МКНКh. В табл. 2 приведены результаты численных экспериментов в случае применения p-МКНКh (4 ячейки) и p-МКНК$_{1}$ при использовании различных вариантов расстановки точек коллокации (фиг. 3). В первом столбце для указания способа расстановки точек коллокации в скобках приведены номера рисунков, на которых изображены соответствующие им расстановки.

Таблица 1.  

Результаты численных экспериментов (пример п. 3.1)

Chen и др., см. [3] h-МКНК, решение (2.4) при $K = 4$
${{N}_{1}} \times {{N}_{2}}$ ${{\left\| {{{E}_{a}}} \right\|}_{\infty }}$ $R$ ${{N}_{1}} \times {{N}_{2}}$ ${{\left\| {{{E}_{a}}} \right\|}_{\infty }}$ $R$ ${{N}_{{{\text{iter}}}}}$
64 × 64 8.15e–04 – 10 × 10 1.46e–05 – 109
128 × 128 2.08e–04 1.97 20 × 20 1.70e–06 3.10 121
256 × 256 8.07e–05 1.37 40 × 40 2.17e–07 2.97 205
512 × 512 1.83e–05 2.13 80 × 80 2.96e–08 2.87 361
Таблица 2.  

Значения погрешности ${{\left\| {{{E}_{a}}} \right\|}_{\infty }}$ для численных экспериментов (пример п. 3.1)

${{K}_{1}}$-1, ${{K}_{2}}$-1 или $K$ 5 7 9 11 13
Shao и др., см. [8] 1.74e–04 1.80e–06 2.69e–08 5.26e–11 3.09e–12
p-МКНК1 (фиг. 3а), решение (2.3) 1.04e–04 5.87e–07 2.62e–09 5.72e–12 1.77e–14
p-МКНК1 (фиг. 3а), решение (2.4) 4.54e–04 4.08e–06 3.42e–08 1.73e–10 5.16e–13
p-МКНК1 (фиг. 3б), решение (2.3) 1.02e–04 2.65e–07 7.04e–10 1.04e–12 1.15e–14
p-МКНК1 (фиг. 3б), решение (2.4) 4.67e–04 4.07e–06 3.16e–08 1.80e–10 5.20e–13
p-МКНК1 (фиг. 3в), решение (2.3) 9.99e–05 2.52e–07 8.63e–10 2.63e–12 1.42e–14
p-МКНК1 (фиг. 3в), решение (2.4) 4.72e–04 3.41e–06 2.45e–08 1.87e–10 3.25e–13
p-МКНКh, решение (2.4) 2.78e–04 1.47e–06 2.35e–09 5.44e–12 3.6e–14

В h-МКНК, p-МКНК1 и p-МКНКh дополнительно подсчитаны значения чисел обусловленности локальных и глобальных матриц. В первом случае без использования диагонального предобуславливателя в зависимости от размера сетки и типов ячеек значения порядка чисел обусловленности варьировались в пределах от ${{10}^{2}}$ до ${{10}^{3}}$ (см. табл. 3). Во втором случае с представлением решения в виде (2.4) порядок числа обусловленности менялся от ${{10}^{1}}$ до ${{10}^{4}}$ в зависимости от $K$, а при представлении решения в виде (7) – от ${{10}^{2}}$ до ${{10}^{8}}$ в зависимости от ${{K}_{1}} - 1$, ${{K}_{2}} - 1$ (см. табл. 4). Наконец, в последнем случае при увеличении степени $K$ максимальное значение числа обусловленности среди всех матриц локальных СЛАУ достигало порядка ${{10}^{7}}$ (см. табл. 4). Здесь число обусловленности вычислялось как отношение максимального сингулярного числа к минимальному, т.е. в спектральной норме. Отметим, что в [8] не приведены значения чисел обусловленности для этого примера, а только для случая прямоугольной области, где они были в диапазоне от $6.29 \times {{10}^{1}}$ до $6.51 \times {{10}^{2}}$ при изменении ${{K}_{1}} - 1$ и ${{K}_{2}} - 1$ в пределах от 7 до 13 соответственно.

Таблица 3.  

${{N}_{{{\text{cond}}}}}$ – количество ячеек в h-МКНК, в которых матрицы локальных СЛАУ имеют числа обусловленности порядка ${{10}^{2}}$ и ${{10}^{3}}$, condmax – максимальное число обусловленности среди всех матриц локальных СЛАУ и condin – число обусловленности для матриц СЛАУ внутренних ячеек

${{N}_{1}} \times {{N}_{2}}$ ${{N}_{{{\text{cond}}}}}$ condmax condin
cond порядка ${{10}^{2}}$ cond порядка ${{10}^{3}}$
10 × 10 50 0 3.34e+2 1.31e+2
20 × 20 192 0 6.99e+2 1.61e+2
40 × 40 746 30 1.67e+3 2.26e+2
80 × 80 2992 110 4.53e+3 3.91e+2
Таблица 4.  

Значения чисел обусловленности матриц в p-МКНК1 и их максимальные значения в p-МКНКh среди всех ячеек области

${{K}_{1}}$-1, ${{K}_{2}}$-1 или $K$ 5 7 9 11 13
Без использования диагонального предобуславливателя
p-МКНК1 (фиг. 3а), решение (2.3) 3.34e+2 8.79e+3 1.30e+5 1.23e+6 1.18e+7
p-МКНК1 (фиг. 3а), решение (2.4) 5.91e+1 2.73e+2 1.26e+3 6.93e+3 4.42e+4
p-МКНК1 (фиг. 3б), решение (2.3) 2.88e+2 1.60e+3 6.93e+3 6.71e+4 4.51e+5
p-МКНК1 (фиг. 3б), решение (2.4) 6.34e+1 2.75e+2 1.33e+3 6.58e+3 3.12e+4
p-МКНК1 (фиг. 3в), решение (2.3) 7.76e+2 1.20e+5 6.17e+5 2.81e+7 1.67e+8
p-МКНК1 (фиг. 3в), решение (2.4) 6.52e+1 3.58e+2 2.08e+3 1.14e+4 7.83e+4
p-МКНКh, решение (2.4) 1.12e+3 8.85e+3 9.85e+4 1.14e+6 1.23e+7
C использованием диагонального предобуславливателя
p-МКНК1 (фиг. 3а), решение (2.3) 4.81e+1 1.86e+2 8.54e+2 4.71e+3 2.51e+4
p-МКНК1 (фиг. 3а), решение (2.4) 3.51e+1 1.01e+2 3.44e+2 1.76e+3 1.23e+4
p-МКНК1 (фиг. 3б), решение (2.3) 8.60e+1 3.30e+2 9.10e+2 4.17e+3 1.90e+4
p-МКНК1 (фиг. 3б), решение (2.4) 3.95e+1 1.13e+2 6.98e+2 4.15e+3 1.88e+4
p-МКНК1 (фиг. 3в), решение (2.3) 2.60e+2 3.29e+3 1.01e+4 2.54e+6 1.13e+7
p-МКНК1 (фиг. 3в), решение (2.4) 4.09e+1 1.59e+2 7.30e+2 2.75e+3 1.70e+4

При дополнительном использовании диагонального предобуславливателя в p-МКНК1 число обусловленности здесь уменьшалось по-разному. В случае представления решения в виде (2.4) оно уменьшалось в среднем приблизительно в 2–4 раза, в виде (2.3) – на 1–3 десятичных порядка в зависимости от расположения точек коллокации. Например, при ${{K}_{1}} - 1 = {{K}_{2}} - 1 = 13$ и расстановкой точек как на фиг. 3a оно изменилось с 1.18e+7 на 2.51e+4, как на фиг. 3б – с 4.51e+5 на 1.90e+4, как на фиг. 3в – с 1.67e+8 на 1.13e+7. Увеличение числа обусловленности при применении базиса из полиномов Чебышёва прежде всего связано с тем, что в этом случае число неизвестных коэффициентов значительно больше, чем при использовании базиса из мономов при фиксированных значениях ${{K}_{1}} - 1 = {{K}_{2}} - 1 = K$, вследствие чего точек записи уравнений приближенной задачи берется больше по алгоритму, предложенному в данной работе. К тому же звездная область является невыпуклой, точки коллокации в случае, когда не используется законтурная часть для записи уравнений, будут расположены достаточно близко друг к другу внутри исходной области (а не в прямоугольнике). Это приводит к тому, что обусловленность матрицы СЛАУ ухудшается, хотя точность расчетов остается повышенной благодаря ортогональности используемых полиномов. В p-МКНК1 в расчетах использованы следующие значения множителей: ${{k}_{c}}{{ = 10}^{{ - 3}}}$, ${{k}_{{{{b}_{0}}}}} = {{k}_{{{{b}_{1}}}}} = 1$.

Из табл. 1 и 2 следует, что результаты, полученные применением МКНК, не только не уступают результатам, полученным другими авторами применением МКР повышенного порядка аппроксимации (см. [3]) и CTMMID (англ. Chebyshev tau meshless method based on the integration–differentiation) (см. [8]), но и превосходят их во многих случаях в этом примере. В [8] приближенное решение и его производные искались в виде (2.3). Из табл. 2 видно, что лучшие результаты достигнуты с применением p-МКНК$_{1}$.

3.2. Задача Дирихле в односвязной полигональной области

Рассмотрим в невыпуклой области (фиг. 4а) задачу Дирихле (1.1)–(1.3) для бигармонического уравнения с тестовым решением $u({{x}_{1}},{{x}_{2}}) = \sin (\pi {{x}_{1}})\sin (\pi {{x}_{2}})$. В [11] такая область разбивалась на пять подобластей и решение в каждой подобласти искалось отдельно в виде (2.3). В табл. 5 приведены результаты численных экспериментов, полученных с применением p-МКНК1 и p-МКНКh и CTMMID-DDM (англ. Chebyshev tau meshless method based on the integration-differentiation with domain decomposition method) (см. [11]). Поскольку в [11] представлены значения погрешности ${{\left\| {{{E}_{a}}} \right\|}_{\infty }}$ в пяти выбранных авторами некоторых подобластях (${{\Omega }_{1}},{{\Omega }_{2}},{{\Omega }_{3}},{{\Omega }_{4}},{{\Omega }_{5}}$), то для сравнения с полученными результатами и экономии места приведем в табл. 5 их осредненные значения ${{\left\| {{{E}_{a}}} \right\|}_{\infty }} = \left( {{{{\left\| {{{E}_{{a,{{\Omega }_{1}}}}}} \right\|}}_{\infty }} + {{{\left\| {{{E}_{{a,{{\Omega }_{2}}}}}} \right\|}}_{\infty }} + {{{\left\| {{{E}_{{a,{{\Omega }_{3}}}}}} \right\|}}_{\infty }} + {{{\left\| {{{E}_{{a,{{\Omega }_{4}}}}}} \right\|}}_{\infty }} + {{{\left\| {{{E}_{{a,{{\Omega }_{5}}}}}} \right\|}}_{\infty }}} \right){\text{/}}5$.

Фиг. 4.

Область решения задачи в примере 3.2 с координатами вершин многоугольника: $(1,1)$, $(0.5 - {\text{ctg}}70^\circ ,1)$, $(0.5,0)$, $( - 0.5,0)$, $( - 0.5 + {\text{ctg}}70^\circ ,0)$, $( - 1,1)$, $( - 1, - 1)$, $(1, - 1)$ (а). Область решения задачи в примере 3.3 с координатами вершин многоугольника: $(0,0)$, $(0, - 1)$, $(1, - 1)$, $(1,1)$, $( - 7{\text{/}}12,1)$, $( - 1,7{\text{/}}12)$, $( - 1,0)$. Круглое отверстие радиуса $1{\text{/}}3$ с центром в $(1{\text{/}}2, - 1{\text{/}}2)$, квадратное отверстие $[1{\text{/}}6,1{\text{/}}6]$$ \times $$[5{\text{/}}6,5{\text{/}}6]$ (б). Область решения задачи в примере 3.5 с координатами вершин: $(1,0)$, $(0,1)$, $(0,0.5)$, $(0.5,0)$ (в).

Из табл. 5 видно, что при применении p-МКНКh и p-МКНК1 удается получить решение, не уступающее по точности решению по сравнению со случаем применения CTMMID-DDM (см. [11]). Кроме того, при применении p-МКНКh удается достичь значение абсолютной погрешности ${{\left\| {{{E}_{a}}} \right\|}_{\infty }}$ порядка e–11. В p-МКНКh и p-МКНК1 при расчетах использовались значения множителей из предыдущего примера.

Таблица 5.  

Результаты численных экспериментов (пример п. 3.2)

Shao, Wu, см. [11]
${{K}_{1}}$-1, ${{K}_{2}}$-1 8 10 12 14
${{\left\| {{{E}_{a}}} \right\|}_{\infty }}$ 1.71e–03 4.34e–05 3.91e–06 2.41e–07
p-МКНК1 (фиг. 3в), решение (2.3)
${{K}_{1}}$-1, ${{K}_{2}}$-1 8 10 12 14
${{\left\| {{{E}_{a}}} \right\|}_{\infty }}$ 4.89e–03 1.31e–04 2.92e–06 4.03e–08
p-МКНКh, решение (2.4), ${{N}_{1}} \times {{N}_{2}} = 8 \times 8$
$K$ 5 7 9 11
${{\left\| {{{E}_{a}}} \right\|}_{\infty }}$ 6.34e–02 5.81e–05 2.02e–06 3.47e–09
p-МКНКh, решение (2.4), ${{N}_{1}} \times {{N}_{2}} = 16 \times 16$
$K$ 5 7 9 11
${{\left\| {{{E}_{a}}} \right\|}_{\infty }}$ 1.46e–02 2.61e–06 3.30e–08 3.88e–11

3.3. Задача Дирихле в многосвязной области

Рассмотрим в многосвязной области (фиг. 4б) задачу Дирихле (1.1)–(1.4) для бигармонического уравнения с тестовым решением $u({{x}_{1}},{{x}_{2}}) = \frac{1}{{{{\pi }^{4}}}}(1 + \cos (\pi {{x}_{1}}))(1 + \cos (\pi {{x}_{2}}))$. Здесь условие (1.4) выполнено почти всюду на границе квадратного отверстия. В табл. 6 приведены результаты численных экспериментов. В [10] исходная область разбивалась на три подобласти, в каждой из которых приближенное решение искалось отдельно. При этом в [10] не приведены значения погрешности при ${{K}_{1}},{{K}_{2}} > 14$. В [8] задача сразу решалась во всей области и там было достигнуто значение погрешности порядка ~e–13. Более того, при увеличении степени полиномов в [8] точность решения падала  при ${{K}_{1}} - 1$, ${{K}_{2}} - 1 > 20$.  Однако из табл. 6 видно,  что  с  помощью p-МКНК1 удается достигнуть погрешности порядка ~e–15. В данном примере использовался p-МКНК1 с представлением приближенного решения в виде (2.3) при расстановке точек коллокации, изображенной на фиг. 3а. В качестве весовых множителей были взяты ${{k}_{c}}{{ = 10}^{{ - 3}}}$, ${{k}_{{{{b}_{0}}}}} = {{k}_{{{{b}_{1}}}}} = 1$.

Таблица 6.  

Результаты численных экспериментов (пример п. 3.3)

Mai-Duy и др., см. [10]
${{K}_{1}}$-1, ${{K}_{2}}$-1 7 9 11 13
${{\left\| {{{E}_{r}}} \right\|}_{2}}$ 4.8e–05 2.5e–07 1.9e–09 8.7e–12
Shao и др., см. [8]
${{K}_{1}}$-1, ${{K}_{2}}$-1 8 12 16 20
${{\left\| {{{E}_{r}}} \right\|}_{2}}$ 8.06e–04 1.34e–05 2.58e–09 2.45e–13
p-МКНК$_{1}$ (фиг. 3а), решение (2.3)
${{K}_{1}}$-1, ${{K}_{2}}$-1 8 12 16 20
${{\left\| {{{E}_{r}}} \right\|}_{2}}$ 3.43e–05 2.71e–07 2.42e–11 5.13e–15

3.4. Расчет изгиба шарнирно-закрепленной прямоугольной пластины

Прогиб изотропной упругой пластины в рамках классической теории Кирхгофа–Лява (см. [29]) определяется из решения уравнения (1.1), где $w = u$ – прогиб срединной поверхности пластины, $f = q{\text{/}}D$, $q$ – интенсивность нагрузки, $D = E{{h}^{3}}{\text{/}}(12(1 - {{\nu }^{2}}))$ – жесткость пластины при цилиндрическом изгибе, $E,\nu $ – модуль Юнга и коэффициент Пуассона изотропного материала пластины, $h$ – толщина пластины. Рассмотрим следующие краевые условия закрепления, которые могут быть заданы на каждой части границ ${{\gamma }_{i}}$, $i = \overline {0,{{N}_{\gamma }}} $, пластины:

(3.5)
$w = 0,\;\;\frac{{\partial w}}{{\partial n}} = 0\; - \;{\text{защемленный}}\;{\text{край,}}$
(3.6)
$w = 0,\;\;{{M}_{n}} = 0\; - \;{\text{шарнирно - зекрепленный край}},$
(3.7)
${{M}_{n}} = 0,\;\;{{V}_{n}} = {{Q}_{n}} - \frac{{\partial {{M}_{{n\tau }}}}}{{\partial \tau }} = 0{\kern 1pt} \; - \;{\text{свободный край}},$
где ${\mathbf{n}}$ и $\tau $ – внешняя нормаль и касательный вектор к границам области ${{\gamma }_{i}}$, $i = \overline {0,{{N}_{\gamma }}} $, ${{M}_{n}}$ – изгибающий момент, ${{V}_{n}}$ – обобщенная поперечная сила, ${{Q}_{n}}$ – перерезывающая сила, ${{M}_{{n\tau }}}$ – крутящий момент (см. [29]). Здесь
${{M}_{n}} = {{M}_{{{{x}_{1}}}}}{{\cos }^{2}}\alpha + {{M}_{{{{x}_{2}}}}}{{\sin }^{2}}\alpha - 2{{M}_{{{{x}_{1}}{{x}_{2}}}}}\sin \alpha \cos \alpha ,$
${{M}_{{n\tau }}} = {{M}_{{{{x}_{1}}{{x}_{2}}}}}({{\cos }^{2}}\alpha - {{\sin }^{2}}\alpha ) + ({{M}_{{{{x}_{1}}}}} - {{M}_{{{{x}_{2}}}}})\sin \alpha \cos \alpha ,$
${{Q}_{n}} = {{Q}_{{{{x}_{1}}}}}\cos \alpha + {{Q}_{{{{x}_{2}}}}}\sin \alpha ,$
${{M}_{{{{x}_{1}}}}} = - D\left( {\frac{{{{\partial }^{2}}w}}{{\partial x_{1}^{2}}} + \nu \frac{{{{\partial }^{2}}w}}{{\partial x_{2}^{2}}}} \right),\quad {{M}_{{{{x}_{2}}}}} = - D\left( {\frac{{{{\partial }^{2}}w}}{{\partial x_{2}^{2}}} + \nu \frac{{{{\partial }^{2}}w}}{{\partial x_{1}^{2}}}} \right),\quad {{M}_{{{{x}_{1}}{{x}_{2}}}}} = D(1 - \nu )\frac{{{{\partial }^{2}}w}}{{\partial {{x}_{1}}\partial {{x}_{2}}}},$
${{Q}_{{{{x}_{1}}}}} = - D\frac{\partial }{{\partial {{x}_{1}}}}(\Delta w),\quad {{Q}_{{{{x}_{2}}}}} = - D\frac{\partial }{{\partial {{x}_{2}}}}(\Delta w),$
где $\alpha $ – угол между вектором нормали и осью ${{x}_{1}}$, отсчитываемый от ${{x}_{1}}$ против хода часовой стрелки, $\Delta $ – оператор Лапласа.

Пусть на шарнирно закрепленную (3.6) по всему контуру прямоугольную пластину размера ${{d}_{1}}$$ \times $${{d}_{2}}$ действует синусоидальная нагрузка $q = {{10}^{5}}\sin (\pi {{x}_{1}}{\text{/}}{{d}_{1}})\sin (\pi {{x}_{2}}{\text{/}}{{d}_{2}})$ [Па]. В этом случае существует точное аналитическое решение краевой задачи (1.1), (3.6) (см. [29])

(3.8)
$w({{x}_{1}},{{x}_{2}}) = \frac{{qd_{1}^{4}d_{2}^{4}}}{{{{\pi }^{4}}D{{{\left( {d_{1}^{2} + d_{2}^{2}} \right)}}^{2}}}}.$

В табл. 7 приведены результаты численных экспериментов, полученных h-МКНК с использованием операции продолжения на многосеточном комплексе. В четвертом столбце в скобках указано количество итераций при комбинированном применении операции продолжения и ускорения по Крылову. Здесь в расчетах взяты следующие значения параметров: ${{d}_{1}} = 10$ м, ${{d}_{2}} = 10$ м, $h = 0.1$ м, $E = 200$ ГПа, $\nu = 0.28$, ${{k}_{c}} = {{k}_{{{{m}_{0}}}}} = {{k}_{{{{m}_{2}}}}} = {{k}_{{{{b}_{0}}}}} = {{k}_{{{{b}_{2}}}}} = {{k}_{{{{m}_{1}}}}} = {{k}_{{{{m}_{3}}}}} = 1$, $\epsilon {{ = 10}^{{ - 14}}}$ в h-МКНК и p-МКНКh. Здесь ${{k}_{{{{b}_{2}}}}}$ – весовой множитель условия ${{M}_{n}} = {{g}_{2}}({{x}_{1}},{{x}_{2}})$ в переопределенной СЛАУ приближенной задачи, ${{d}_{{{\text{coeff}}}}}$ – максимальная разница между коэффициентами приближенного решения после первой итерации и коэффициентами начального приближения. Расположение точек коллокации в каждой ячейке соответствовало расположению точек, изображенных на фиг. 3б. Из табл. 7 следует, что применение операции продолжения позволяет сократить количество итераций приблизительно в два раза.

Таблица 7.  

Результаты численных экспериментов в прямоугольной области (пример п. 3.4)

${{N}_{1}}$ × ${{N}_{2}}$ Операция продолжения Без операции продолжения
${{\left\| {{{E}_{r}}} \right\|}_{\infty }}$ ${{d}_{{{\text{coeff}}}}}$ ${{N}_{{{\text{iter}}}}}$ ${{\left\| {{{E}_{r}}} \right\|}_{\infty }}$ ${{d}_{{{\text{coeff}}}}}$ ${{N}_{{{\text{iter}}}}}$
h-МКНК, реш. (2.3) при ${{K}_{1}} = {{K}_{2}} = 7$
8 × 8 4.76e–06 1.68e–06 90 (41) 4.76e–06 2.73e+3 166
16 × 16 1.11e–07 3.03e–08 142 (81) 1.11e–07 8.00e+3 329
32 × 32 1.91e–08 4.23e–10 218 (161) 1.91e–08 1.95e+3 682
h-МКНК, реш. (2.4) при $K = 6$
8 × 8 6.91e–05 1.06e–05 99 (41) 6.91e–05 1.55e+2 160
16 × 16 4.14e–06 2.27e–07 164 (58) 4.14e–06 9.42e+2 322
32 × 32 2.56e–07 8.18e–09 343 (164) 2.56e–07 5.93e+3 688

3.5. Краевая задача с граничными условиями, содержащими производные второго порядка

Рассмотрим в нерегулярной области (фиг. 1) краевую задачу для бигармонического уравнения с тестовым решением $u = \sin (\pi {{x}_{1}}{\text{/}}{{d}_{1}})\sin (\pi {{x}_{2}}{\text{/}}{{d}_{2}})$, где в численных экспериментах заданы на соответствующих границах значения функции $u$ и действующего на нее дифференциального оператора ${{M}_{n}} = {{g}_{{2i}}}({{x}_{1}},{{x}_{2}})$, $i = \overline {0,3} $, вычисленные из точного решения. Внешней границей ${{\gamma }_{0}}$ многосвязной области является плоская кривая “улитка Паскаля” (фиг. 1), заданная соотношениями

$\begin{array}{*{20}{l}} {{{x}_{1}} = 7 - 4\sin (t) + \frac{{43}}{{25}}\sin (2t),} \\ {{{x}_{2}} = 7.8 + 4\cos (t) - \frac{{43}}{{25}}\cos (2t),} \end{array}$
уравнения отверстий ${{({{x}_{1}} - 5.5)}^{2}} + {{({{x}_{2}} - 7)}^{2}} \leqslant {{0.6}^{2}}$, ${{({{x}_{1}} - 8.5)}^{2}} + {{({{x}_{2}} - 7.8)}^{2}} \leqslant {{0.5}^{2}}$, ${{({{x}_{1}} - 7)}^{2}} + ({{x}_{2}} - $ $ - \;{{4.5)}^{2}} \leqslant {{0.4}^{2}}$. В табл. 8 приведены результаты численных экспериментов, где в h-МКНК приближенное решение искалось в виде (2.3). Расположение точек коллокации в каждой ячейке соответствовало расположению точек, изображенных на фиг. 3а, т.е. в корнях ортогональных полиномов Чебышёва. В качестве точек согласования брались точки с координатами $( \pm 1,\alpha _{{{{i}_{2}}}}^{2})$, $(\alpha _{{{{i}_{1}}}}^{1}, \pm 1)$, ${{i}_{1}} = \overline {1,{{K}_{1}}} $, ${{i}_{2}} = \overline {1,{{K}_{2}}} $, $\alpha _{{{{i}_{1}}}}^{1}$ и $\alpha _{{{{i}_{2}}}}^{2}$ – корни многочленов Чебышёва I рода степеней ${{K}_{1}}$ и  ${{K}_{2}}$ соответственно. В расчетах взяты следующие значения параметров: ${{d}_{1}} = 10.6667$, ${{d}_{2}} = 9.64107$, ${{k}_{c}} = {{k}_{{{{m}_{0}}}}} = {{k}_{{{{m}_{2}}}}} = {{k}_{{{{b}_{0}}}}} = {{k}_{{{{b}_{2}}}}} = 1$, ${{k}_{{{{m}_{1}}}}} = 0.12,$ ${{k}_{{{{m}_{3}}}}} = 0.125$, $\epsilon {{ = 10}^{{ - 12}}}$ в h-МКНК и p-МКНК$_{h}$, ${{k}_{c}}{{ = 10}^{{ - 3}}}$, ${{k}_{{{{b}_{0}}}}} = {{k}_{{{{b}_{2}}}}} = 1$ в p-МКНК1.

Таблица 8.  

Результаты численных экспериментов в нерегулярной области (пример п. 3.5)

h-МКНК, решение (2.3) при ${{K}_{1}} = {{K}_{2}} = 5$ p-МКНКh, решение (2.3) p-МКНК1, решение (2.3)
${{N}_{1}}$ × ${{N}_{2}}$ ${{\left\| {{{E}_{r}}} \right\|}_{\infty }}$ $R$ ${{K}_{1}} = {{K}_{2}}$ ${{\left\| {{{E}_{r}}} \right\|}_{\infty }}$ ${{K}_{1}} = {{K}_{2}}$ ${{\left\| {{{E}_{r}}} \right\|}_{\infty }}$
12 × 12 3.92e–04 – 7 1.06e–05 13 7.88e–10
24× 24 5.02e–05 2.96 8 2.04e–06 15 4.98e–12
48 × 48 5.63e–06 3.15 9 1.85e–08 17 3.24e–13
96 × 96 1.45e–06 1.95 10 7.49e–09 19 5.73e–14

3.6. Краевая задача с различными граничными условиями

Рассмотрим в нерегулярной области (фиг. 4в) краевую задачу для бигармонического уравнения с тестовым решением $u = \sin (\pi {{x}_{1}}{\text{/}}{{d}_{1}})\sin (\pi {{x}_{2}}{\text{/}}{{d}_{2}})$, сформулированную в [8]. На дуге окружности заданы значения действующих на функцию $u$ дифференциальных операторов ${{M}_{n}}$ и ${{V}_{n}}$, на двух прилежащих к дуге прямолинейных краях заданы значения функции $u$ и действующего на нее дифференциального оператора ${{M}_{n}}$, на оставшемся прямолинейном крае заданы значения функции $u$ и $\frac{{\partial u}}{{\partial n}}$, вычисленные из точного решения.

В табл. 9 приведены результаты численного эксперимента в случае применения p-МКНК1 с представлением приближенного решения в виде (2.3) при расстановке точек коллокации по схеме, изображенной на фиг. 3б. В этой задаче в краевые условия входят производные второго и третьего порядков. В данной работе использованы следующие значения: ${{k}_{{{{b}_{0}}}}} = {{k}_{{{{b}_{1}}}}} = 1$, ${{k}_{c}}$, ${{k}_{{{{b}_{2}}}}}$ и ${{k}_{{{{b}_{3}}}}}$ варьировались в диапазоне 1e–3 ~ 1e–2 при фиксированной степени полинома. Здесь ${{k}_{{{{b}_{3}}}}}$ – весовой множитель перед условием ${{V}_{n}} = {{g}_{2}}({{x}_{1}},{{x}_{2}})$.

Таблица 9.  

Результаты численных экспериментов (пример 3.6)

Shao и др., см. [8]
${{K}_{1}}$-1, ${{K}_{2}}$-1 4 6 8 10 12
${{\left\| {{{E}_{r}}} \right\|}_{2}}$ 1.75e–02 1.47e–04 2.61e–06 1.29e–07 7.10e–10
p-МКНК1 (фиг. 3б), решение (2.3)
${{K}_{1}}$-1, ${{K}_{2}}$-1 4 6 8 10 12
${{\left\| {{{E}_{r}}} \right\|}_{2}}$ 2.80e–02 1.24e–04 3.09e–06 4.74e–08 2.99e–10

3.7. Задача с особенностью. Симметричный изгиб кольцевой пластины

Рассмотрим задачу об изгибе равномерной поперечной нагрузкой $q$ круглой пластины радиуса ${{r}_{1}}$, имеющей центральное круглое отверстие радиуса ${{r}_{0}}$. Поскольку нагрузка распределена симметрично относительно центра, то прогиб $w$ в полярных координатах имеет вид (см. [29])

(3.9)
$w = \frac{{q{{r}^{4}}}}{{64D}} + {{C}_{1}}{{r}^{2}}(\log r - 1) + {{C}_{2}}\frac{{{{r}^{2}}}}{4} + {{C}_{3}}\log r + {{C}_{4}},$
где $r = \sqrt {x_{1}^{2} + x_{2}^{2}} $ – полярный радиус, ${{C}_{1}}$,…,${{C}_{4}}$ – неизвестные константы, которые определяются из четырех граничных условий.

Рассмотрим сначала кольцевую пластину с защемленными (3.5) внешним и внутренним контурами. В этом случае граничные условия имеют вид

${{\left. w \right|}_{{r = {{r}_{1}}}}} = {{\left. {\frac{{\partial w}}{{\partial n}}} \right|}_{{r = {{r}_{1}}}}} = {{\left. w \right|}_{{r = {{r}_{0}}}}} = {{\left. {\frac{{\partial w}}{{\partial n}}} \right|}_{{r = {{r}_{0}}}}} = 0.$
Пусть $q{{ = 10}^{6}}$ Па, $\nu = 0.28$, $E = 200$ ГПа, $h = 0.1$ м, ${{r}_{1}} = 5$ м, ${{r}_{0}} = {{r}_{1}}{\text{/}}k$ м, $k$ – коэффициент отношения, ${{d}_{1}} = 10$ м, ${{d}_{2}} = 10$ м. Все значения весовых множителей были взяты равными единице. В табл. 10 приведены результаты численных экспериментов в случае применения h-МКНК при $K = 4$ и p-МКНК$_{h}$ на сетке 24 $ \times $ 24 с представлением приближенного решения в виде (2.4). Отметим, что в случае достаточно малого отверстия все точки коллокации и согласования располагаются вне отверстия (см. п. 2.2). Для учета граничных условий в такой ячейке на ее внешней стороне располагалась одна точка записи краевых условий.

Таблица 10.  

Результаты численных экспериментов (пример п. 3.7)

${{r}_{0}} = {{r}_{1}}{\text{/}}5$ ${{r}_{0}} = {{r}_{1}}{\text{/}}50$ ${{r}_{0}} = {{r}_{1}}{\text{/}}500$ ${{r}_{0}} = {{r}_{1}}{\text{/}}5000$
h-МКНК, $K = 4$
${{N}_{1}}$ × ${{N}_{2}}$ ${{\left\| {{{E}_{r}}} \right\|}_{\infty }}$ ${{N}_{1}}$ × ${{N}_{2}}$ ${{\left\| {{{E}_{r}}} \right\|}_{\infty }}$ ${{N}_{1}}$ × ${{N}_{2}}$ ${{\left\| {{{E}_{r}}} \right\|}_{\infty }}$ ${{N}_{1}}$ × ${{N}_{2}}$ ${{\left\| {{{E}_{r}}} \right\|}_{\infty }}$
8 × 8 3.60e–1 8 × 8 1.79e–1 8 × 8 2.78e–1 8 × 8 2.92e–1
16 × 16 6.40e–2 16 × 16 6.93e–2 16 × 16 7.60e–2 16 × 16 9.25e–2
32 × 32 5.20e–3 32 × 32 3.96e–2 32 × 32 1.72e–2 32 × 32 2.77e–2
64 × 64 1.23e–3 64 × 64 7.45e–3 64 × 64 4.24e–3 64 × 64 5.22e–3
p-МКНКh, ${{N}_{1}}$ × ${{N}_{2}}$ = 24 × 24, ${{r}_{0}} = {{r}_{1}}{\text{/}}5$
$K$ ${{\left\| {{{E}_{r}}} \right\|}_{\infty }}$ $K$ ${{\left\| {{{E}_{r}}} \right\|}_{\infty }}$ $K$ ${{\left\| {{{E}_{r}}} \right\|}_{\infty }}$ $K$ ${{\left\| {{{E}_{r}}} \right\|}_{\infty }}$
4 1.37e–2 6 2.50e–3 8 2.93e–4 10 2.72e–5
Таблица 11.  

Результаты численных экспериментов (пример п. 3.8)

$\alpha $ = 5/3
${{N}_{1}}$ × ${{N}_{2}}$ h-МКНК при $K = 4$ h-МКНК при $K = 6$   МКЭ [7] адапт. МКЭ [7]
DoF ${{\left\| {{{E}_{a}}} \right\|}_{{{{L}_{2}}}}}$ $R$ DoF ${{\left\| {{{E}_{a}}} \right\|}_{{{{L}_{2}}}}}$ $R$ DoF ${{\left\| {{{E}_{a}}} \right\|}_{{{{L}_{2}}}}}$ ${{\left\| {{{E}_{a}}} \right\|}_{{{{L}_{2}}}}}$
4 × 4 180 6.51e–3 – 336 2.82e–3 – 833 7.15e–3 8.88e–4
8 × 8 720 3.64e–3 0.84 1344 8.21e–4 1.78 3201 2.72e–3 2.16e–4
16 × 16 2880 9.93e–4 1.87 5376 1.67e–4 2.29 12 545 1.08e–3 5.34e–5
32 × 32 11 520 3.76e–4 1.40 21 504 1.17e–4 0.51 49 665 4.40e–4 1.32e–5
64 × 64 46 080 1.74e–4 1.11 86 016 7.52e–5 0.64 197 633 1.82e–4 3.28e–6
Таблица 12.  

Результаты численных экспериментов (пример п. 3.8) при различных $\alpha $

  h-МКНК при $K = 4$ h-МКНК при $K = 6$ h-МКНК при $K = 8$
${{N}_{1}}$ × ${{N}_{2}}$ DoF ${{\left\| {{{E}_{a}}} \right\|}_{{{{L}_{2}}}}}$ R DoF ${{\left\| {{{E}_{a}}} \right\|}_{{{{L}_{2}}}}}$ R DoF ${{\left\| {{{E}_{a}}} \right\|}_{{{{L}_{2}}}}}$ R
$\alpha $ = 7/3
4 × 4 180 6.55e–3 – 336 1.42e–3 – 720 1.78e–3 –
8 × 8 720 1.82e–3 1.85 1344 4.89e–4 1.53 2880 7.45e–4 1.26
16 × 16 2880 2.77e–4 2.72 5376 9.41e–5 2.37 11 520 1.91e–4 1.96
32 × 32 11 520 4.23e–5 2.71 21 504 1.52e–5 2.63 46 080 5.14e–5 1.89
64 × 64 46 080 6.83e–6 2.63 86 016 2.42e–6 2.65 184 320 1.43e–5 1.86
$\alpha $ = 10/3
4 × 4 180 7.49e–3 – 336 7.70e–4 – 720 2.61e–4 –
8 × 8 720 6.64e–4 3.50 1344 1.69e–4 2.19 2880 3.62e–5 2.85
16 × 16 2880 3.49e–5 4.25 5376 2.14e–5 2.98 11 520 5.33e–6 2.76
32 × 32 11 520 5.88e–7 5.89 21 504 2.39e–6 3.16 46 080 6.29e–7 3.08
64 × 64 46 080 4.07e–7 0.53 86 016 2.62e–7 3.19 184 320 1.04e–7 2.60
$\alpha $ = 13/3
4 × 4 180 7.87e–3 – 336 3.63e–4 – 720 7.28e–5 –
8 × 8 720 7.47e–4 3.40 1344 2.59e–5 3.81 2880 8.91e–6 3.03
16 × 16 2880 3.46e–5 4.43 5376 1.50e–6 4.11 11 520 5.89e–7 3.92
32 × 32 11 520 1.62e–6 4.42 21 504 8.05e–8 4.22 46 080 3.20e–8 4.18
64 × 64 46 080 1.23e–7 3.72 86 016 4.48e–9 4.17 184 320 1.84e–9 4.12
$\alpha $ = 16/3
4 × 4 180 5.45e–2 – 336 6.16e–4 – 720 7.13e–5 –
8 × 8 720 7.76e–3 2.81 1344 1.91e–5 5.01 2880 4.12e–6 4.11
16 × 16 2880 3.39e–4 4.51 5376 9.69e–7 4.30 11 520 1.11e–7 5.21
32 × 32 11 520 1.32e–5 4.68 21 504 5.76e–8 4.06 46 080 2.20e–9 5.65
64 × 64 46 080 6.54e–7 4.33 86 016 3.55e–9 4.02 184 320 3.21e–11 6.10

В предложенном здесь p-МКНК$_{1}$ приближенные решения (2.3) или (2.4) представляют собой равномерное приближение искомого решения в прямоугольнике, включающем всю исходную расчетную область и все отверстия в случае многосвязной области (в нашем примере $x_{1}^{2} + x_{2}^{2} \leqslant r_{0}^{2}$). Точное решение (3.9) данной задачи в точке $(0,0)$ имеет особенность в виде разрыва II рода из-за наличия в нем $\log r$. Поэтому решение (3.9) не является непрерывным в данном прямоугольнике и его не удается приблизить с помощью p-МКНК1 (2.3) или (2.4). Однако, зная вид особенности решения и асимптотику его поведения, можно модифицировать базис используемого линейного пространства в p-МКНК1. С другой стороны, как это видно из табл. 10, применение h-МКНК и p-МКНКh позволяет построить сходящееся приближенное решение в задаче с особенностью без изменения базиса применяемого пространства. Эта сходимость достигается благодаря использованию кусочно-полиномиального представления решения, что подтверждает универсальность и гибкость h-МКНК и p-МКНК$_{h}$ при решении краевых задач для УЧП в общем случае.

Для повышения точности аппроксимации и сходимости p-МКНК1 решения рассматриваемой задачи с особенностью необходимо добавить в его представление (2.3) или (2.4) слагаемые вида $\tilde {c}\log r$ и $\tilde {\tilde {c}}{{r}^{2}}\log r$, где $\tilde {c}$ и $\tilde {\tilde {c}}$ – неизвестные константы. В этом случае уже при ${{K}_{1}} = {{K}_{2}} = 5$ или $K = 4$ значение относительной погрешности ${{\left\| {{{E}_{r}}} \right\|}_{\infty }}$ достигается порядка e–13 при любых соотношениях ${{r}_{1}}{\text{/}}{{r}_{2}}$. На фиг. 5а показана форма прогиба кольцевой пластины с защемленными (3.5) внешним и внутренним контурами с ${{r}_{0}} = 1$ м.

Фиг. 5.

Величина прогиба $w$ кольцевой пластины с защемленным внешним контуром и защемленным (а) (свободным (б)) внутренним контуром.

Рассмотрим теперь кольцевую пластину c защемленным внешним (2.5) и со свободным (3.7) внутренним контурами. На фиг. 5б показана форма прогиба при ${{r}_{1}} = 5$ м, ${{r}_{0}} = 0.5$ м, $q{{ = 10}^{5}}$ Па. Значение ${{\left\| {{{E}_{r}}} \right\|}_{\infty }}$ здесь достигается также порядка e–13 при применении модифицированного p-МКНК$_{1}$. При применении МКНК с полиномиальным базисом при $K = 10$ на сетке $12 \times 12$ значение ${{\left\| {{{E}_{r}}} \right\|}_{\infty }}$ имеет величину порядка e–3.

3.8. Задача с особенностью в виде обращения в бесконечность производных решения

Рассмотрим задачу (1.1)–(1.3) в L-образной области $\Omega = $ $( - 1,1)$ $ \times $ $( - 1,1){{\backslash }}(0,1)$ $ \times $ $( - 1,0)$ с тестовым решением $u = {{r}^{\alpha }}\sin (\alpha \theta )$, где $(r,\theta )$ – полярные координаты, $\alpha > 0$ – const. Здесь $f = 0$, на границе заданы значения $u$ и $\partial u{\text{/}}\partial n$, вычисленные из точного решения.

В табл. 11 приведены результаты численного эксперимента с $\alpha = 53$. В этом случае вторые производные в точке $(0,0)$ обращаются в бесконечность. Известно, что решение уравнений эллиптического типа в области зависит от его значений на границе. Отсутствие гладкости даже в одной точке на ${{\gamma }_{0}}$ может сказываться на характере поведения решения в $\Omega $ (см. [30]). Здесь использовался МКНК (2.4) при расстановке точек коллокации по схеме, изображенной на фиг. 2. В табл. 11 DoF обозначает число степеней свободы, что в МКНК эквивалентно общему количеству неизвестных коэффициентов. Отметим, что в МКЭ при решении этой задачи использовалась триангуляция с адаптацией к особенности и без нее (см. табл. 11, столбцы 8–10) (см. [7]). В МКНК были взяты следующие значения весовых множителей ${{k}_{c}} = {{k}_{{{{m}_{0}}}}} = {{k}_{{{{m}_{2}}}}} = {{k}_{{{{b}_{0}}}}} = {{k}_{{{{b}_{1}}}}} = 1$, ${{k}_{{{{m}_{1}}}}} = 0.12,\;{{k}_{{{{m}_{3}}}}} = 0.125$ и $\epsilon {{ = 10}^{{ - 12}}}$.

Из представленных результатов следует, что в МКЭ прежде всего за счет адаптации удается построить более точное решение по сравнению с МКНК и МКЭ с использованием равномерной сетки (более подробно о порядках сходимости при решении рассматриваемой задачи МКЭ см. в [7]). Такая возможность повышения точности решения задач с особенностями в МКНК показана в [20]. Однако отметим, что МКНК здесь применялся к решению краевой задачи для уравнения высокого порядка в исходном виде (1.1)–(1.3), в то время как в МКЭ дополнительно приходится работать со слабыми постановками. При этом надо иметь в виду, что многими численными методами при решении уравнений низкого порядка точность и порядок сходимости будут выше при аналогичном способе аппроксимации уравнений высокого порядка, например, одинаковыми полиномами или конечно-разностными схемами на одинаковых шаблонах. Это было продемонстрировано при решении задачи изгиба шарнирно закрепленной прямоугольной пластины МКНК при $K = 4$ (см. [23]). В этом случае решение краевой задачи для бигармонического уравнения можно легко свести к последовательному решению двух уравнений Пуассона (см. [29]), которые имеют меньший порядок, что приводит к лучше обусловленной задаче по сравнению с исходной. При сравнении с решением бигармонического уравнения в исходном виде этот способ позволяет достичь существенно лучшей точности и более высокого порядка сходимости.

Кроме того, повышение гладкости решения задачи в МКЭ не увеличит порядка сходимости в случае применения линейных элементов (см. [7]) в отличие от МКНК, где он возрастет и, следовательно, погрешность будет уменьшаться быстрее при увеличении степеней базисных элементов полиномиального пространства для приближения искомого решения в ячейке и уменьшении шага сетки. Заметим также, что в МКНК достаточно просто увеличить степень аппроксимирующих полиномов по методике, предложенной в этой и предыдущих работах. В табл. 12 приведены результаты численных экспериментов с различными значениями $\alpha = 7{\text{/}}3,\;10{\text{/}}3,\;13{\text{/}}3,\;16{\text{/}}3$, при этом в точке $(0,0)$ в бесконечность обращаются производные третьего, четвертого, пятого и шестого порядков соответственно. При увеличении гладкости решения его точность и порядок сходимости $R$ возрастают, но $R$ в целом не превышает порядка гладкости решения.

4. ЗАКЛЮЧЕНИЕ

Разработаны новые h-, p- и hp-МКНК повышенной точности для численного решения неоднородных бигармонических уравнений в нерегулярных областях. Возможности предложенных здесь вариантов метода проверены на решении тестовых задач путем вычисления погрешности и порядка сходимости приближенного решения на последовательности сеток. Показано, что новые варианты МКНК во многих случаях не уступают по точности высокоточным известным результатам других авторов, полученным другими методами.

Успешное применение разработанных p-МКНКh и p-МКНК1 позволило обеспечить экспоненциальный порядок сходимости, поскольку в рассмотренных примерах пп. 3.1–3.6 решения были достаточно гладкими функциями. Однако решения реальных, практически важных задач (пример п. 3.7) зачастую не являются достаточно гладкими (см. также пример п. 3.8). Вследствие этого методы, основанные на полиномиальной аппроксимации высокого порядка, в том числе и предложенный здесь p-МКНК1, становятся менее эффективными и гибкими, например, по сравнению с h-МКНК и p-МКНКh.

Авторы вырaжают признательность С.В. Идимешеву, В.И. Исаеву и В.А. Колотилову за полезное обсуждение результатов работы.

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

  1. Mayo A. The fast solution of Poisson’s and the biharmonic equations on irregular regions // SIAM J. Numer. Anal. 1984. V. 21. № 2. P. 265–299.

  2. Lai M.C., Liu H.C. Fast direct solver for the biharmonic equation on a disk and its application to incompressible flows // App. Math. Comput. 2005. V. 164. № 3. P. 679–695.

  3. Chen G., Li Z., Lin P. A fast finite difference method for biharmonic equations on irregular domains and its application to an incompressible Stokes flow // Adv. Comput. Math. 2007. V. 29. № 2. P. 113–133.

  4. Ben-Artzi M., Chorev I., Croisille J.P., Fishelov D. A compact difference scheme for the biharmonic equation in planar irregular domains // SIAM J. Numer. Anal. 2009. V. 47. № 4. P. 3087–3108.

  5. Davini C., Pitacco I. An unconstrained mixed method for the biharmonic problem // SIAM J. Numer. Anal. 2000. V. 38. № 3. P. 820–836.

  6. Eymard R., Gallouët T., Herbin R., Linke A. Finite volume schemes for the biharmonic problem on general meshes // Math. Comp. 2012. V. 812. № 280. P. 2019–2048.

  7. Guo H., Zhang Z., Zou Q. A  ${{C}^{0}}$ Linear Finite Element Method for Biharmonic Problems // J. Sci. Comput. 2018. V. 74. P. 1397–1422.

  8. Shao W., Wu X., Chen S. Chebyshev tau meshless method based on the integration-differentiation for biharmonic-type equations on irregular domain // Engin. Anal. Bound. Elem. 2012. V. 36. № 12. P. 1787–1798.

  9. Jiang Y., Wang B., Yuesheng X. A fast Fourier–Galerkin method solving a boundary integral equation for the biharmonic equation // SIAM J. Numer. Anal. 2014. V. 52. № 5. P. 2530–2554.

  10. Mai-Duy N., See H., Tran-Cong T. A spectral collocation technique based on integrated Chebyshev polynomials for biharmonic problems in irregular domains // Appl. Math. Model. 2009. V. 33. № 1. P. 284–299.

  11. Shao W., Wu X. An effective Chebyshev tau meshless domain decomposition method based on the integration-differentiation for solving fourth order equations // Appl. Math. Model. 2015. V. 39. № 9. P. 2554–2569.

  12. Сорокин С.Б. Переобусловливание при численном решении задачи Дирихле для бигармонического уравнения // Сиб. журн. вычисл. матем. 2011. Т. 14. № 2. С. 205–213.

  13. Greenbaum A., Greengard L., Mayo A. On the numerical solution of the biharmonic equation in the plane // Physica D. 1992. V. 60. № 1–4. P. 216–225.

  14. Bialecki B. A fast solver for the orthogonal spline collocation solution of the biharmonic Dirichlet problem on rectangles // J. Comput. Phys. 2003. V. 191. № 2. P. 601–621.

  15. Gomez-Polanco A., Guevara-Jordan J.M., Molina B. A mimetic iterative scheme for solving biharmonic equations // Math. Comput. Model. 2013. V. 57. № 9–10. P. 2132–2139.

  16. Беляев В.А., Шапеев В.П. Варианты метода коллокации и наименьших невязок для решения задач математической физики в трапециевидных областях // Вычисл. технологии. 2017. Т. 22. № 4. С. 22–42.

  17. Шапеев В.П., Беляев В.A. Решение с повышенной точностью бигармонического уравнения в нерегулярных областях методом коллокации и наименьших квадратов // Вычисл. методы и программирование. 2018. Т. 19. № 4. С. 340–355.

  18. Shapeev V.P., Belyaev V.A., Golushko S.K., Idimeshev S.V. New possibilities and applications of the least squares collocation method // EPJ Web of Conferences. 2018. V. 173. № 01012. P. 01012-1–01012-8.

  19. Идимешев С.В. Модифицированный метод коллокаций и наименьших невязок и его приложение в механике многослойных композитных балок и пластин: Дис. канд. физ.-мат. наук. Новосибирск: Ин-т вычисл. технологий СО РАН, 2016.

  20. Исаев В.И., Шапеев В.П. Варианты метода коллокаций и наименьших квадратов повышенной точности для численного решения уравнений Навье–Стокса // Ж. вычисл. матем. и матем. физ. 2010. Т. 50. № 10. С. 1758–1770.

  21. Исаев В.И., Шапеев В.П., Еремин С.А. Исследование свойств метода коллокации и наименьших квадратов решения краевых задач для уравнения Пуассона и уравнений Навье–Стокса // Вычисл. технологии. 2007. Т. 12. № 3. С. 53–70.

  22. Голушко С.К., Идимешев С.В., Шапеев В.П. Метод коллокаций и наименьших невязок в приложении к задачам механики изотропных пластин // Вычисл. технологии. 2013. Т. 18. № 6. С. 31–43.

  23. Шапеев В.П., Брындин Л.С., Беляев В.A. Решение эллиптических уравнений в полигональных областях методом коллокации и наименьших квадратов // Вестник ЮУрГУ ММП. 2019. Т. 22. № 3. С. 140–152.

  24. Vorozhtsov E.V., Shapeev V.P. On the efficiency of combining different methods for acceleration of iterations at the solution of PDEs by the method of collocations and least residuals // Appl. Math. Comput. 2019. V. 363. № 124644. P. 1–19.

  25. Лаевский Ю.М. Метод конечных элементов (основы теории, задачи). Новосибирск: Новосиб. гос. ун-т, 1999. 166 с.

  26. Федоренко Р.П. Введение в вычислительную физику. М.: МФТИ, 1994. 527 с.

  27. Ramšak M., Škerget L. A subdomain boundary element method for high-Reynolds laminar flow using stream function–vorticity formulation // Int. J. Numer. Meth. Fluids. 2004. V. 46. № 8. P. 815–847.

  28. Saad Y. Numerical methods for large eigenvalue problems. Manchester Univ. Press, 1992. 346 p.

  29. Тимошенко С.П., Войновский-Кригер С. Пластины и оболочки. М.: Физматгиз, 1963. 636 с.

  30. Шапеев В.П., Шапеев А.В. Решение эллиптических задач с особенностями по схемам высокого порядка аппроксимации // Вычисл. технологии. 2006. Т. 11. Часть 2. Специальный выпуск. С. 84–91.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал