CC BY
9
Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки. 2022. Т. 26, № 3. С. 556-572_
ISSN: 2310-7081 (online), 1991-8615 (print) d https://doi.org/10.14498/vsgtu1936
EDN: ETBROJ
УДК 519.635.1
hp-Вариант метода коллокации и наименьших квадратов с интегральными коллокациями решения бигармонического уравнения
В. П. Шапеев1'2, Л. С. Брындин1'2, В. А. Беляев1
1 Институт теоретической и прикладной механики им. С. А. Христиановича СО РАН, Россия, 630090, Новосибирск, ул. Институтская, 4/1.
2 Новосибирский государственный университет, Россия, 630090, Новосибирск, ул. Пирогова, 1.
Аннотация
Разработан новый алгоритм численного решения бигармонического уравнения. Он основан на впервые реализованном hp-варианте метода коллокации и наименьших квадратов (hp-МКНК) с интегральными коллокациями для эллиптического уравнения четвертого порядка в комбинации с современными способами ускорения итерационных процессов решения систем линейных алгебраических уравнений (СЛАУ). В hp-МКНК использовались его возможности измельчать шаги расчетной сетки (h-подход) и увеличивать степень базисных аппроксимирующих полиномов до произвольного порядка (p-подход). На примере численного моделирования изгиба шарнирно закрепленной изотропной пластины проведен анализ сходимости приближенных решений, полученных реализованным вариантом метода. Показано достижение высокой точности и повышенного порядка сходимости решений при применении полиномов высоких вплоть до десятой степеней в hp-МКНК.
Исследована эффективность комбинированного применения сочетающихся с МКНК алгоритмов ускорения итерационных процессов решения СЛАУ. Применены предобуславливание матриц СЛАУ; алгоритм
Математическое моделирование, численные методы и комплексы программ Научная статья
© Коллектив авторов, 2022 © СамГТУ, 2022 (составление, дизайн, макет)
3 ©® Контент публикуется на условиях лицензии Creative Commons Attribution 4.0 International (https://creativecommons.org/licenses/by/4.0/deed.ru) Образец для цитирования
Шапеев В. П., Брындин Л. С., Беляев В. А. hp-Вариант метода коллокации и наименьших квадратов с интегральными коллокациями решения бигармонического уравнения // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2022. Т. 26, № 3. С. 556-572. EDN: ETBROJ. DOI:10.14498/vsgtu1936. Сведения об авторах
Василий Павлович Шапеев © https://orcid.org/0000-0001-6761-7273 доктор физико-математических наук, профессор; главный научный сотрудник; лаб. термомеханики и прочности новых материалов1; профессор; каф. математического моделирования2; e-mail: [email protected]
Лука Сергеевич Брындин & https://orcid.org/0000-0002-0211-5800 младший научный сотрудник; лаб. термомеханики и прочности новых материалов1 ; аспирант; каф. математического моделирования2; e-mail: [email protected] Василий Алексеевич Беляев © https://orcid.org/0000-0001-5901-2880 младший научный сотрудник; лаб. термомеханики и прочности новых материалов1; e-mail: belyaevasily@mail .ru
ускорения итераций, основанный на подпространствах Крылова; операция продолжения на многосеточном комплексе; распараллеливание вычислительной программы с помощью ОрепМР; модифицированный алгоритм решения локальных СЛАУ, определяющих решение задачи в каждой ячейке сетки. Последний, применимый в случае решения линейного дифференциального уравнения, позволяет более эффективно решать переопределенные СЛАУ в МКНК, реализуемом итерациями по подобластям, в которых вид матриц локальных СЛАУ не изменяется на каждой итерации. Комбинированное применение всех перечисленных способов ускорения уменьшило время расчетов на персональном компьютере более чем в 350 раз по сравнению со случаем, когда использовалось только предобуславливание.
Ключевые слова: метод коллокации и наименьших квадратов, интегральные уравнения коллокации, бигармоническое уравнение, изгиб пластины, ускорение итерационных процессов, предобуславливание, подпространства Крылова, многосеточный алгоритм, распараллеливание.
Получение: 21 июня 2022 г. / Исправление: 9 сентября 2022 г. / Принятие: 13 сентября 2022 г. / Публикация онлайн: 29 сентября 2022 г.
Введение. Известно, что численное решение краевых задач для уравнений бигармонического типа, встречающихся в различных разделах механики сплошных сред, например, в теории тонких пластин [1,2], в гидродинамике при малых числах Рейнольдса [3] и др., вызывает ряд трудностей. Это связано с наличием в эллиптическом уравнении производных четвертого порядка, существенно сказывающихся на обусловленности исходной дифференциальной задачи. При этом цель многих исследователей заключается в разработке новых вариантов высокоточных численных методов, по развитию которых в настоящее время ведется активная работа [3—11]. Одним из таких перспективных методов является проекционно-сеточный метод коллокации и наименьших квадратов (МКНК) [9—11]. В [11] при решении бигармонического уравнения в нерегулярных областях применен вариант МКНК с «дифференциальными» коллокациями, аппроксимирующими уравнение в точках внутри ячеек. Там же продемонстрированы преимущества и достоинства предложенных подходов при сравнении полученных МКНК результатов с результатами других авторов, использовавших метод конечных разностей (МКР) [3], метод конечных элементов (МКЭ) [8] и спектральные методы [5-7].
Под Ир-вариантом МКНК (Ир-МКНК) подразумевается такой вариант метода, который позволяет как измельчать шаги расчетной сетки (И-подход), так и увеличивать степени аппроксимирующих полиномов (р-подход). Понятие введено по аналогии с Ир-МКЭ [12].
Данная работа посвящена реализации и верификации Ир-МКНК с «интегральными коллокациями», аппроксимирующими закон сохранения, записанный в интегральном виде, впервые примененного к решению бигармониче-ского уравнения. Область влияния уравнения интегральной коллокации вся ячейка, по площади которой проводится интегрирование при его записи. Область влияния дифференциальной коллокации — малая окрестность точки, в которой она аппроксимирует дифференциальное уравнение задачи. В [10] бигармоническое уравнение было решено МКНК с применением интеграль-
ных уравнений коллокации, взвешенных наименьших квадратов и полиномов только четвертой степени. Отмечалось, что в этом случае наблюдается сходимость МКНК в более широком диапазоне значений весовых множителей по сравнению с дифференциальным вариантом МКНК, в котором их приходится подбирать более тщательно [9]. В МКНК на весовые множители до-множаются уравнения коллокации, условия согласования и краевые условия, а их объединение составляет переопределенную систему линейных алгебраических уравнений (СЛАУ) в каждой ячейке расчетной области. Глобальная СЛАУ, являющаяся совокупностью локальных СЛАУ, в МКНК чаще всего решается с помощью метода итераций по подобластям в комбинации с современными способами ускорения итерационных процессов.
Сочетание МКНК с предобуславливанием (диагональным [13] и многопараметрическим, связанным с выбором значений весовых множителей), ускорением, основанным на методе подпространств Крылова [14], и операцией продолжения, используемой для задания начального приближения на многосеточном комплексе [15] при переходе с грубой сетки на более подробную, показало уменьшение времени расчетов до 362 раз по сравнению со случаем, когда ни один из алгоритмов не использовался на примере решения уравнений Навье—Стокса [16]. При решении краевой задачи Дирихле для уравнения Пуассона в [17] продемонстрировано сокращение времени расчетов практически в 2000 раз и количества итераций более чем в 250 раз, где помимо вышеуказанных методов также применялось распараллеливание вычислительной программы с помощью ОрепМР в сочетании с обходом области, основанным на красно-черном упорядочивании [18], а также свойство локальной системы координат в МКНК, которое имеет место в случае решения линейного дифференциального уравнения.
Однако комплексное исследование влияния алгоритмов ускорения для решения уравнений с частными производными (УЧП) четвертого порядка с использованием полиномов различных степеней в МКНК пока не проводилось, чему отчасти и посвящена настоящая работа. При этом стоит отметить, что ускорение по Крылову в некоторых ситуациях способно сделать расходящийся итерационный процесс сходящимся [19], что является особо важным при решении плохо обусловленных задач.
Подчеркнем, что тематика ускорения итерационных процессов сама по себе является весьма актуальной и активно развивается. В качестве примера здесь упомянем работу [20], где экспериментально изучается влияние размеров пересечений (перехлеста) подобластей на скорость сходимости методов декомпозиции. Также в [21] можно познакомиться с одним из примеров двухуровневых предобуславливателей для методов Крылова — методом дефляции. В этой же работе [21] проводится достаточно подробное сравнение между методами дефляции, декомпозиции областей и многосеточных.
1. Постановка задачи. Рассмотрим задачу статического изгиба изотропной прямоугольной пластины О С К2 с границей 5О в рамках теории Кирхгофа—Лява. В этом случае прогиб пластины -ш(ж1,ж2) определяется из решения бигармонического уравнения (уравнения Софи Жермен—Лагран-жа) [1]:
дх\ дх2дх2 дх2 О'
где q(x\,x2) — поперечная нагрузка, D = Et3/[12(1 — г/2)] —жесткость пластины при изгибе, t — толщина пластины, E — модуль Юнга, v — коэффициент Пуассона.
Пусть пластина шарнирно закреплена, тогда краевые условия имеют следующий вид [1]:
w = 0, Мп = 0, (2)
где Мп — изгибающий момент, который определяется по формуле
Мп = МХ1 cos2 a + МХ2 sin2 a — 2МХ1Х2 sin a cos a,
где
/d2w d 2w \ _ /d2w d2w \ _ . d2w
М- =—n^2 +г/щ),М- = +ищ), М- = D(1—);
(cosa, sin a) —компоненты внешней единичной нормали к <Ш.
Краевая задача (1), (2) плохо обусловлена, так как малые изменения во входных данных задачи приводят к значительным изменениям ее решения. Необходимо получить эффективное численное решение данной плохо обусловленной краевой задачи. Для достижения этой цели авторами были поставлены и решены следующие задачи:
1) для получения надежных результатов и более устойчивого счета были использованы интегральные уравнения коллокации;
2) для эффективного решения соответствующей краевой задачи для би-гармонического уравнения впервые был реализован и верифицирован hp-МКНК с интегральными уравнениями коллокации в комбинации с современными алгоритмами ускорения итерационных процессов;
3) было проведено комплексное исследование влияние алгоритмов ускорения итерационных процессов в МКНК при решении дифференциального уравнения четвертого порядка.
2. Интегральный закон сохранения. Домножим (1) на D, проинтегрируем левую и правую часть дважды с использованием формулы Грина и получим интегральный закон сохранения для бигармонического уравнения
j) F2dx2 — F2dx2 = jj q(x2,x2)dx2dx2, (3)
где
v
^ ( !, 2) ídx3 + dxtdx2, dx'i, + dx2dx<^) ( ' ^x2^
. дх\ дх\дх2 дх3 дх2дх
V С О, £ — граница подобласти V С О, QX1 и QX2 — перерезывающие силы. Уравнение (3) имеет определенный физический смысл: действие внешней поперечной нагрузки ц на элемент пластины V уравновешивается перерезывающими силами, действующими на боковых площадках этого элемента. Таким образом, отыскание обобщенных решений задачи (1), (2) свелось к решению краевой задачи (3), (2).
Отметим, что подобный подход решения дифференциального уравнения, исходно заданного в дивергентном виде, эффективно применялся при решении задач с разрывными параметрами при построении различных консервативных численных методов с целью нахождения обобщенных решений [22].
Замечание. Непосредственно можно аппроксимировать как исходное би-гармоническое уравнение (1) [5-7,9,11], так и систему УЧП, к которой оно может сводиться. В некоторых случаях это сделать достаточно просто, например, когда рассматривается изгиб шарнирно закрепленной пластины полигональной формы [1]. В случае сложных краевых условий и нерегулярной области это может вызывать существенные трудности. Но зачастую решать задачу для системы УЧП более низкого порядка проще из-за ее лучшей обусловленности по сравнению с задачей для уравнения высокого порядка. Например, в [3, 4] вводится дополнительная неизвестная переменная и решение исходной задачи Дирихле (на границе заданы значения и> и производной по нормали дтл/дп) для бигармонического уравнения сводится к последовательному решению двух уравнений Пуассона.
3. Описание hp-МКНК. Покроем О регулярной сеткой, состоящей из N1 х N2 прямоугольных ячеек. Обозначим через ^сеПв = N^2 общее количество ячеек в расчетной области. Введем локальную систему координат в каждой 8-й ячейке, 5 = 1, 2,..., ^ец8, с линейной зависимостью от глобальной:
У1 = (Х1 -Х1э)/к1, У2 = (Х2 -Х2з)/Ь2,
где (х1з,х2з) — координаты центра ячейки в глобальной системе, 2^1х2^2 — размер ячейки в направлении X] и Х2 соответственно.
Приближенное решение задачи в каждой ячейке ищем в виде
к к—1
унз(Уъ у 2) = ^ У!У22. (4)
¿1=0 г2=0
Неизвестные коэффициенты Сг1г2,8 определим из решения переопределенной СЛАУ, полученной аппроксимацией задачи (3), (2) полиномами (4) и состоящей из уравнений коллокации, условий согласования и краевых условий, если ячейка является граничной. Обозначим Уь(у1(х1 ),у2(х2)) = ,ш(х1,х2).
Определим потоки вектора Р через стороны прямоугольных ячеек сетки (ъ, з), ъ = 1, 2,... ,N1, ] = 1, 2,..., N2, вычисленные на приближенном решении ьь,:
±
шо Г"' +
'Х1-,
V дх2 дх2дх2
Х1=Х1
№
±
шв Г'+' (^ + ^
1X2,-
V дх\ дх1дх2.
Х2=Х2,]
с1х2.
(5)
(6)
Здесь значения с полуцелыми индексами обозначают потоки сохраняемых величин через границы расчетной ячейки.
С учетом (5), (6) закон сохранения (3), записанный для ячейки с номером (г, ]), имеет вид
Г%1,1+1 Г х2,] + 1
К+1л+1/2-^3+1/2 + К~+1/2,3+1-]¥Ш/2,3 = / Ч(х1, х2)^х1^х2. (7)
Предлагается каждую прямоугольную ячейку сетки разделить на К2 одинаковых прямоугольных ячеек и записать в них уравнения коллокации (7) в локальных координатах (у 1,^2). При этом уравнения (7) домножим на весовой множитель рс.
На каждой стороне 8-й прямоугольной ячейки равномерно расставим К точек согласования, в которых по аналогии с [9-11] запишем условия согласования:
дь ы „. дун
Рто ЪНз + Рт1- = Рто Ьн + Ртг^~, (8)
дп8 дп8
д2 V из + д 3у из = д 2ун + д Чн (9)
Рт2 дп2 + Ртз дп3 = Рт2 дп2 + Ртз дп\ , (9)
где п8 — внешняя единичная нормаль к соответствующей границе 8-й ячейки; Ъиз и "Ьн — пределы приближенного решения задачи при стремлении изнутри и извне к границе в-й ячейки; Рто, Рт1, Рт2, Ртз — положительные весовые множители.
Краевые условия для шарнирно закрепленного края (2) записываются на прямолинейной стороне ячейки, являющейся частью 5О, в К равномерно распределенных на ней точках:
Рь0Укз = о, Рь2Мп = 0, (10)
где Рь0, Рь2 —положительные весовые множители.
Объединяя уравнения (7)-(10) в каждой ячейке относительно неизвестных Сг1г2,8 получим локальную переопределенную СЛАУ
Ах = Ь, (11)
где А — прямоугольная вещественная матрица размера тх I, х — вектор неизвестных, Ь — вектор правой части. При описанном выше способе записи уравнений имеем относительно I = ( К + 1)(К + 2)/2 неизвестных в каждой ячейке т = К2 + 4К линейных алгебраических уравнений.
Глобальная СЛАУ, являющаяся объединением локальных СЛАУ, здесь решается с помощью метода итераций по подобластям. Одна глобальная итерация заключается в последовательном решении локальных СЛАУ (11) в каждой ячейке области. Для решения линейной задачи наименьших квадратов использовалась Q_R-декомпозиция матрицы А, найденная методом отражений Хаусхолдера [23]. Итерационный процесс продолжался до тех пор, пока не было выполнено следующее условие его остановки:
тах 1- ^2,31 < £, 1 2 1 2 1 2
где — номер итерации, — наперед заданная малая константа, называемая псевдопогрешностью решения. Значение зависит от решаемой задачи, размера ячейки сетки 2к]х2к2, старшей степени полиномов К в (4). В МКНК она выбирается экспериментально так, чтобы погрешность решения глобальной СЛАУ была по возможности существенно меньше ожидаемой погрешности решения, исходя из значений погрешностей на более грубых сетках. В приведенной ниже таблице значения е варьировались от 10-13 до 10-14.
Следует отметить, что от выбора значений весовых множителей могут сильно зависеть скорость сходимости итерационного процесса, обусловленность СЛАУ и точность приближенного решения [9,20,24].
4. Эффективное использование свойства матрицы локальной
СЛАУ. В МКНК при решении линейной задачи вид матрицы локальной СЛАУ в каждой ячейке одинаковый и не зависит от номера итерации. Также для каждой ячейки в итерационном процессе компоненты вектора правой части Ь, соответствующие уравнениям коллокации и краевым условиям (если ячейка граничная), не изменяются.
Пусть Ь — вектор размера т, первые компоненты которого в количестве взятых условий согласования совпадают с компонентами вектора , являющимися правыми частями уравнений согласования. Остальные компоненты
вектора Ь являются нулями. Обозначим Ь = Ь — Ь, 0Т —транспонированная матрица размера 1хт для матрицы 0, К-1 — обратная для верхнетреугольной матрицы К размера 1х1 для каждой локальной матрицы А. В связи
с этим достаточно только один раз вычислить и запомнить К-10Т и К-10ТЪ
и на каждой итерации для каждой ячейки положить х = К-10ТЪ + К-10ТЬ, т.к. в итерационном процессе изменяется только правая часть условий согласования. Отметим также, что при решении линейной задачи с постоянными коэффициентами для внутренних прямоугольных ячеек вид матриц К-1 и 0Т одинаковый, поэтому их можно вычислить один раз и записать в памяти компьютера по одному экземпляру для всех таких ячеек, чтобы сэкономить память.
Ранее, например, при решении уравнения Пуассона МКНК [17] предлагалось на каждой итерации умножать на записанную в памяти компьютера матрицу 0Т, а затем делать обратный ход, как и в методе Гаусса. Очевидно, что реализованная здесь модификация алгоритма решения СЛАУ еще сильнее уменьшает время расчетов по сравнению с методикой, описанной в [17].
5. Результаты расчетов и их обсуждение. Рассмотрим изгиб шар-нирно закрепленной (2) пластины размера (-хг^ м, на которую действует синусоидальная нагрузка д = 105 ё'т(ттх1/(11) ё'т(ттх2/(12) Па. В этом случае известно точное решение задачи (1), (2), см., например, [1]:
■(х1,х2) =
1(2
■ж40((12 + (11)2'
В численных расчетах полагалось = (2 = 10 м, £ = 0.1 м, Е = 200 ГПа, V = 0.28 и ^хЩ = NхN, Н1 = к2 = к.
4
В таблице приведены результаты численных экспериментов, где относительная погрешность приближенного решения определялась по формуле
max ( max lw(xlm,x2m) - vhs(xlm,X2m)\)
lip II _ s=1,...,Ncells m=1,...,Qs_
r ^ max ( max \w(xim,x2m)\) ,
S=1,...,Ncells m=1,...,Qs
где Qs _ 100 — количество равномерно распределенных контрольных точек (xim,X2m) для подсчета в них погрешности. Порядок сходимости погрешности приближенного решения определялся величиной
\og2(EN/2/EN),
где En/2 — погрешность \\Erна сетке, состоящей из N/2xN/2 ячеек, En — погрешность \\Er на сетке, состоящей из NxN ячеек.
В таблице приведено количество итераций Niter, время выполнения итерационного процесса в секундах tso\ и их отношения AFj и AF соответственно для случая «с ускорением» для итерационного процесса, в котором комбинированно использовались:
- диагональный предобуславливатель,
- Рс = h2, Pm0 = Pm1 = Pb0 = 1, Pm2 = Pm3 = Pb2 = h2,
- ускорение по Крылову,
- операция продолжения,
- распараллеливание с помощью OpenMP с обходом области, основанным на красно-черном упорядочивании,
- модифицированный алгоритм решения локальных СЛАУ
по отношению к случаю «без ускорения», когда применялся только диагональный предобуславливатель, фиксировались аналогичные значения весовых множителей и был реализован последовательный обход области.
Вычисления проводились на компьютере Intel Core i5-8265UCPU 1.6 GHz, 4Cores, DIMM DDR4-2400 1200MHz 8Gb. Для распараллеливания с помощью OpenMP использовалось 4 потока.
В расчетах в итерационном процессе на самой грубой сетке (и во всех случаях без использования операции продолжения) в начальном приближении решения взяты Ci1i2,j = 0.4. Для демонстрации представления чисел здесь используется экспоненциальный формат их записи.
Проведенное исследование показало следующее.
1. Порядок сходимости в среднем не хуже 0(hK-p+2-Kmod2), что хорошо согласуется со многими результатами решения различных задач МКНК [10,17,25], где р — порядок уравнения (в данном случае 4).
2. Наличие сходимости МКНК для достаточно высоких степеней полиномов К, что позволяет в случае достаточно гладких решений получать высокоточное решение с малым количеством степеней свободы (DOF, от англ. degrees of freedom), например, чтобы достичь точность порядка 10-10, необходимо DOF = 57600 при К = 7, DOF = 18000 при К = 8, DOF = 14080 при К = 9, наконец, DOF = 4224 при К = 10. В МКНК DOF — общее количество неизвестных коэффициентов в представлении приближенного решения (4).
<1
<1
£
£
£
s;
м о
£
X £
оо со 4
г- г-
оо ^
^о ^ ю оо со оо оо
сч ^ о со
CN О CN t^ © НП CS О CN Ю CN О О О СО
оо оо CN ю
оо ^
о оо оо
оо оо оо
НСО N
Í—1 CN 6
oo cn i—i
со 00
4 6 oo o
H H H H
5 CO oo oo
í—1 í—1 t^ o¿
ю o 1—1 o o
X X X X
ю o 1—1 o o
5
оо iv
46
со к ^
со оо со оо оо со ^ оо
N ^f CS Сй
оо со t^ оо
IV СО ю о
О CN СО О
о о о IV
^ ^ ^ оо
16 51
51
оо CN ю оо
CN 00 ^н о оо о ^оо оо
oo í—1 oo
IV CD 00
co Ю
H H H н
Ю ^ о ^
CN CS О OO
00 СЧ О ^
o o o
ю 1—1
X X X X
Ю o o o
1—1
СО СО IV toiMra
5
оо о ю i—i оо
оо :
о оо оо
CO N Ю
о
CN
о
CN
ю оо оо о
M СЧ Ю (О О О ^н Ю
444
с^ ^ ^
I««
СО о
г- о
CN CN
^ оо CN оооо и со оо
oo 4 o
oo CN СО
Ю со со
H H н н о- о- о о
гН н N
oó ю
o o o
Ю 1—1
X X X X
Ю o o o
1—1
5
¡с оо «э
IV
о
CN
COCOIM «
1—1 О СО о
оо о о о
444
1-1 ^ ю œ со N О! о ^ ^
ОЗ СО
4
66 со оо i—i оо
œt-N СО IV ^н
H H H н
со IV ю о
©Ю M tJ
CN oó oó ^
СО
С^ ^ 00 ^ч X X X X с^ ^ оо со
IV £
со оо IV
^Н ^ ОЗ
о IV IV оо
оо сч CN со IV ю оо о ^^ IV
í— О i— i—I
444
IV со ^ 2
^ О СО .
о ^ ^ ю о
4
оо ; оо ^
о IV
4
IV СО
64 i—i оо
í—1 4
p
o í—1
í—1 í—1
H H H H ю to f-oo o es ю
^H CO CN t^
CO
C^ ^ 00
X X X X
c^ ^ оо со
3. Сходимость численного решения без применения алгоритмов ускорения итерационного процесса решения СЛАУ, что говорит о хороших свойствах аппроксимации в МКНК при решении УЧП высокого порядка с применением полиномов высоких степеней и использованием интегральных уравнений коллокации.
4. Заметное сокращение времени решения задачи (более 350 раз) и количества итераций (более 20 раз) в случае комбинирования МКНК с другими вычислительными алгоритмами. Коэффициенты ускорения А^ и А^ увеличиваются по мере роста К и уменьшения шагов сетки при решении УЧП четвертого порядка.
5. Применение впервые реализованного в данной работе модифицированного алгоритма решения локальных СЛАУ в МКНК позволило сократить время расчетов в 7.42 раза при К = 6, в 7.93 раза при К = 7, в 9.77 раз при К = 8, в 11.73 раза при К = 9 ив 12.41 раз при К = 10 относительно времен, указанных в пятой колонке таблицы. Данные величины были рассчитаны как среднее арифметическое от А^, вычисленного на четырех указанных сетках в таблице. Как и следовало ожидать, этот эффект возрастает при увеличении степени полиномов К.
Исходя из своего прежнего опыта решения бигармонического уравнения для эффективного построения поправок приближенных решений по Крылову авторы использовали количество невязок N =40. Это оказалось достаточным для достижения заданной точности, поэтому практически всегда = 41 в таблице. Однако для каждого случая оптимальное значение величины N может быть другим. Например, если использовать N = 50 на сетке 40х40 при К = 7, то = 51, а = 4.921 с, т.е. в этом случае коэффициенты ускорения станут равными А^ = 15.70 и АЕ4 = 269.07. Максимальные числа обусловленности матриц локальных СЛАУ увеличивались от С- ■ 102 до С2 ■ 105 при увеличении К и уменьшении шагов сетки в приведенных в таблице случаях, где С- и С2 — константы.
Сравнение с опубликованным ранее дифференциальным Ир-МКНК для этой задачи [11], когда уравнения коллокации записывались в равномерно распределенных точках коллокации и имели вид
(1 д4Уьз д4уН8 д4уНз \ =
Рс\к4 ду4 + к4ду2ду\ + к4 ду4 ) = РсО'
с применением базиса (4) и рс = к4 (остальные весовые множители аналогичны тем, которые использовались в интегральном варианте, если не оговорено иное) показало следующее.
1. Итерационный процесс для дифференциального МКНК разошелся при К = 4 по сравнению с интегральным МКНК без использования ускорения по Крылову на сетках 40х40, 80х80 при указанных выше весовых множителях и на сетке 80х80, когда рт1 = 0.12, ртз = 0.125к2.
2. Точность при применении дифференциального МКНК упала по сравнению с интегральным, результаты которого приведены в таблице. Дифференциальный вариант МКНК позволил достичь точности \\ЕГ= = 1.07Е-9 при К = 7 на сетке 40х40, а интегральный— \\ЕГ= = 4.34Е-10. При использовании полиномов более высоких степеней наблюдалась аналогичная картина. Например, с помощью дифферен-
циального МКНК при К = 8 на сетке 40x40 получено значение погрешности \\Er IU = 3.82 E-11, а при К = 10 на сетке 16x16 — \\Er = = 1.36 E-12. Применение интегрального МКНК позволяет достичь WEr= 7.29E-12 и W\Er= 7.58E-13 соответственно при К = 8 и К = 10.
Стоит отметить, что многие результаты у двух сравниваемых здесь вариантов МКНК оказались близкими, что говорит об их хороших аппроксимаци-онных свойствах. Однако, как это показано выше, применение интегрального варианта позволяет улучшить свойства метода, что является особо принципиальным при решении плохо обусловленных задач.
Заключение. В данной работе на примере решения бигармонического уравнения в приложении к моделированию изгиба тонкой пластины показана возможность достижения высокой точности hp-МКНК с интегральными кол-локациями и проведена его верификация с использованием полиномов вплоть до десятой степени. Продемонстрировано существенное сокращение количества итераций и времени расчетов при использовании различных способов ускорения итерационных процессов в МКНК.
Как видно из описанного алгоритма, hp-МКНК с представлением решения в виде (4) привлекателен прежде всего своей возможностью автоматической реализации на компьютере с точки зрения применения произвольных степеней полиномов для достижения повышенной точности. Для других численных методов это может вызывать более существенные трудности. Например, в МКР в [26] для этой цели авторы дополнительно используют систему компьютерной алгебры.
Таким образом, показанная здесь возможность сокращения времени расчетов в несколько сотен раз и высокоточного решения плохо обусловленной задачи открывает все более широкие перспективы эффективного численного моделирования напряженно-деформированного состояния тонких пластин с помощью МКНК, в том числе нерегулярных форм и с особенностями [11]. Разрабатываемый комплекс программ со временем может быть успешно применен для моделирования и оптимизации композитных пластин достаточно произвольных форм в рамках различных теорий пластин [27,28].
Конкурирующие интересы. Заявляем, что в отношении авторства и публикации этой статьи конфликта интересов не имеем.
Авторский вклад и ответственность. Все авторы принимали участие в разработке концепции статьи и в написании рукописи. Авторы несут полную ответственность за предоставление окончательной рукописи в печать. Окончательная версия рукописи была одобрена всеми авторами.
Финансирование. Институт теоретической и прикладной механики им. С. А. Хри-стиановича СО РАН. Работа выполнена в рамках государственного задания (№ госрегистрации 121030500137-5).
Библиографический список
1. Timoshenko S., Woinowsky-Krieger S. Theory of Plates and Shells. New York, St. Louis:
McGraw-Hill Book Comp., 1959. xiv+580 pp.
2. Reddy J. N. Mechanics of Laminated Composite Plates and Shells. Theory and Analysis.
Boca Raton: CRC Press, 2011. xxiv+832 pp. DOI: https://doi.org/10.1201/b12409.
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., 2008. vol.29, no. 2. pp. 113-133. DOI: https://doi.org/10.1007/s10444-007-9043-6.
4. McLaurin G. W. A general coupled equation approach for solving the biharmonic boundary value problem// SIAM J. Numer. Anal., 1974. vol. 11, no. 1. pp. 14-33. DOI: https://doi. org/10.1137/0711003.
5. 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. vol. 33, no. 1. pp. 284-299. DOI: https://doi.Org/10.1016/j.apm.2007.11.002.
6. 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. vol. 36, no. 12. pp. 1787-1798. DOI: https://doi.org/10.1016/j.enganabound.2012. 06.005.
7. 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. vol.39, no.9. pp. 2554-2569. DOI: https://doi.org/10.1016/j.apm.2014.10.048.
8. Guo H., Zhang Z., Zou Q. A С0 linear finite element method for biharmonic problems// J. Sci. Comput., 2018. vol.74, no. 3. pp. 1397-1422. DOI: https://doi.org/10. 1007/s10915-017-0501-0.
9. Голушко С. К., Идимешев С. В., Шапеев В. П. Метод коллокаций и наименьших невязок в приложении к задачам механики изотропных пластин // Вычисл. технол., 2013. Т. 18, №6. С. 31-43. EDN: RUDTZF.
10. Belyaev V. A., Bryndin L. S., Shapeev V. P., Golushko S. K. The least squares collocation method with the integral form of collocation equations for bending analysis of isotropic and orthotropic thin plates// AIP Conf. Proc., 2020. vol.2288, 030084. EDN: FWKXNV. DOI: https://doi.org/10.1063/5.0028298.
11. Беляев В. А., Брындин Л. С., Голушко С. К., Семисалов Б. В., Шапеев В. П. h-, p- и hp-варианты метода коллокации и наименьших квадратов для решения краевых задач для бигармонического уравнения в нерегулярных областях и их приложения // Ж. вычисл. матем. и матем. физ., 2022. Т. 62, №4. С. 531-552. EDN: MSMVOK. DOI: https://doi.org/ 10.31857/S0044466922040020.
12. Babuska I., Guo B. Q. The h, p and h-p version of the finite element method; basis theory and applications// Adv. Eng. Softw., 1992. vol. 15, no. 3-4. pp. 159-174. DOI: https://doi. org/10.1016/0965-9978(92)90097-Y.
13. Ramsak M., Skerget L. A subdomain boundary element method for high-Reynolds laminar flow using stream function-vorticity formulation// Int. J. Numer. Meth. Fluids, 2004. vol. 46, no. 8. pp. 815-847. DOI: https://doi.org/10.1002/fld.776.
14. Saad Y. Numerical Methods for Large Eigenvalue Problems. Philadelphia: SIAM, 2011. xvi+276 pp. DOI: https://doi.org/10.1137Z1.9781611970739.
15. Федоренко Р. П. Введение в вычислительную физику. М.: МФТИ, 1994. 528 с. EDN: QJUAEP.
16. Ворожцов Е. В., Шапеев В. П. Бездивергентный метод коллокаций и наименьших квадратов для расчета течений несжимаемой жидкости и его эффективная реализация // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2020. Т. 24, №3. С. 542-573. EDN: OSPEBL. DOI: https://doi.org/10.14498/vsgtu1758.
17. Беляев В. А. Об эффективной реализации и возможностях метода коллокации и наименьших квадратов решения эллиптических уравнений второго порядка // Выч. мет. программирование, 2021. Т.22, №3. С. 211-229. EDN: VBQEPE. DOI: https://doi.org/ 10.26089/NumMet.v22r313.
18. Ortega J. M. Introduction to Parallel and Vector Solution of Linear Systems. New York: Springer, 1988. xi+305 pp. DOI: https://doi.org/10.1007/978-1-4899-2112-3.
19. Karamyshev V. B., Kovenya V. M., Sleptsov A. G., Cherny S. G. Variational method of accelerating linear iterations and its applications// Comput. Fluids., 1996. vol.25, no. 5. pp. 467-484. EDN: LDTBIH. DOI: https://doi.org/10.1016/0045-7930(96)00011-4.
20. Ильин В. П., Кныш Д. В. Параллельные методы декомпозиции в пространствах следов// Выч. мет. программирование, 2011. Т. 12, №1. С. 110-119. EDN: OJAYQV.
21. Tang J. M., Nabben R., Vuik C., Erlangga Y. A. Comparison of two-level preconditioners derived from deflation, domain decomposition and multigrid methods // J. Sci. Comput., 2009. vol.39, no. 3. pp. 340-370. DOI: https://doi.org/10.1007/s10915-009-9272-6.
22. Годунов С. К., Рябенький В. С. Разностные схемы. Введение в теорию. М.: Наука, 1977. 440 с.
23. Demmel J. W. Applied Numerical Linear Algebra. Philadelphia: SIAM, 1997. xi+416 pp. DOI: https://doi.Org/10.1137/1.9781611971446.
24. Sleptsov A. G. Grid-projection solution of an elliptic problem for an irregular grid // Russ. J. Numer. Anal. Math. Model., 1993. vol.8, no. 6. pp. 519-543. DOI: https://doi.org/ 10.1515/rnam.1993.8.6.519.
25. Исаев В. И., Шапеев В. П. Варианты метода коллокаций и наименьших квадратов повышенной точности для численного решения уравнений Навье-Стокса // Ж. вычисл. матем. и матем. физ, 2010. Т. 50, №10. С. 1758-1770. EDN: MVSONV.
26. Drozdov G. M., Shapeev V. P. CAS application to the construction of high-order difference schemes for solving Poisson equation / Computer Algebra in Scientific Computing/ Lecture Notes in Computer Science, 8660; eds. V. P. Gerdt, W. Koepf, W. M. Seiler, E. V. Vorozhtsov. Cham: Springer, 2014. pp. 99-110. EDN: UEYEZH. DOI: https://doi.org/ 10.1007/978-3-319-10515-4_8.
27. Голушко С. К., Немировский Ю. В. Прямые и обратные задачи механики упругих композитных пластин и оболочек вращения. М.: Физматлит, 2008. 432 с. EDN: MUWRZD.
28. Идимешев С. В. Модифицированный метод коллокаций и наименьших невязок и его приложение в механике многослойных композитных балок и пластин: Дис. ... канд. физ.-мат. наук: 05.13.18. Новосибирск, 2016. 179 с. EDN: VPHODY.
Vestn. Samar. Gos. Tekhn. Univ., Ser. Fiz.-Mat. Nauki
[J. Samara State Tech. Univ., Ser. Phys. Math. Sci.], 2022, vol. 26, no. 3, pp. 556-572
d https://doi.org/10.14498/vsgtu1936
ISSN: 2310-7081 (online), 1991-8615 (print)
MSC: 65N35
The hp-version of the least-squares collocation method with integral collocation for solving a biharmonic equation
V. P. Shapeev1'2, L. S. Bryndin1'2, V. A. Belyaev1
1 Khristianovich Institute of Theoretical and Applied Mechanics, Siberian Branch of the Russian Academy of Sciences,
4/1, Institutskaya st., Novosibirsk, 630090, Russian Federation.
2 Novosibirsk State University,
1, Pirogova st., Novosibirsk, 630090, Russian Federation.
Abstract
A new algorithm for the numerical solution of the biharmonic equation is developed. It is based on the first implemented hp-version of the least-squares collocation method (hp-LSCM) with integral collocations for a fourth-order elliptic equation in combination with modern methods of accelerating iterative processes for solving systems of linear algebraic equations (SLAE). The hp-LSCM provides the possibilities to refine the grid (h-version) and increase the degree of polynomials to the arbitrary order (p-approach). The convergence of approximate solutions obtained by the implemented version of the method is analyzed using an example of a numerical simulation of the bending of a hinged isotropic plate. The high accuracy and the increased order of convergence using polynomials up to the tenth order in the hp-LSCM are shown.
The effectiveness of the combined application of algorithms for accelerating iterative processes to solve SLAE that are combined with LSCM is investigated. In this paper, we consider the application of the following algorithms: preconditioning of SLAE matrices; the iteration acceleration algorithm based on Krylov subspaces; the prolongation operation on a multigrid complex; parallelization using OpenMP; a modified algorithm for solving local
Mathematical Modeling, Numerical Methods and Software Complexes Research Article
© Authors, 2022
© Samara State Technical University, 2022 (Compilation, Design, and Layout)
3 ©® The content is published under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/) Please cite this article in press as:
Shapeev V. P., Bryndin L. S., Belyaev V. A. The hp-version of the least-squares collocation method with integral collocation for solving a biharmonic equation, Vestn. Samar. Gos. Tekhn. Univ., Ser. Fiz.-Mat. Nauki [J. Samara State Tech. Univ., Ser. Phys. Math. Sci.], 2022, vol. 26, no. 3, pp. 556-572. EDN: ETBROJ. DOI: 10.14498/vsgtu1936 (In Russian). Authors' Details:
Vasily P. Shapeev © https://orcid.org/0000-0001-6761-7273
Dr. Phys. & Math. Sci., Professor; Chief Researcher; Lab. of Thermomechanics and Strength of New Materials1; Professor; Dept. of Mathematical Modeling2; e-mail: [email protected] Luka S. Bryndin https://orcid.org/0000-0002-0211-5800
Junior Research Scientist; Lab. of Thermomechanics and Strength of New Materials1; Postgraduate Student; Dept. of Mathematical Modeling2; e-mail: [email protected] Vasily A. Belyaev © https://orcid.org/0000-0001-5901-2880
Junior Research Scientist; Lab. of Thermomechanics and Strength of New Materials1; e-mail: belyaevasily@mail .ru
SLAEs. The latter is implemented with iterations over subdomains (which are cells) and makes it possible to more effectively solve overdetermined SLAEs in the LSCM in the case of solving a linear differential equation. The form of the matrices does not change at each iteration. Only the elements of the vectors of their right parts corresponding to the matching conditions are modified. The calculation time on a personal computer is reduced by more than 350 times with the combined use of all acceleration techniques compared to the case when only preconditioning was used.
Keywords: least-squares collocation method, integral collocation equation, biharmonic equation, plate bending, acceleration of iterative processes, preconditioning, Krylov subspaces, multigrid algorithms, parallelization.
Received: 21st June, 2022 / Revised: 9th September, 2022 / Accepted: 13th September, 2022 / First online: 29th September, 2022
Competing interests. We declare that we have no conflicts of interest with respect to the authorship and publication of this article.
Authors' contributions and responsibilities. Each author has participated in the development of the concept of the article and in the writing of the manuscript. The authors are absolutely responsible for submitting the final manuscript in print. Each author has approved the final version of the manuscript.
Funding. Khristianovich Institute of Theoretical and Applied Mechanics, Siberian Branch of the Russian Academy of Sciences. The research was carried out within the state assignment of Ministry of Science and Higher Education of the Russian Federation (project no. 121030500137-5).
References
1. Timoshenko S., Woinowsky-Krieger S. Theory of Plates and Shells. New York, St. Louis, McGraw-Hill Book Comp., 1959, xiv+580 pp.
2. Reddy J. N. Mechanics of Laminated Composite Plates and Shells. Theory and Analysis. Boca Raton, CRC Press, 2011, xxiv+832 pp. DOI: https://doi.org/10.1201/b12409.
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., 2008, vol.29, no. 2, pp. 113-133. DOI: https://doi.org/10.1007/s10444-007-9043-6.
4. McLaurin G. W. A general coupled equation approach for solving the biharmonic boundary value problem, SIAM J. Numer. Anal., 1974, vol.11, no. 1, pp. 14-33. DOI: https://doi. org/10.1137/0711003.
5. 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, vol. 33, no. 1, pp. 284-299. DOI: https://doi.org/10.1016Zj.apm.2007.11.002.
6. 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, vol. 36, no. 12, pp. 1787-1798. DOI: https://doi.org/10.1016/j.enganabound.2012. 06.005.
7. 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, vol.39, no.9, pp. 2554-2569. DOI: https://doi.org/10.1016/j.apm.2014.10.048.
8. Guo H., Zhang Z., Zou Q. A C0 linear finite element method for biharmonic problems, J. Sci. Comput., 2018, vol.74, no. 3, pp. 1397-1422. DOI: https://doi.org/10.1007/ s10915-017-0501-0.
9. Golushko S. K., Idimeshev S. V., Shapeev V. P. Application of collocations and least residuals method to problems of the isotropic plates theory, Vychisl. Tekhnol., 2013, vol. 18, no. 6, pp. 31-43 (In Russian). EDN: RUDTZF.
10. Belyaev V. A., Bryndin L. S., Shapeev V. P., Golushko S. K. The least squares collocation method with the integral form of collocation equations for bending analysis of isotropic and orthotropic thin plates, AIP Conf. Proc., 2020, vol. 2288, 030084. EDN: FWKXNV. DOI: https://doi.org/10.1063/5.0028298.
11. Belyaev V. A., Bryndin L. S., Golushko S. K., Semisalov B. V., Shapeev V. P. h-, p-, and hp-versions of the least-squares collocation method for solving boundary value problems for biharmonic equation in irregular domains and their applications, Comput. Math. Math. Phys., 2022, vol.62, no. 4, pp. 517-537. EDN: VXTBZN. DOI: https://doi.org/10. 1134/S0965542522040029.
12. Babuska I., Guo B. Q. The h, p and h-p version of the finite element method; basis theory and applications, Adv. Eng. Softw., 1992, vol.15, no. 3-4, pp. 159-174. DOI: https://doi. org/10.1016/0965-9978(92)90097-Y.
13. Ramsak M., Skerget L. A subdomain boundary element method for high-Reynolds laminar flow using stream function-vorticity formulation, Int. J. Numer. Meth. Fluids, 2004, vol.46, no. 8, pp. 815-847. DOI: https://doi.org/10.1002/fld.776.
14. Saad Y. Numerical Methods for Large Eigenvalue Problems. Philadelphia, SIAM, 2011, xvi+276 pp. DOI: https://doi.org/10.1137Z1.9781611970739.
15. Fedorenko R. P. Vvedenie v vychislitel'nuiu fiziku [Introduction to Computational Physics]. Moscow, Moscow Institute of Physics and Technology, 1994, 528 pp. (In Russian). EDN: QJUAEP.
16. Vorozhtsov E. V., Shapeev V. P. A divergence-free method of collocations and least squares for the computation of incompressible fluid flows and its efficient implementation, Vestn. Samar. Gos. Tekhn. Univ., Ser. Fiz.-Mat. Nauki [J. Samara State Tech. Univ., Ser. Phys. Math. Sci.], 2020, vol. 24, no. 3, pp. 542-573 (In Russian). EDN: OSPEBL. DOI: https://doi. org/10.14498/vsgtu1758.
17. Belyaev V. A. On the effective implementation and capabilities of the least-squares collocation method for solving second-order elliptic equations, Vychisl. Meth. Prog. [Num. Meth. Prog.], 2021, vol.22, no. 3, pp. 211-229 (In Russian). EDN: VBQEPE. DOI: https://doi.org/ 10.26089/NumMet.v22r313.
18. Ortega J. M. Introduction to Parallel and Vector Solution of Linear Systems. New York, Springer, 1988, xi+305 pp. DOI: https://doi.org/10.1007/978-1-4899-2112-3.
19. Karamyshev V. B., Kovenya V. M., Sleptsov A. G., Cherny S. G. Variational method of accelerating linear iterations and its applications, Comput. Fluids., 1996, vol. 25, no. 5, pp. 467484. EDN: LDTBIH. DOI: https://doi.org/10.1016/0045-7930(96)00011-4.
20. Il'in V. P., Knysh D. V. Parallel decomposition methods in the trace spaces, Vychisl. Meth. Prog. [Num. Meth. Prog.], 2011, vol. 12, no. 1, pp. 110-119 (In Russian). EDN: OJAYQV.
21. Tang J. M., Nabben R., Vuik C., Erlangga Y. A. Comparison of two-level preconditioners derived from deflation, domain decomposition and multigrid methods, J. Sci. Comput., 2009, vol.39, no. 3, pp. 340-370. DOI: https://doi.org/10.1007/s10915-009-9272-6.
22. Godunov S. K., Ryabenkii V. S Difference Schemes: An Introduction to the Underlying Theory, Studies in Mathematics and its Applications, vol. 19. Amsterdam, North-Holland, 1987, xvii+489 pp.
23. Demmel J. W. Applied Numerical Linear Algebra. Philadelphia, SIAM, 1997, xi+416 pp. DOI: https://doi.org/10.1137/1.9781611971446.
24. Sleptsov A. G. Grid-projection solution of an elliptic problem for an irregular grid, Russ. J. Numer. Anal. Math. Model., 1993, vol.8, no. 6, pp. 519-543. DOI: https://doi.org/ 10.1515/rnam.1993.8.6.519.
25. 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. Math. Phys., 2010, vol.50, no.10, pp. 1670-1681. EDN: OHNDSJ. DOI: https://doi.org/10.1134/ S0965542510100040.
26. Drozdov G. M., Shapeev V. P. CAS application to the construction of high-order difference schemes for solving Poisson equation, In: Computer Algebra in Scientific Computing, Lecture Notes in Computer Science, 8660; eds. V. P. Gerdt, W. Koepf, W. M. Seiler, E. V. Vorozhtsov. Cham, Springer, 2014, pp. 99-110. EDN: UEYEZH. DOI: https://doi.org/ 10.1007/978-3-319-10515-4_8.
27. Golushko S. K., Nemirovskii Yu. V. Priamye i obratnye zadachi mekhaniki uprugikh kom-pozitnykh plastin i obolochek vrashcheniia [Direct and Inverse Problems in the Mechanics of Composite Plates and Shells]. Moscow, Fizmatlit, 2008, 432 pp. (In Russian). EDN: MUWRZD
28. Idimeshev S. V. Modified method of collocations and least residuals and its application in the mechanics of multilayer composite beams and plates, Thesis of Dissertation (Cand. Phys. & Math. Sci.). Novosibirsk, 2016, 179 pp. (In Russian). EDN: VPHODY.
