CC BY
14
Вычислительные технологии
Том 22, № 4, 2017
Варианты метода коллокации и наименьших невязок для решения задач математической физики в трапециевидных областях
В. А. Беляев1'2'*, В. П. Шапеев1'2
1 Новосибирский национальный исследовательский государственный университет, Россия 2Институт теоретической и прикладной механики СО РАН, Новосибирск, Россия *Контактный e-mail: belyaevasily@mail.ru
Предложены и реализованы новые варианты метода коллокации и наименьших невязок (коллокации и наименьших квадратов) повышенной точности для численного решения краевых задач для уравнений с частными производными в трапециевидных областях. Указаны достоинства предложенного подхода к решению с повышенной точностью задач для уравнений с частными производными в непрямоугольных областях. Реализация этого подхода и численные эксперименты выполнены на примерах решения уравнений Пуассона и бигармонического уравнения. Решение второго уравнения использовано для моделирования напряженно-деформированного состояния изотропной пластины, имеющей форму трапеции и находящейся под действием поперечной нагрузки.
Ключевые слова: метод коллокации и наименьших невязок, краевая задача, трапециевидная область, повышенный порядок аппроксимации, уравнение Пуассона, бигармоническое уравнение.
Введение
Существует немало задач численного моделирования, для которых при ограниченных вычислительных ресурсах можно добиться такой малой погрешности их численного решения на ЭВМ методом повышенной точности (МПТ), которая недостижима методом низкой точности (МНТ). И зачастую, когда требуется найти решение задачи с некоторой заданной погрешностью, применением МПТ за счет выбора параметров метода это можно сделать за существенно меньшее время работы ЭВМ, чем с использованием МНТ [1, 2]. Конечно, для многих задач техники не требуется повышенная точность расчета. Однако в истории астрономических открытий и в настоящее время в задачах моделирования физических процессов можно назвать ряд примеров, когда требовалось изобретать численные МПТ, чтобы решить конкретную задачу с требуемой точностью. Собственно, в самой вычислительной математике методы повышенной точности позволяют строить решения задач, которые могут применяться в качестве эталонов для выяснения возможностей существующих и вновь создаваемых методов [3-6].
Исследования по развитию МПТ направлены на расширение возможностей численных методов для решения практических задач. Применение МПТ может дополнить
© ИВТ СО РАН, 2017
и усилить быстро развивающееся для этих целей современное магистральное направление в области численного решения задач, основанное на распараллеливании алгоритмов. Многие численные методы, дающие высокую точность решения краевых задач для уравнений с частными производными (РВЕ) в прямоугольных областях, не имеют аналогов в случае областей иной формы. Однако большое количество практических задач требуется решать в областях, отличных от прямоугольных.
В данной работе краевые задачи для РВЕ решаются с повышенной точностью проекционным сеточным методом КНН (коллокации и наименьших невязок — вариантом метода коллокации и наименьших квадратов) [5, 7] в трапециевидных областях. В методе КНН задаче для РВЕ путем проектирования в конечномерное линейное функциональное пространство ставится в соответствие приближенная задача, решение которой сводится к решению системы линейных алгебраических уравнений (СЛАУ). Решение последней определяет приближенное решение дифференциальной задачи. Поскольку в методе КНН приближенное решение кусочно-аналитическое (в данной работе кусочно-полиномиальное), в нем относительно просто реализуются варианты метода в областях различной формы и на сетках с различной формой ячеек, в том числе варианты повышенной точности и на адаптивных сетках [5, 8, 9].
Здесь реализованы три разных варианта построения расчетной сетки в трапециевидных областях. Их возможности проверены на решении краевых задач для уравнения Пуассона и бигармонического уравнения, которое описывает напряженно-деформированное состояние (НДС) под различными нагрузками изотропной трапециевидной пластины. Задача о НДС прямоугольной пластины методом КНН с повышенной точностью решена в [10]. В настоящей работе на решении тестовой задачи для пластины в форме трапеции показана высокая точность предложенного варианта метода КНН. Показан высокий порядок сходимости численного решения по методу КНН при решении задачи о НДС трапециевидной пластины и в случае, когда неизвестно ее точное решение.
Описание, характеристика, история возникновения, некоторые возможности метода КНН и библиография по нему приведены в [5,6,9,11, 12-14]. Здесь коротко отметим, что в методе КНН для нахождения решения приближенной задачи в каждой ячейке выписывается переопределенная система уравнений, после подстановки в которую искомого приближенного решения из нее получается переопределенная СЛАУ для нахождения коэффициентов представления в линейном функциональном пространстве решения приближенной задачи. Решение переопределенной СЛАУ ищется из требования минимизации функционала невязки, составленного из невязок всех ее уравнений. В качестве функционала обычно берется сумма квадратов невязок. Его минимум может отыскиваться разными способами. В случае переопределенной СЛАУ он достигается на ее псевдорешении, построенном методом наименьших квадратов или ортогональным методом [9, 13]. Детальное описание метода КНН дано ниже.
1. Задача Дирихле для уравнения Пуассона
1.1. Постановка задачи и описание метода
Рассмотрим задачу Дирихле для уравнения Пуассона
Аь = f (х\, х2), (х\, х2) € П, Нгп = 9{хь х2)
в трапециевидной области О с границей 5О (рис. 1, а), где V(хх, х2) — искомая функция; /(х\, х2) и д(х\, х2) — заданные функции. Реализация рассматриваемого метода решения других краевых задач для уравнения Пуассона осуществляется аналогично алгоритму, описанному далее.
Пусть Нп — высота трапеции О, ап и Ьп — длины ее верхнего и нижнего оснований, при этом значения всех трех величин могут сильно различаться между собой. Чтобы покрыть О равномерной регулярной сеткой с трапециевидными ячейками, разобьем высоту трапеции на тх равных частей. Через концы отрезков на высоте трапеции проведем параллельные прямые до пересечения с ее боковыми сторонами. Верхнее и нижнее основания трапеции равномерно разобьем точками на т2 отрезков с длинами а = ап/т2 и Ь = Ьп/т2 соответственно. Числа тх и т2 выбираем такими, чтобы размеры ячеек по высоте Н = Нп/т1 и полусумме длин оснований сильно не различались. Пронумеруем полученные точки на нижнем и верхнем основаниях в возрастающей последовательности. Соединим отрезками прямых точки с одинаковыми номерами верхнего и нижнего оснований. Полученные трапециевидные ячейки, которые полностью расположены внутри области, назовем внутренними. Ячейки, стороны которых расположены около границы области 6О, назовем граничными. Для удобства реализации метода в каждой ячейке области введем локальные координаты
У1
Х\ — Хц
К '
У2
Х2 — Х2^
К '
(2)
Здесь (х^, хц) — центр у-й ячейки, у = , N = тх х т2 — количество ячеек;
2 К = д/Н(а + Ь)/2 — характерный линейный размер ячейки сетки, Н — высота ячейки; V(х\, х2) = и(у\, у2). Координаты центра ячейки в данном случае определяются как средние арифметические величины от соответствующих координат вершин трапеции. Задача (1) после замены (2) в локальных переменных примет вид
1 д2 и 1 д2 и
+
К2 дух
К2 ду22
/ (х1 (ух, У2 ),Х2 (у1, У2 )), (У1, У2 ) е Оj с О,
(3)
и = д(хх (ух ,У2 ),Х2 (У1 ,У2 )), (У1 ,У2) е 6 О п 6 О^, у = 1,...,И.
В каждой у-й ячейке сетки приближенное решение и^ задачи (3) ищем в виде линейной комбинации с неопределенными коэффициентами базисных элементов простран-
Рис. 1. Область решения задачи (а). Схема расположения точек коллокации •, точек записи условий согласования х и краевых условий □ в ячейке сетки (б)
ства полиномов четвертой степени от двух переменных. Здесь в качестве базисных элементов фг (г = 1,..., 15) взяты мономы:
15
иН] (У1, у2) = ^ & = ЪИ + у1 + ^ у2 + ^ у12 + Ь5] У^ +
г=1
О ООО о о о
+Ь&] У1У2 + Ь71У1 У2 + Ь^ У1У2 + Ц' У1 У2 + Ьюз У1 + ЬЦ1У2 + &12^ У1 У2 +
У1У23 + Ьщ У14 + bl5j У24.
Неизвестные коэффициенты Ь^ в методе КНН находятся из переопределенной системы уравнений, состоящей в каждой ячейке из уравнений коллокации, условий согласования (глобального решения задачи, искомого в виде кусков, связанных с центром соответствующей ячейки) на общих сторонах, принадлежащих двум соседним ячейкам, и краевых условий на 8П, если ячейка является граничной [5, 6, 9, 14].
Определим по следующему алгоритму расстановку точек коллокации, записи условий согласования в любой ячейке или краевых условий для граничной ячейки. Каждую ее сторону поделим на четыре равных ("частичных") отрезка. В середине каждого частичного отрезка запишем условия согласования или, соответственно, краевые условия, если взятая сторона принадлежит внешней границе.
Для расстановки точек коллокации в любой ячейке соединим отрезками прямых, параллельными оси у1, расположенные на ее левой стороне точки записи условий согласования (или краевых условий в случае граничной ячейки) с соответствующими точками на ее правой стороне. Затем каждый такой полученный отрезок поделим на четыре равных частичных отрезка. В середине каждого частичного отрезка запишем уравнения коллокации (рис. 1, б). Сформулированный алгоритм однозначно определяет в каждой ячейке координаты у1,у2 точек коллокации, запись условий согласования и краевых условий.
Можно отметить, что существует много различных способов расстановки точек для записи уравнений коллокации, условий согласования и краевых условий. При этом во многих случаях полученные приближенные решения будут с высокой точностью совпадать друг с другом. Расстановка точек для записи перечисленных уравнений, в некоторой мере близкая к симметричной и равномерной, в методе КНН дает хорошие результаты.
На рис. 2, а изображена расчетная область, покрытая регулярной сеткой с трапециевидными ячейками, в которых по описанному выше правилу расставлены точки записи уравнений коллокации, условий согласования и краевых условий. Для краткости изложения этот алгоритм построения сетки и расстановки точек, в которых записываются уравнения приближенной задачи, назовем вариантом А.
Уравнения коллокации в каждой ячейке выписываются в 16 точках коллокации и имеют вид
1 д 2и(У1С,У2с) 1 д2 и{у1С,у2С)
Ь2-д^2-+ V-- = f (Х1 (У1с ,У2с), (у 1с ,У2С)),
где (у1с,у2с), с =1,..., 16, — точки коллокации.
В качестве условия согласования решения в четырех указанных выше точках согласования на каждой общей стороне между соседними ячейками требуем непрерывность линейной комбинации с весами функции и и ее производной по нормали
Рис. 2. Расчетная область (размер сетки 3 х 3) (а). Здесь • — точки записи уравнений колло-кации, х — точки записи условий согласования, □ — точки записи краевых условий. Трапециевидная область с вытянутыми ячейками (б)
ди+ + ди к-\—--+ к2и+ = к\—--+ к2и
дпз
дщ
Здесь п^ — внешняя нормаль к границе ]-й ячейки; и+ и и~ — пределы функции при стремлении изнутри и извне к границе ]-й ячейки соответственно; к1 ,к2 — весовые параметры.
Если сторона ячейки совпадает с некоторым куском границы области, в четырех указанных выше точках на ней выпишем краевые условия
к3и = к3д(х1, х2),
где к3 — весовой параметр краевого условия в переопределенной системе уравнений.
Как показано в [15], значения весовых параметров влияют в некоторых пределах на обусловленность СЛАУ, из которой находятся коэффициенты приближенного решения, скорость сходимости итераций и величина погрешности приближенного решения задачи. Но точность приближенного решения определяется прежде всего порядком аппроксимации дифференциальной задачи приближенной задачей и, соответственно, зависит от значения старшей степени полиномов в базисе используемого для этого пространства [6]. Если решается хорошо обусловленная задача (например, задача Дирихле для уравнения Пуассона, краевая задача для уравнений Навье — Стокса при относительно небольшом числе Рейнольдса), то в этом случае можно положить к1 = 1, к2 = 1, к3 = 1 без необходимости поиска их оптимальных значений и получить методом КНН ее решение, сходящееся на последовательности сеток.
В каждой ячейке, объединяя уравнения коллокации, условия согласования и краевые условия (в случае граничных ячеек), относительно неизвестных коэффициентов Ь^ получим переопределенную СЛАУ вида
15
= Я
I = 1,..., 32,
3 = 1,...,^
(4)
г=1
которую для краткости назовем "локальной". Для того чтобы определить, что понимается под решением этой системы, рассмотрим функционал
32/15 \ 2
ф = ЕЕ — я
1=1 \г=1 )
ИЩ] —
Он представляет собой сумму квадратов невязок всех уравнений в ячейке на приближенном решении задачи. Решение системы (4) на каждой итерации находится из условия минимума функционала Ф по коэффициентам Ь^, г = 1,..., 15 [5]. Система уравнений, полученная объединением уравнений во всех ячейках расчетной области (глобальная СЛАУ), решается в итерационном процессе Гаусса — Зейделя. В ней одна "глобальная итерация" состоит из последовательного решения локальных СЛАУ (4) во всех ячейках области. Здесь и далее в условиях согласования значения решения в соседних ячейках брались с предыдущей итерации, если они не были сосчитаны на текущей итерации, или в противном случае — их значения с текущей итерации. При построении решения в каждой ячейке матрица системы из 32 уравнений приводится к верхнетреугольному виду ортогональным методом Гивенса [9, 13]. При этом в строках с номерами больше 15 в прямоугольной матрице локальной системы образуются нулевые строки, а компоненты столбца правых частей имеют значения невязок соответствующих им уравнений на решении подсистемы с полученной верхнетреугольной матрицей. В случае краевых задач для линейных РВЕ с достаточно гладкими решениями невязки уравнений приближенной задачи определяются прежде всего величиной погрешности аппроксимации решений дифференциальной задачи решениями приближенной задачи. В сходящемся процессе решения приближенной задачи на последовательности сеток при стремлении шага сетки к нулю невязки уравнений на ее решении стремятся к нулю, если имеет место аппроксимация с ненулевым порядком уравнений дифференциальной задачи уравнениями приближенной задачи. Поэтому при построении решения приближенной задачи уравнения локальной переопределенной системы с номерами больше 15 после приведения ее матрицы к верхнетреугольному виду не учитываются.
Итерационный процесс продолжается до тех пор, пока не выполнится условие
тах 1ЬцП+1 — Ь^а1 < с, ч
где 6™- — г-й (г = 1,..., 15) коэффициент полинома, аппроксимирующего решение в ячейке с номером ] на п-й итерации. Величина е — заданная константа, называемая псевдопогрешностью решения. Ее значение выбирается таким, чтобы погрешность решения глобальной СЛАУ была существенно меньше погрешности аппроксимации исходной задачи.
Замечания. Известно, что при аппроксимации в численных методах дифференциальных задач с применением сеток с сильно вытянутыми ячейками получается плохо обусловленная приближенная задача. Если верхнее основание трапеции О существенно меньше нижнего (рис. 2, б), то стремление избежать появления вытянутых ячеек может привести к неоправданно мелким ячейкам сетки (и соответственно большему их количеству) около верхнего основания. Например, это может быть, если там решение не имеет существенных особенностей. Также в этом варианте метода имеет место небольшая возможность добиваться некоторого "сгущения" сетки внутри области О за счет неравномерного расположения ее узлов на основаниях и сторонах трапеции, чтобы
при необходимости адаптировать сетку в некоторых подобластях О. В случае адаптации сетки к особенностям решения задачи возможно применение более эффективного способа адаптации, который применен в [8] для прямоугольных ячеек. Поэтому во многих случаях в целях получения решений задач в трапециевидных областях могут быть предпочтительнее изложенные ниже варианты построения сеток.
1.2. Численные эксперименты и результаты расчетов
В численных экспериментах на сходимость метода на последовательности сеток при измельчении шагов сетки вдвое рассматривалась задача (1) в трапециевидных областях. В качестве тестовых решений брались аналитические функции, подстановкой которых в уравнения подбиралась их правая часть (см. разд. 2) для того, чтобы эти функции удовлетворяли рассматриваемым дифференциальным уравнениям. Краевые условия брались из значений тестовых функций и их производных на 80. Во всех представленных ниже таблицах в этом разделе приведены значения погрешности приближенного решения в равномерной норме
II U - uh ||с= max max Iu(xls,x2s) - uhi(xis,X2S)l
i=1,...N s=l,...Li
где N — количество расчетных ячеек в области; Li — количество равномерно распределенных точек (xls,x2s), взятых в г-й ячейке для подсчета в них погрешности; и — точное решение задачи (1); Uh — приближенное решение в области 0; Uhi — приближенное решение в г-й ячейке. Порядок сходимости определяется следующим образом:
ENl/2,N2/2 . , lo§2 7?-, (5)
где — значение абсолютной погрешности на сетке размера N = Nl х N2,
En1/2,n2/2 — значение абсолютной погрешности на сетке ^l/2 х N2/2. В численных экспериментах в одном случае взяты следующие значения параметров: Li = 100, е = 0.5 • 10-14, kl = 1, к2 = 5h, к3 = 1.5. Порядок сходимости такой же, как и при kl = 1, к2 = 1, к3 = 1, но погрешность во втором случае незначительно больше, чем в первом. Во всех численных экспериментах здесь в начальном приближении решения взяты bij = 0.4, где г =1,..., 15, j = 1,..., N. Такое задание начальных значений bij дало некоторое уменьшение числа итераций по сравнению со случаем начальных значений bij = 0. В расчетах использовано комбинированное применение операции продолжения на многосеточном комплексе в методе Федоренко и метода ускорения сходимости итерационного процесса, основанного на подпространствах Крылова [16-18].
Из табл. 1-4 видно, что численное решение имеет порядок сходимости не хуже четвертого, за исключением случая перехода с сетки 40 х 40 на сетку 80 х 80 в численном эксперименте с результатами, приведенными в табл. 3. Здесь при стремлении шага сетки h ^ 0 появляются сильно вытянутые ячейки и ухудшается обусловленность приближенной задачи.
1.3. Алгоритм с перестройкой ячеек сетки
Рассмотрим другие два варианта алгоритма реализации метода КНН решения задач для уравнений с частными производными в трапециевидных областях в пространстве
Таблица 1. Результаты численного эксперимента с тестовым решением и(х1, Ж2) = 8т(ж1) + вт(ж2) в равнобедренной трапециевидной области с вершинами (0, 0), (1, 0), (0.9, 1), (0.1, 1)
Размер сетки Погрешность Количество итераций Порядок сходимости
10 х 10 1.83 ■ 10-9 88 —
20 х 20 9.55 ■ 10-11 111 4.26
40 х 40 5.28 ■ 10-12 272 4.17
80 х 80 2.41 ■ 10-13 625 4.45
Таблица 2. Результаты численного эксперимента с тестовым решением и(х1, Ж2) = ео8(ж1) + еов(ж2) в "скошенной" трапециевидной области с вершинами (0.4, 0), (0.7, 0), (1.9, 2), (1.8, 2)
Размер сетки Погрешность Количество итераций Порядок сходимости
10 х 5 1.68 ■ 10-8 87 —
20 х 10 5.78 ■ 10-10 117 4.86
40 х 20 1.88 ■ 10-11 171 4.94
80 х 40 7.76 ■ 10-13 441 4.59
Таблица 3. Результаты численного эксперимента с тестовым решением и(ж1, Ж2) = ех1 + ех2 в трапециевидной области с вершинами (0, 0), (1, 0), (0.51, 1), (0.49, 1), у которой длины верхнего и нижнего оснований сильно различаются
Размер сетки Погрешность Количество итераций Порядок сходимости
10 > 10 4.86 ■ 10-9 281 —
20 х 20 1.92 ■ 10-10 521 4.66
40 х 40 7.44 ■ 10-12 781 4.68
80 х 80 7.61 ■ 10-13 2231 3.28
Таблица 4. Результаты численного эксперимента с тестовым решением и(х1, Ж2) = е10ж1 + е10ж2 + е5ж1+5ж2, имеющим большие градиенты, в трапециевидной области с вершинами (0, 0), (1, 0), (1, 1), (0.4, 1). Здесь е = 0.5 ■ 10-10
Размер сетки Погрешность Количество итераций Порядок сходимости
10 х 10 1.61 ■ 10-1 141 —
20 х 20 1.00 ■ 10-2 229 4.00
40 х 40 6.23 ■ 10-4 845 4.00
80 х 80 3.38 ■ 10-5 4061 4.20
полиномов четвертой степени. Аналогично можно реализовать эти варианты метода КНН в пространстве полиномов любой степени.
При появлении вытянутых, а также малых ячеек рядом с достаточно крупными по сравнению с ними ухудшается обусловленность глобальной СЛАУ. В итоге увеличивается (иногда довольно сильно) количество итераций, необходимых для построения
приближенного решения. Например, это наблюдается в трапециевидной области (в случае построения сетки в варианте А) с сильно отличающимися длинами оснований по отношению друг к другу (рис. 2, б). В подобных случаях (и не только) необходим более эффективный вариант метода.
Впишем расчетную область О в прямоугольник, две стороны которого параллельны оси х1 и проходят через верхнее и нижнее основания трапеции, а две другие стороны, параллельные оси х2, проходят через принадлежащие 8О точки с минимальной и максимальной координатами по оси х1. Полученную прямоугольную область покроем регулярной сеткой с прямоугольными ячейками размера 2Н1 х 2к2. Ячейки, пересеченные границей области, для краткости здесь назовем граничными. Часть граничной ячейки, отсеченную границей и лежащую внутри области, назовем нерегулярной ячейкой. Часть граничной ячейки, лежащую вне области, назовем "законтурной". Прямоугольную ячейку сетки, от которой границей области отсечена нерегулярная ячейка, назовем коротко материнской. Назовем внешней стороной нерегулярной ячейки часть границы области, которая оказалась внутри граничной ячейки. Стороны и отсеченные границей части сторон граничной ячейки, расположенные внутри области, назовем внутренними сторонами нерегулярной ячейки. Решение в прямоугольных ячейках, которые не содержат часть области О, естественно, не строим.
Во внутренних и граничных ячейках введем локальную систему координат с началом в центре ячейки
Х1 — Х^ Х2 — Х2]
У1 = —т-, У2 =-т-,
П1 П2
где (х^, х2^) — центр ]-й ячейки, ] = 1,..., N, N — количество ячеек. Решаемое дифференциальное уравнение запишем в локальных координатах.
В следующем варианте алгоритма реализации метода КНН в трапециевидной области при формулировке СЛАУ приближенной задачи в рассеченных границей области ячейках регулярной сетки используются расположенные в области их нерегулярные части вместе с их законтурными частями. Предполагается, что решение задачи не имеет особенностей на границе и в некоторой ее малой окрестности и дифференциальное уравнение задачи справедливо в используемой окрестности границы. В данном варианте аналогично тому, как это делалось в [8, 9], преследуется цель избежать необходимости построения решения в малых нерегулярных ячейках, что зачастую ухудшает обусловленность СЛАУ приближенной задачи и точность ее решения.
Рис. 3. Типы ячеек: а — внутренняя ячейка; б, в и г — типы граничных ячеек. Здесь • — точки записи уравнений коллокации, х — точки записи условий согласования, □ — точки записи краевых условий
Каждую сторону ячейки, общую с соседними прямоугольными ячейками, делим на четыре равных отрезка. В серединах этих отрезков расставим точки записи условий согласования (рис. 3). Отрезок границы области 80 внутри каждой рассеченной граничной прямоугольной ячейки делим либо на четыре равных частичных отрезка, если соответствующая нерегулярная ячейка (часть рассеченной ячейки, расположенная внутри области) — трапеция (рис. 3, б), либо на восемь равных частичных отрезков, если соответствующая нерегулярная ячейка — треугольник (рис. 3, г). В серединах частичных отрезков расставим точки записи краевых условий. Пусть нерегулярная угловая в области ячейка является трапецией и длина ее основания на границе области больше, чем шах(^/4, К2/4), тогда на этом основании расставим четыре точки записи граничных условий (рис. 3, в). Если длина рассматриваемого основания меньше, чем шаах(к1/4,к2/4), то на боковой стороне угловой ячейки на границе области, как в случае нерегулярной треугольной ячейки, расставим восемь точек записи краевых условий (рис. 3, г).
Левую и правую стороны всех прямоугольных ячеек (в случае граничной ячейки — левую и правую стороны материнской ячейки) поделим на четыре равных отрезка. Каждый отрезок прямой, соединяющий соответствующие середины этих отрезков, поделим на четыре равных частичных отрезка. Середину каждого такого частичного отрезка назовем точкой коллокации (всего 16 точек коллокации). В каждой точке колло-кации запишем уравнение коллокации (рис. 3).
Для удобства этот алгоритм построения сетки и решения задачи назовем вариантом Б. Будем считать несамостоятельными нерегулярные ячейки, расположенные в таких граничных ячейках, в которых начало локальной системы координат находится вне расчетной области. Треугольную ячейку всегда считаем несамостоятельной. В следующем варианте алгоритма (вариант В) несамостоятельная ячейка треугольного вида всегда присоединяется к соседней с ней, расположенной внутри области, или к самостоятельной граничной ячейке, именно к той, с которой имеет общую сторону наибольшей длины.
Вместо четырех условий согласования, записанных в точках на общей стороне объединяемых ячеек, в локальную систему объединенной ячейки добавляются столько краевых условий, записанных на внешней стороне присоединяемой несамостоятельной ячейки, чтобы число уравнений в системе всегда было 32. Далее решение строится во всех ячейках внутри области и в самостоятельных околограничных ячейках в их локальной системе координат. В несамостоятельные ячейки продолжаются решения, записанные в локальных системах координат ячеек, к которым присоединены несамостоятельные ячейки. То есть в объединенной ячейке — единые локальная система координат и формула решения. При этом в самостоятельной ячейке, соседней с несамостоятельной, в точках, расположенных на их общей стороне, решение согласуется с решением в объединенной ячейке.
В рамках численного эксперимента с использованием варианта В метода КНН, аналогичного эксперименту, проведенному с применением варианта А, результаты которого приведены в табл. 3, на сетке 80 х 80 достигнута точность решения 2.23 ■ 10-12 за 787 итераций. Решение в случае применения варианта А на сетке 80 х 80 получено с точностью 7.61 ■ 10-13, но для ее достижения потребовалась 2231 итерация.
Также были реализованы различные варианты метода КНН с использованием пространства полиномов второй степени. Одним из них был вариант А. В другом любая несамостоятельная ячейка присоединялась к самостоятельной по правилу, описанному
выше. Но при этом ее законтурная часть отбрасывалась и никак не использовалась. В результате численного эксперимента с применением этого варианта метода с тестовым решением и областью, указанными в заголовке табл. 3, достигнута несколько лучшая точность решения за существенно меньшее количество итераций, чем в случае применения варианта А. Из-за ограничения на объем статьи таблицы с этими результатами здесь не приводим.
Следует отметить, что приближенная задача Дирихле для уравнения Пуассона в методе КНН достаточно хорошо обусловлена [15]. В результате многих численных экспериментов с вариантом А установлено, что в случае хорошо обусловленной приближенной задачи наблюдается сходимость численного решения даже при появлении таких вытянутых ячеек, наличие которых в случае плохо обусловленной задачи приводит к расходимости численного решения. Однако в численных экспериментах в последнем случае нередко по варианту В метода КНН эффективно строилось решение приближенной задачи. Здесь в качестве плохо обусловленных рассматривались приближенные задачи для бигармонического уравнения. Некоторые характерные сравнительные результаты численных экспериментов приведены далее.
2. Бигармоническое уравнение 2.1. Постановка задачи и описание метода
Пусть требуется определить напряженно-деформированное состояние упругой изотропной трапециевидной пластины, закрепленной на краях и находящейся под действием поперечных нагрузок произвольного вида. При этом нагрузки такие, что материал деформированной пластины остается упругим. Пластина предполагается тонкой, что позволяет свести трехмерную задачу теории упругости к двумерной задаче теории тонких пластин. Следуя Тимошенко [19], разрешающую систему уравнений запишем в перемещениях, определив которые, можно полностью восстановить НДС пластины. В этом случае величины прогибов пластины определяются из решения бигармонического уравнения
д4-ю(х1 ,Х2) д41^(х1,х2) д4 т(х1,х2) д(х1,Х2) . . ^ ..
+2 ~Щдх2~ + аг4 , (Х1,Х2) е 0 (6)
Здесь т(х1,х2) — прогиб срединной поверхности пластины; д(ж1 ,х2) — поперечная нагрузка; И = Ек3/12(1 — V2) — жесткость пластины при изгибе, К — толщина пластины, Е,и — модуль Юнга и коэффициент Пуассона изотропного материала пластины; трапециевидная область 0 — проекция срединной поверхности пластины на плоскость (ж1,ж2). На каждом участке границы (края) пластины может быть задано одно из следующих условий закрепления:
дт
т = 0,^— = 0 — защемленный край; (7)
оп
т = 0, Мпчи = 0 — шарнирно закрепленный край; (8)
Мпчи = 0, = 0 — свободный край. (9)
Здесь и далее п — внешняя нормаль к границе области д0; Мп — дифференциальный оператор второго порядка; Уп — дифференциальный оператор третьего порядка. Здесь дифференциальные операторы Мп, Уп взяты из [19]. Условие (9) может быть задано
только на некоторой части границы, а на оставшейся ее части должны быть заданы другие условия, иначе задача некорректна.
Область О покрывается расчетной сеткой с прямоугольными ячейками О^ (г = 1,..,Ы, N — число расчетных ячеек сетки). Далее применим к поставленной задаче алгоритм, аналогичный описанному выше варианту В с использованием пространства полиномов четвертой степени. Отличия будут в записи уравнений дискретной задачи.
В каждой ячейке О^ делается переход к локальным переменным (у1,у2). Для определения решения и^ в О^ на к-й итерации выписываются: уравнения коллокации
Ь2 д4тк К1 д4тк фъЫ
-А + 2 + ти^т = , (У1,У2) е О
К2 дУ4 ду2ду22 К2 ду4 1 2 Б
— требования удовлетворения приближенным решением уравнению (6) в заданных точках коллокации;
условия согласования
к 1 дтк 1
™к + 'Р1Т--= ШаЛ] + Р1Т--о-, (У1,У2) е дОг\дО,
Кп оп Кп оп
д21ик 1 д ътк д21иаЛ] 1 д3таЛ]
-¿Т + Р2Т о з~ = О 2 + Р2Т—^^ (У1,У2) е дОi\дО, оп2 Кп оп3 ОП2 Кп оп3
— условия склейки решения тг с решением в соседней ячейке в точках согласования на общей с соседней ячейкой стороне, где Кп = К1, если направление п совпадает с направлением оси у1, и Кп = К2, если направление п совпадает с направлением оси у2;
краевые условия в указанных выше точках записи краевых условий.
В данном случае в каждой ячейке выписывается переопределенная система из 48 уравнений: 16 уравнений коллокации и по два уравнения в каждой точке для записи условий согласования или краевых условий.
Отметим, что в методе КНН точность решения приближенной задачи зависит прежде всего от размера базиса используемого пространства полиномов, в которую проектируем исходную дифференциальную задачу. Далее точность решения зависит от способов расстановки точек для записи уравнений коллокации, условий согласования и краевых условий, а также от количества уравнений в переопределенной СЛАУ приближенной задачи. В связи с этим возникает естественный вопрос, насколько должна быть переопределена СЛАУ [15]?
В качестве примера рассмотрим шарнирно закрепленную на краях прямоугольную пластину размера 10 х 10, находящуюся под действием нагрузки д (Па) вида
105 яп
( ЖхЛ I ПХ2\ т - вт - .
I10/ ^
Покроем область прямоугольной сеткой, запишем уравнения коллокации, условия согласования, краевые условия (в случае граничных ячеек) способом, описанным выше, и в каждой ячейке вычислим число обусловленности матрицы локальной СЛАУ (как отношение ее максимального сингулярного числа к минимальному сингулярному числу). Затем выберем среди всех полученных сингулярных чисел максимальное. В численных
Ч
Таблица 5. Зависимость числа обусловленности от количества уравнений
Количество Число
уравнений обусловленности
48 3142.24
40 3231.52
38 8036.48
36 9366.90
34 12 910.44
32 те
экспериментах берем сетку размера 80 х 80. При этом каждый раз число уравнений кол-локации остается неизменным (16 уравнений), число условий согласования и краевых условий (в случае граничных ячеек) меняем (табл. 5). Уменьшение числа уравнений в локальных СЛАУ, начиная с 38, сначала привело к значительному росту числа обусловленности, увеличению погрешности решения задачи и, наконец, при 32 уравнениях к расходимости итерационного процесса.
Были проведены аналогичные вычисления числа обусловленности и для приближенной задачи Дирихле в случае уравнения Пуассона. Число обусловленности оказалось на два десятичных порядка меньше, чем для приближенной задачи для бигармониче-ского уравнения. Имеет место такое же поведение обусловленности глобальной СЛАУ, а также обусловленности СЛАУ приближенной задачи при решении уравнений На-вье — Стокса и других уравнений с частными производными [15]. Последнее установлено в результате аналогичных численных экспериментов вне рамок данной работы. Этот результат показывает, какую роль в методе КНН играют переопределение локальных СЛАУ и минимизация функционала невязки приближенной задачи, а также почему метод КНН более эффективен для решения дифференциальных уравнений, чем метод коллокации.
2.2. Численные эксперименты
2.2.1. Пластина под специальной нагрузкой
Рассмотрим шарнирно закрепленную на краях прямоугольную пластину размера ^ х , находящуюся под действием нагрузки д (Па) вида
4 = 1°5 -¿ту 81Ч М
В этом случае задача имеет точное аналитическое решение [17]
дс1\ ¿А (^ + ¿2)
гш(х1,х2) = ^ / 1-(10)
Тест 1. Пусть проекция срединной поверхности пластины на плоскость (ж1, х2) имеет форму трапециевидной области О с вершинами (0, 0), (10, 0), (9, 10), (1, 10). Для того чтобы воспользоваться приведенным аналитическим решением в качестве теста, в численном эксперименте на границе 8О зададим точные значения функции т (10)
и изгибающего момента Ранее в [10] методом КНН была решена задача для пря-
моугольной пластины. В расчетах использовались следующие параметры: с11 = 10 м, <12 = 10 м, И = 0.1 м, Е = 200 ГПа, V = 0.28, е = 10-10. Значения параметров р1 и р2 в условиях согласования взяты из работы [10]. Для определения порядка сходимости проведены численные эксперименты на последовательности сеток с уменьшением шага сетки вдвое. Для расчета погрешности в каждой расчетной ячейке равномерно распределялись 100 точек. Относительная погрешность численного решения во всей области вычислялась в этих точках по формуле
Е = шахХ1,Х2 1т(х1,х2) - и)н(х1 ,х2)1 М1,М2 тахх1 ,х2 Щхьх2)| ,
где тк(х1 ,х2) — приближенное решение, полученное методом КНН; N = N х N — размер сетки. Порядок сходимости решения в данной работе формально вычислялся по формуле (5), хотя в случае нерегулярной сетки, строго говоря, определение порядка сходимости отсутствует. В расчетах использовано комбинированное применение операции продолжения на многосеточном комплексе в методе Федоренко и метода ускорения сходимости итерационного процесса, основанного на подпространствах Крылова [16].
Тест 2. Проведен численный эксперимент в случае пластины, проекция срединной поверхности которой на плоскость (х1,х2) имеет форму "скошенной" трапециевидной области О с вершинами (0, 0), (10, 0), (3, 10), (1, 10). В этом тесте взяты данные и значения параметров, которые использовались в тесте 1. В табл. 6 приведены результаты численных экспериментов, полученные применением вариантов А, Б и В метода КНН.
Т а б лица 6 . Результаты численных экспериментов
Вариант Размер сетки Погрешность Количество итераций Порядок сходимости
Тест 1
20 20 1.38 10-3 321 —
А 40 > 40 2.38 10-4 537 2.53
80 80 5.15 10-5 1321 2.20
20 > 20 2.05 10-3 321 —
Б 40 40 1.59 10-4 561 3.68
80 > 80 2.82 10-5 845 2.82
20 20 2.25 10-3 610 —
В 40 > 40 1.42 10-4 813 3.98
80 80 1.02 10-5 1452 3.79
Тест 2
20 > 20 2.4 ■ 10-3 394 —
А 40 40 5.47 10-4 832 2.13
80 > 80 1.16 10-4 2591 2.23
20 20 6.15 10-4 241 —
Б 40 > 40 Расходимость — —
80 80 Расходимость — —
20 > 20 5.18 10-3 415 —
В 40 40 3.37 10-5 862 3.94
80 > 80 4.59 10-6 1848 2.87
Вычислим количество ячеек, в которых числа обусловленности локальных переопределенных СЛАУ имеют порядки 10,102,103 и т.д. Для каждого такого порядка количество соответствующих ячеек обозначим через Мс. В табл. 7 приведены значения Мс на последовательности сеток для вариантов А, Б и В в случае теста 2. Из полученных результатов видно, что сходимость решения начинается уже на достаточно грубых сетках. При этом, когда трапециевидная область близка к прямоугольной (тест 1), результаты численных экспериментов во всех вариантах А, Б и В близки друг к другу. В "скошенной" трапециевидной области (тест 2) при измельчении шагов сетки в случае варианта А появляются сильно вытянутые ячейки, из-за чего ухудшается обусловленность глобальной СЛАУ. В итоге уменьшается точность решения и увеличивается количество итераций, необходимых для его построения.
В варианте Б при стремлении шагов сетки к1 ^ 0, к2 ^ 0 появляются маленькие нерегулярные треугольные ячейки и соответствующие им большие законтурные части ячеек. Ввиду того что у таких ячеек слишком маленькие внешние стороны, при записи на них краевых условий получаются уравнения, близкие к линейно зависимым. Поэтому глобальная СЛАУ приближенной задачи становится плохо обусловленной. В этих случаях для решения задач целесообразно использовать вариант В метода КНН. Хотя, как это следует из приведенных выше примеров, нередко могут иметь место задачи, для решения которых варианты А и Б могут быть почти столь же эффективны, как и вариант В. Естественно, что могут быть предложены какие-то другие варианты метода КНН, в которых несамостоятельные ячейки любого вида присоединяются к соседним самостоятельным ячейкам с целью уменьшения обусловленности приближенной задачи [9]. Скорее всего, алгоритмы построения расчетной сетки в таких вариантах окажутся сложнее по сравнению с рассмотренными в данной работе.
На рис. 4, а показана форма прогиба пластины. Для иллюстрации напряженного состояния пластины (рис. 4, б) используем функцию распределения интенсивности напряжений
Таблица 7. Количество ячеек в трех вариантах метода КНН, локальные матрицы которых имеют разные числа обусловленности
Размер с обусловленностью порядка
сетки 10 102 103 104 105 1010
Вариант А
10X10 5 95 0 0 0 0
20X20 0 398 2 0 0 0
40X40 0 1475 125 0 0 0
80X80 0 6084 185 131 0 0
Вариант Б
10X10 1 67 4 1 0 0
20X20 0 255 9 2 0 0
40X40 0 976 21 15 0 1
80X80 0 3728 175 35 7 1
Вариант В
10X10 3 62 0 0 0 0
20X20 2 248 0 0 0 0
40X40 0 955 25 0 0 0
80x80 0 3660 205 15 0 0
Рис. 4. Величина прогиба и> деформированной пластины (а) и интенсивность напряжений I в случае использования варианта В (б)
I = (ах - оу)2 + о2х + ау + 6аху2
_ _ _ _ Е
ах = -Е (wxx + uwyy), Оу = -Е (wyy + = -Е (wxy), Е = --у.
1 — v
Отметим, что одним важным преимуществом предложенного варианта метода КНН является то, что без особого труда можно вычислить любые производные полученного численного решения, так как это решение является кусочно-полиномиальным. Расчеты проводились на одном вычислительном узле персонального компьютера, имеющего частоту 2.80 ГГц.
Во всех расчетах, приведенных ниже, использован вариант В метода КНН.
2.2.2. Пластина, защемленная по трем краям
Рассмотрим задачу, аналитическое решение которой неизвестно. Пусть трапециевидная изотропная пластина находится под действием равномерной поперечной нагрузки q = const. Вершины трапеции в точках (0, 0), (10, 0), (8, 8), (3, 8). Три края пластины защемлены (7), а оставшийся край свободен (9). В расчетах использована сетка 20x20. На рис. 5 показаны форма прогиба w и интенсивность напряжений I при параметрах задачи h = 0.1 м, Е = 200 ГПа, v = 0.28, q = 0.1 МПа. Для оценки порядка сходимости
Таблица 8. Результат численного эксперимента
У WNl,N2 - WNl/2,N2/2 Ус Величина нормы Порядок сходимости
|| Wi0,i0 - W55 ||с 1.80 ■ 10-2 —
У W20,20 — ^10,10 Ус 4.05 ■ 10-3 2.15
У W40,40 — W20,20 |с 8.66 ■ 10-4 2.22
|| Ws0,80 — W40,40 Ус 1.98 ■ 10-4 2.12
численного решения в равномерной норме на последовательности сеток вычислялась величина
|| WNl,N2 -WNi/2,N2/2 ||с= max lwNl,N2 (XXi,X2i) - WNl/2,N2 /2(Xii,X2i)l,
i=l,...L
где L — количество равномерно распределенных точек (xii,x2i), взятых для простоты подсчета в прямоугольнике, включающем рассматриваемую трапециевидную область (в точках, лежащих вне исходной трапециевидной области, погрешность полагалась равной нулю); wNl,N2 — численное решение, полученное на сетке Ni х N2; wNl/2,N2/2 — численное решение на сетке Ni/2 х N2/2. В данном примере L = 6400 (табл. 8). Порядок сходимости решения формально определялся по формуле
, 11 WNI,N2 — WNI/2,N2/2 11 с
bg2,-Г,
|| wN1/2,N2/2 — wNl/4,N2/4 ||с
где wNl/4,N2/4 — численное решение на сетке Ni/4 х N2/4.
2.2.3. Пластина, шарнирно закрепленная по двум краям
Рассмотрим задачу, аналитическое решение которой неизвестно. Пусть трапециевидная изотропная пластина находится под действием равномерной поперечной нагрузки q = const. Вершины трапеции: (0, 0), (10, 0), (8, 9), (1, 9). Два края пластины шарнирно закреплены (8), а два других края свободны (9). В расчетах использована сетка 20 х 20. На рис. 6 показаны форма прогиба w и распределение на ней интенсивности напряжений при параметрах задачи h = 0.1 м, Е = 200 ГПа, v = 0.28, q = 50 КПа.
а б
L ГПа
Заключение
В работе реализовано несколько различных подходов к построению вариантов метода КНН решения краевых задач для PDE в трапециевидных областях. Показано, что их приближенные решения, полученные c применением этих вариантов, сходятся с повышенным порядком, а в случаях тестовых численных экспериментов они с высокой точностью совпадают с известными решениями. Важным обстоятельством является наличие кусочно-аналитического решения задачи в методе КНН, которое не столь жестко, как в некоторых других методах, связано с формой ячейки расчетной сетки и конфигурацией ее узлов. Это позволяет внутри области брать регулярную сетку, а на границе относительно просто использовать "одинарный слой" нерегулярных ячеек. Кроме того, регулярность сетки внутри области существенно упрощает программную реализацию решения приближенной задачи, например, в сравнении со случаем, когда нужно запрограммировать и реализовать для каждой ячейки поиск соседних ячеек на нерегулярной сетке. Последнее приводит к увеличению кода программы и времени решения задачи на ЭВМ.
Установлено, что обусловленность приближенной задачи зависит прежде всего от обусловленности дифференциальной задачи, которую она аппроксимирует, от геометрии области и существенно от способа построения расчетной сетки. В конечном итоге от указанных факторов зависит возможность эффективного численного решения методом КНН краевых задач для уравнений с частными производными. Как показали численные эксперименты, с применением алгоритма В за счет нового подхода построения сетки и присоединения несамостоятельных ячеек к самостоятельным почти всегда удавалось в большей степени уменьшить обусловленность СЛАУ приближенной задачи, чем в случаях применения алгоритмов А и Б. Поэтому для хорошо обусловленных задач предпочтительным является вариант А ввиду его простоты, а для плохо обусловленных задач рекомендуется применять вариант В.
Показанная здесь возможность решать задачи теории упругости с повышенной точностью в неканонических областях имеет самостоятельное прикладное значение. С другой стороны, проделанная работа может рассматриваться как методическая с точки зрения поиска эффективных и достаточно просто реализуемых алгоритмов решения с повышенной точностью PDE в неканонических областях с более сложными границами, состоящими из сочетания криволинейных и полигональных участков. Скорее всего, в этом случае будет эффективно сочетание разных подходов, которые проверялись в [9], в данной работе и, возможно, появятся в ходе дальнейшего поиска.
Благодарности. Работа выполнена при поддержке РФФИ (грант № 13-01-00227-a). Список литературы / References
[1] Липавский М.В., Толстых А.И. Об одной мультиоператорной схеме десятого порядка и ее применении в прямом численном моделировании // Журн. вычисл. математики и матем. физики. 2013. Т. 53, № 4. C. 600-614.
Lipavskii, M.V., Tolstykh, A.I. Tenth-order accurate multioperator scheme and its application in direct numerical simulation // Comput. Math. and Math. Phys. 2013. Vol. 53, No. 4. P. 455-468.
[2] Ш^апеев А.В., ШШапеев В.П. Разностные схемы повышенной точности для решения эллиптических уравнений в области с криволинейной границей // Журн. вычисл. математики и матем. физики. 2000. Т. 40, № 2. С. 223-232.
Shapeev, A.V., Shapeev, V.P. High-order accurate difference schemes for elliptic equations in a domain with a curvilinear boundary // Comput. Math. and Math. Phys. 2000. Vol. 40, No. 2. P. 213-221.
[3] Botella, O., Peyret, R. Benchmark spectral results on the lid-driven cavity flow // Computers & Fluids. 1998. Vol. 27, No. 4. P. 421-433.
[4] Shapeev, A.V., Lin, P. An asymptotic fitting finite element method with exponential mesh refinement for accurate computation of corner eddies in viscous flows // SIAM J. Sci. Comput. 2009. No. 31(3). P. 1874-1900.
[5] Слепцов А.Г. Коллокационно-сеточное решение эллиптических краевых задач // Моделирование в механике: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ВЦ; ИТПМ. 1991. Т. 5(22), № 2. С. 101-126.
Sleptsov, A.G. Collocation-grid solution of elliptic boundary-value problems // Modelirovanie v Mekhanike. 1991. Vol. 5(22), No. 2. P. 101-126. (In Russ.)
[6] Исаев В.И., Ш^апеев В.П. Варианты метода коллокаций и наименьших квадратов повышенной точности для численного решения уравнений Навье — Стокса // Журн. вычисл. математики и матем. физики. 2010. Т. 50, № 10. С. 1758-1770.
Isaev, V.I., Shapeev, V.P. High-accuracy versions of the collocations and least squares method for the numerical solution of the Navier —Stokes equations // Comput. Math. and Math. Phys. 2010. Vol. 50, No. 10. P. 1670-1681.
[7] Слепцов А.Г., Шокин Ю.И. Адаптивный проекционно-сеточный метод для эллиптических задач // Журн. вычисл. математики и матем. физики. 1997. Т. 37, № 5. С. 572-586. Sleptsov, A.G., Shokin, Yu.I. An adaptiv grid-projection method for elliptic problems // Comput. Math. and Math. Phys. 1997. Vol. 37, No. 5. P. 558-571.
[8] Беляев В.В., ШШапеев В.П. Метод коллокаций и наименьших квадратов на адаптивных сетках в области с криволинейной границей // Вычисл. технологии. 2000. Т. 5, № 4. С. 13-21.
Belyaev, V.V., Shapeev, V.P. The collocation and least squares method on adaptive grids in a domain with a curvilinear boundary // Comput. Technologies. 2000. Vol. 5, No. 4. P. 13-21. (In Russ.)
[9] ШШапеев В.П., Беляев В^. Варианты метода коллокации и наименьших невязок повышенной точности в области с криволинейной границей // Вычисл. технологии. 2016. Т. 21, № 5. С. 95-110.
Shapeev, V.P., Belyaev, V.A. Versions of high order accuracy collocation and least residuals method in the domain with a curvilinear boundary // Comput. Technologies. 2016. Vol. 21, No. 5. P. 95-110. (In Russ.)
[10] Голушко С.К., Идимешев С.В., ШШапеев В.П. Метод коллокаций и наименьших невязок в приложении к задачам механики изотропных пластин // Вычисл. технологии. 2013. Т. 18, № 6. С. 31-43.
Golushko, S.K., Idimeshev, S.V., Shapeev, V.P. Application of collocations and least residuals method to problems of the isotropic plates theory // Comput. Technologies. 2013. Vol. 18, No. 6. P. 31-43. (In Russ.)
[11] Семин Л.Г., Слепцов А.Г., ШШапеев В.П. Метод коллокаций-наименьших квадратов для уравнений Стокса // Вычисл. технологии. 1996. Т. 1, № 2. С. 90-98.
Semin, L.G., Sleptsov, A.G., Shapeev, V.P. Collocation and least-squares method for Stokes equations // Comput. Technologies. 1996. Vol. 1, No. 2. P. 90-98. (In Russ.)
[12] Исаев В.И., ШШапеев В.П. Развитие метода коллокаций и наименьших квадратов // Труды ИММ УрО РАН. 2008. Т. 14, № 1. С. 41-60.
Isaev, V.I., Shapeev, V.P. Development of the collocations and least squares method // Trudy IMM UrO RAN. 2008. Vol. 14, No. 1. P. 41-60. (In Russ.)
[13] Шапеев В.П., Ворожцов Е.В., Исаев В.И., Идимешев С.В. Метод коллокаций и наименьших невязок для трехмерных уравнений Навье — Стокса // Вычисл. математика и программирование. 2013. Т. 14, № 1. С. 306-322.
Shapeev, V.P., Vorozhtsov, E.V., Isaev, V.I., Idimeshev, S.V. The method of collocations and least residuals for three-dimensional Navier —Stokes equations // Numerical Methods and Programming. 2013. Vol. 14, No. 1. P. 306-322. (In Russ.)
[14] Shapeev, V. Collocation and least residuals method and its applications // EPJ Web of Conf. 2016. Vol. 108. P. 01009-1-01009-12. D0I:10.1051/epjconf/201610801009.
[15] Исаев В.И., ШШапеев В.П., Еремин С.А. Исследование свойств метода коллокаций и наименьших квадратов решения краевых задач для уравнения Пуассона и уравнений Навье —Стокса // Вычисл. технологии. 2007. Т. 12, № 3. С. 53-70.
Isaev, V.I., Shapeev, V.P., Eremin, S.A. An investigation of the collocation and the least squares method for solution of boundary value problems for the Navier —Stokes and Poisson equations // Computat. Technologies. 2007. Vol. 12, No. 3. P. 53-70. (In Russ.)
[16] Ворожцов Е.В., ШШапеев В.П. Об ускорении итерационных процессов решения краевых задач комбинированием методов Крылова и Федоренко // Символ науки. 2015. № 10-2. С. 24-43.
Vorozhtsov, E.V., Shapeev, V.P. On acceleration of iterative processes for solving boundary value problems by combining Krylov and Fedorenko methods // Simvol Nauki. 2015. No. 10-2. P. 24-43. (In Russ.)
[17] Слепцов А.Г. Об ускорении сходимости линейных итераций. II // Моделирование в механике: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ВЦ; ИТПМ. 1989. Т. 3, № 5. C. 118-125. Sleptsov, A.G. On convergence acceleration of linear iterations. II // Modelirovanie v Mekhanike. 1989. Vol. 3, No. 5. P. 118-125. (In Russ.)
[18] Saad Y. Numerical methods for large eigenvalue problems. Manchester Univ. Press, 1991. 358 p.
[19] Тимошенко С.П., Войновский-Кригер С. Пластины и оболочки. М.: Физмат-гиз, 1963. 536 c.
Timoshenko, S.P., Woinowsky-Krieger, S. Theory of plates and shells. 2 ed. New York, Toronto, London: McGraw-Hill, 1959. 580 p.
Поступила в 'редакцию 21 декабря 2016 г., с доработки — 22 февраля 2017 г.
The versions of collocation and least residuals method for solving problems of mathematical physics in the trapezoidal domains
Belyaev, Vasily A.1,2'*, Shapeev, Vasily P.1,2
1 Novosibirsk State University, Novosibirsk, 630090, Russia 2Khristianovich Institute of Theoretical and Applied Mechanics SB RAS, Novosibirsk, 630090, Russia * Corresponding author: Belyaev, Vasily A., e-mail: belyaevasily@mail.ru
The paper addresses new versions of collocations and least residuals (CLR) method are proposed and implemented for the numerical solution of boundary value problems for PDE in trapezoidal domains. Their implementation and numerical experiments are
© ICT SB RAS, 2017
42
B. А. BexaeB, B. n. WlaneeB
performed on the examples of the biharmonic and Poisson equations. The solution of the biharmonic equation is used for simulation of the stress-strain state of an isotropic plate under the action of the transverse load.
Three different variants of the construction of the computational grid in the trapezoidal cells are implemented in the present study. First of all, the original domain is covered by a regular grid with trapezoidal cells. In the second variant of the method we apply the idea of using parts of the cells of a regular grid outside the domain, which are cut off by the boundary for the constructing the CLR methods. It is assumed that the solution has no singularities at the boundary and in a certain small neighborhood of it. The differential equation for the problem is true not only in the computational domain, but also in a small neighborhood of its boundary. Then we use the idea of joining the "small" irregular cells to the adjacent cells in order to reduce the condition number for the global system of linear algebraic equations.
It is shown that the approximate solutions, obtained by CLR, converge with high order of accuracy, thus accurately match the analytical solutions of test problems.
Keywords: collocations and least residuals method, boundary value problem, trapezoidal domain, high order approximation, Poisson's equation, biharmonic equation.
Acknowledgements. This research was supported by RFBR (grant No. 13-01-00227).
Received 21 December 2016 Received in revised form 22 February 2017
