CC BY
18
УДК 519.612:632.4
H. Д. Золотарёв^, Е. С. Николаев2
СРАВНИТЕЛЬНЫЙ АНАЛИЗ АДАПТИВНЫХ АЛГОРИТМОВ В МЕТОДЕ КОНЕЧНЫХ ЭЛЕМЕНТОВ РЕШЕНИЯ КРАЕВОЙ ЗАДАЧИ ДЛЯ СТАЦИОНАРНОГО УРАВНЕНИЯ РЕАКЦИИ-ДИФФУЗИИ
Предложен новый адаптивный алгоритм построения сетки в /ф-версии метода конечных элементов с кусочно-полиномиальными базисными функциями, ориентированный на нахождение имеющего локальные особенности решения краевой задачи для одномерного уравнения реакции-диффузии. Алгоритм обеспечивает сглаживание сеточного решения путем адаптивного исключения и добавления узлов сетки. Проведено сравнение с предложенным ранее алгоритмом, в котором для адаптивного измельчения сетки и удаления узлов используется оценка локального эффекта от пробного добавления новых и удаления старых базисных функций. Приведены результаты численных экспериментов для оценки эффективности предложенного алгоритма на модельной сингулярно возмущенной задаче с гладким решением.
Ключевые слова: метод конечных элементов, адаптивные методы, индикатор коррекции, стационарное уравнение реакции-диффузии, сингулярно возмущенные задачи.
I. Введение. Для нахождения имеющих локальные особенности решений краевых задач для дифференциальных уравнений наиболее эффективными являются адаптивные варианты метода конечных элементов (МКЭ) или конечных разностей (МКР). Целью адаптации является построение приближенного решения с относительной погрешностью, равномерно распределенной по области интегрирования и не превышающей заданного значения. В контексте /i-версии МКЭ и МКР для достижения этой цели используется последовательное измельчение элементов сетки в тех частях области, где решение имеет особенность, и их укрупнение там, где особенностей нет. В р-версии МКЭ с кусочно-полиномиальными базисными функциями сетка фиксирована и адаптация осуществляется путем увеличения степени полиномов, ассоциированных с ячейками сетки, где решение имеет особенность.
Адаптивные алгоритмы разрабатываются для автоматического выделения указанных частей сетки с малыми вычислительными затратами и основываются на анализе информации о получаемом приближенном решении. Одно из направлений построения таких алгоритмов связано с использованием индикаторов коррекции. В работах [1-8] индикаторы коррекции — это величины, априорно оценивающие изменение приближенного решения, вызываемое пробным расширением или сужением конечноэлементного подпространства, в котором решение ищется, путем добавления новых или удаления старых базисных функций. Алгоритмы, построенные на основе вычисления указанных индикаторов коррекции, будем называть адаптивными алгоритмами пробных базисных функций (ПБФ).
Использовать индикаторы коррекции для построения адаптивной /i-версии МКЭ предложили авторы работы [1]. Этот подход с адаптивным измельчением сеток был реализован в [4] для систем линейных обыкновенных дифференциальных уравнений (ОДУ) 2-го порядка, в [5] — для ОДУ 2-го порядка с измельчением и укрупнением сеток, в [6] — для эллиптического уравнения 2-го порядка в полигональной области, в [7] — для ОДУ 2-го и 4-го порядков при использовании эрмитовых базисных функций. Построению адаптивной иерархической р-версии МКЭ посвящены работы [8, 9]. В [9] исследована особенность р-версии МКЭ, проявляющаяся в том, что приближенные решения могут не меняться на конечном числе последовательно вложенных конечномерных подмножеств увеличивающейся размерности. Найдены необходимые и достаточные условия ее возникновения и предложена модифицированная стратегия, учитывающая эту особенность. В [10] был построен
1 Факультет ВМК МГУ, науч. сотр., к.ф.-м.н., e-mail: nzolQcs.msu.su
2 Факультет ВМК МГУ, вед. науч. сотр., к.ф.-м.н., e-mail: enikQcs.msu.su
адаптивный вариант /ф-версии МКЭ нахождения решения краевой задачи для уравнения реакции-диффузии, исследованию эффективности которого на примере решения сингулярно возмущенных краевых задач посвящена работа [11].
В данной работе для /¡.-версии МКЭ, используемой на первом этапе /ф-версии МКЭ, построен новый адаптивный алгоритм сглаживания решения (СР). В нем, в отличие от основанного на предсказании поведения решения на следующей итерации алгоритма ПБФ, используется информация о решении, полученном на текущей итерации. Он не требует, чтобы приближенное решение было решением минимизационной задачи в конечномерном подпространстве, что делает область его применения более широкой по сравнению с алгоритмом ПБФ.
В предложенном алгоритме индикатор коррекции есть величина, оценивающая величину излома полученного приближенного решения в узлах сетки. Выбор такого индикатора мотивирован тем, что для задач с решением, имеющим непрерывную производную, излома точного решения нет. Кроме того, вычисление такого индикатора коррекции реализуется по простым формулам с малыми затратами. Несложно определить аналог этого индикатора коррекции и для задач, решаемых в многомерных областях.
2. Общая схема /ф-версии МКЭ. Рассмотрим краевую задачу для ОДУ 2-го порядка
Ьи = —(ки')' + ци = /(ж), а < х < Ь, и(а) = /¿о, и(Ь) = ¿¿1, (1)
где к = к(х) и д = ц(х). Предполагается, что выполнены условия, гарантирующие существование единственного решения задачи (1). Используемый в работе для нахождения приближенного решения задачи (1) метод состоит из двух этапов, итерационно реализующих /¡.-версию МКЭ для построения сетки и р-версию для нахождения решения с требуемой точностью. Сначала приведем слабую формулировку задачи (1), а затем дадим описание этапов МКЭ, реализующих указанные версии.
о
Определим множество функций //. соответствующее ему линейное подпространство II. функ-
ционалы J(u) и b(v), а также форму а (и, v) следующим образом:
Н = G W21(a, b), и(а) = ¿¿о, u(b) = ¿¿i j, Я={«е W21(a, b), v(a) = 0, v(b) = 0},
ь ь
J(u) = —а (и, и) — b(u), а (и, v) = j(ku'v' + quv) dx, b(v) = j fv dx.
а а
Пусть форма а (и, и) коэрцитивна, а функционал b (и) ограничен, т. е. имеют место неравенства
о
с||«||2 < а (и, и), ^ С\\и\\ Уи Е Н, с > О, С > О,
о
где ||.|| — норма в Н. Сформулируем задачу минимизации: на отрезке [а,Ь] найти функцию и*(х), такую, что
и* = arg min J(u). (2)
иен
Итерационный процесс нахождения приближенного решения состоит в построении последовательности функций г^1), v,(2\ ..., являющихся решениями задачи минимизации функционала J(u) на конечномерных подмножествах // у, • //д .. • • • множества Н. Соответствующая (2) сеточная задача минимизации на шаге к ^ 1 формулируется следующим образом: на отрезке [а, b] найти функцию такую, что
и(к) _ arg mjn (3)
u£HNk
Известно, что решение этой задачи однозначно находится из условия Галеркина
u{k)EHNk-, а (и{к\ v) = b(v) \/уЕНщ, (4)
О
где HN — соответствующее подмножеству линейное подпространство.
Адаптивная Н-вер сия МКЭ. Приведем формулировку /¡.-версии МКЭ, используемой на первом этапе итерационного процесса. Пусть для к ^ 1
— |а — Xq ^ < ж^ ^ < ... < х^ < —
— построенная на к-й итерации первого этапа сетка, е^ = х^ ^ — интервал этой сетки,
а ср\к\х) — непрерывная кусочно-линейная на [а, Ь] и линейная на интервалах сетки жПк узловая базисная функция, такая, что
-i i
i = О,
<Р
f) (xf)) = 1, 0 < i < п + 1, supp [cp¡k>} = <j e¡k) U e¡% 1 < i < n,
г = n + 1.
e(fc)
Предполагается, что точки разрыва функций к, ц и / являются узлами сетки 7гП)Г
Для к ^ 1 положим Л^ = Пк и определим множество // \ л и связанное с ним подпростран-
о
ство //размерности Л^ следующим образом:
^ «¡¡ + 1 ч
//д, = ^ и = иМк\ ио = МО, ««¡¿+1 = к //д, = е //V,. Мо = 0, /¿1 = 0}.
^ г=0 ^
Заметим, что и(х^) = гц, 0 ^ г ^ пк + 1, для любой функции и € //д,..
о
Известно, что условие Галеркина (4) с учетом определения //и //д л • симметрии формы
а (и, г>) и коэрцитивное™ а (и, и) порождает систему уравнений А^и^ = с симметричной трехдиагональной положительно определенной матрицей, где
А« =
1 ci 1 ai+l
nk
!
i= 1
(fc) -
Uv y = И
П i ' ' ' i u'nk
T
U
«¡¡ + 1 ¿=0
fO) _ ) nQa,[ \ /2 \ ..., fnj-п fhkk + !
(fc) _ № _ UQ ~ Mq! ^¡¡ + 1 ~~ Mi!
(fc) ( (fc) (fc) <4 = -a ( ¥>i_i,
(fc) / (fc) (к)л q = a , Vi I ,
Эта система может быть записана в виде сеточной краевой задачи
1 < i < пк,
(к) (к) (к) (к) _ (к) (к) _ Лк)
иг-1 + сг ui ai+lui+l ~ Ji
(к) _ иа — Мсь
(к) _ Unk + l ~ Mi;
г+1 г+1
решение которой находится методом прогонки. Для генерации элементов системы и вычисления значения функционала на решении задачи (3) можно использовать следующие формулы*:
a¿ =
J (kip'i^ip't + q<pi-i<pi) dx, l^i^n + 1,
Xi- 1
J qv>i dx, f~ = J f<pi dx, 1 < i < n + 1, cn+i = c~+1, /п+1 = /~+1,
ai+i + J q<pi dx, ft = J f<Pi dx, 0 < i < n, c0 = cj, /0 = /0+,
J(u) = -
n+1
Mo(co«o - ai«l - /o)+Ml(cn+l«n+l - ttn+l!ín - fn+l) - y]uifi
i=0
*Для упрощения записи зависимость от к узлов сетки, базисных функций, элементов системы уравнений и приближенного решения не указывается.
Положим* Сг = С~ + с^, /г = + /г+ ДЛЯ 1 ^ '1 ^ П.
В адаптивном варианте /¡.-версии МКЭ переход от множества II , к II для к ^ 2 основан на преобразовании сетки жПк_1 путем добавления на выбранные интервалы новых узлов и удалении некоторых старых узлов, так что
Здесь и 1¿7 — множества номеров измельчаемых интервалов и удаляемых узлов. Если — число элементов множества то Л^ = + для к ^ 2.
При переходе в /¡.-версии от текущей итерации к следующей необходимо перестроить систему базисных функций, сгенерировать новую матрицу и правую часть. В адаптивном варианте метода пересчитывается лишь часть элементов системы уравнений, что позволяет существенно сократить вычислительные затраты, которые приходятся, в основном, на вычисление соответствующих интегралов.
Адаптивная р-версия МКЭ. Пусть выполнены к^ итераций первого этапа метода. Полагая п = пки-, определим:
(Рг(х) = (р\кн\х), хг=х\кн\ 0 < I < П + 1,
аг = а\кк\ 1<г<П+1, * = /г = /?Н\
шп=жп, Ж1=п, А22 = [—а<, Сг, -а»+1]"=1 = Г1'
Приведем формулировку иерархической р-версии МКЭ, используемой на втором этапе для нахождения приближенного решения задачи (2) на сетке
шп = {а = ха < х! < ... <хп< хп+1 = Ь}, ен = (xi-1, ж*).
С узлом х% будем связывать узловую базисную функцию <Рг(х), а с интервалом е^ — последовательность кусочно-полиномиальных неузловых базисных функций 9?г2(ж),9?гз(ж), . . ., таких, ЧТО <Рц(х) есть полином степени ] на еi и
зирр-^у} = Сг, з ^ 2, 1 < г < п + 1.
(к)
На к-й итерации интервалу сц будет соответствовать конечный набор из — 1 неузловых базисных
функций 2 ^ 2 ^ Где ^ ^ Здесь + ^ «¿^ и хотя бы для одного i должно быть
выполнено строгое неравенство.
о о
Подмножество IIи связанное с ним Л^-мерное линейное подпространство //С Н определим следующим образом:
{п+1 п+1 4к) }
и = + По = МО, ип+1=Ц1>,
i=Q г=13=2 )
п+1
Нщ = б Нщ, /х0 = 0, /¿1 = О}, = 4к) ~
г=1
так что IIС //д.. С ... С Н. При к = 1 множество IIпорождают только узловые базисные функции, при этом Ы}. = п. В силу определения базисных функций имеем и(хг) = гц, 0 ^ г ^ п + 1, для любой функции и € Ньтк.
о
Для любого к условие Галеркина (4) с учетом определения //и II и свойств формы а (и, у)
порождает систему уравнений А^и^ = с симметричной положительно определенной матрицей. Опуская индекс к, эту систему при приведенном далее упорядочении неизвестных и уравнений можно записать в следующем блочном виде:
Г1
¡2 ,
11 А12
т А
12 22
и1 и2
*При вычислении с» учтено, что + = 1 для ж»_1 ^ х ^ ж».
где А и = 11А1,..., Ап+1 II — блочно-диагональная, А12 — прямоугольная матрица,
4* =
«1
и„+1
А10 =
VI У 2
V«
Уп+1 г
г =
Ь1 - М0У1
Ь2
ь„
Ьп+1 - М1уп+1
и2 = («1,... ,ип)Т,
У* = (а((рг2, (Рг),.. • , у* = (а (еря, ¡¿>¿-1), • •. ¥>г-1))
г
А,- =
{а(<РИ,
'ргт
и5
г,™=2
= («¿2, • • • Ьг = (^Ь (^¿2) , • • •
Г
а А22 и {2 — определенные выше симметричная трехдиагональная матрица и вектор. Заметим, что в отличие от А22 и {2 матрицы Ац, А12 и вектор Г1 зависят от номера итерации к, причем их размеры увеличиваются с ростом к.
В адаптивном иерархическом варианте р-версии МКЭ переход от множества //д,. , к //д,.
(к)
для к ^ 2 основан на изменении значений л по следующему правилу:
(к) _ (к-1) йг ~~ г
1, ¿еХ/г, О,
к> 2,
,(1)
где — множество выделенных интервалов сетки ц)п. Если 1к — число элементов множества то Щ = + 1к для к ^ 2.
В [10] показано, что вектор и2 есть решение редуцированной системы Аи = ?, где симметричная положительно определенная трехдиагональная матрица размера пх п и правая часть системы определяются следующим образом:
Д — Д ^ДТД_1Д — Г — ^ 22 12 11 12 — I.
" О'/1 с?
"Ог+1 ,=
г=1'
^ _
,Т
А12А111г1 = (Л +М0а1' /2, /п-1, /п+М^п-ы)'
а подвекторы и^ вектора и1 для тех г, при которых ^ ^ 2, вычисляются по формулам
Щ = Щ- 12г - «¿"И^г, 1 < % < п + 1.
Редуцированная система содержит для любого к одинаковое число неизвестных, равное п, и может быть записана в виде сеточной краевой задачи
-ЩЩ-1 + ЩЩ - Щ+1'Щ+1 = /г
1 < % < п,
«О — М(Ь
г^п+1 — Мъ
решение которой находится методом прогонки. Используемые в приведенных выше формулах величины определяются следующим образом:
8. = А, г. = А, V,, = А,- Ч
если > 2,
г г < «'г' г < "г
Ог = а% + Ог, Иг), 1 < % < П + 1, Сг = Сг - (Уг, - (у<+1, г^), /» = /»- (V*, §г) - (Уг+1, §¿+1),
1 < г < п,
где (Уг, 2,г) = (у^, = (у*, 2,^) = (у^, = (у^, §г) = 0, если = 1. Здесь использовано обозначение (а, Ь) для скалярного произведения векторов а и Ь.
При переходе в р-версии от текущей итерации к следующей матрица и векторы Ьг, V^ и уг для 'I € X/; перестраиваются следующим образом: матрица окаймляется строкой и столбцом, а к указанным векторам добавляется по одному элементу. Для % ^ соответствующие матрицы и векторы не меняются. В работе [10] для приведенного метода решения системы (5) была построена эффективная реализация, основанная на выражении матриц А"1, векторов ъ, и V/,. а также матрицы А и правой части ? через соответствующие величины, вычисленные на предыдущей итерации.
3. Адаптивные алгоритмы для //^-версии МКЭ. Принципы построения множеств Xи определяют конкретные алгоритмы адаптации. Заметим, что стандартный вариант /г-версии определяется алгоритмом бисекции всех интервалов сетки, в котором = {1, 2,..., Пк + 1}, а I— пустое множество. Стандартный вариант р-версии определяется алгоритмом, в котором = {1,2,... ,п + 1}.
3.1. Адаптивные алгоритмы ПБФ для //- и р-версий. Алгоритмы пробных базисных функций можно кратко описать следующим образом.
В случае /¡.-версии для каждого интервала е}к ^ сетки жПк_1 вычисляется неотрицательная величина — индикатор коррекции. Эта величина оценивает изменение либо приближенного
решения на отрезке ж^1^ х\к ^ , либо значения функционала возникающее в результате
решения задачи минимизации на множестве, расширенном за счет пробного добавления к множеству базисных функций ассоциированной с этим интервалом и квадратичной на нем неузловой базисной функции ^2.
Аналогично для каждого узла х\к ^ сетки ттПк_1 вычисляется индикатор коррекции
(к-1) (к-1) 1 ■hi+l
либо
Эта величина оценивает изменение либо приближенного решения на отрезке значения функционала Л {и) при пробном удалении этого узла и ассоциированной с ним узловой базисной функции ^(ж) и решении на суженном множестве задачи минимизации.
(к)
В случае р-версии для каждого интервала сц сетки шп вычисляется индикатор коррекции Ц . Эта величина оценивает изменение либо приближенного решения на отрезке [жг_1, ж*], либо значения функционала Л {и) при пробном добавлении к множеству базисных функций ассоциированной
(к)
с этим интервалом неузловой базисной функции где ] = , и решении на расширенном множестве задачи минимизации.
Общая стратегия построения множеств Х^ и определяется следующим образом: номер интервала заносится в соответствующее множество, если значение I^^ и соответственно I¡^ относительно велико. Стратегия построения X¿Г определяется так: номер узла заносится в множество, — (к)
если значение ^ относительно мало.
Эффективность алгоритма ПБФ зависит от затрат на вычисление индикаторов коррекции и стратегии построения множеств Х^, Xи!^. В [10] предложены методы вычисления индикаторов коррекции, не требующие решения задачи минимизации на расширенном или суженном пробном множестве функций. Вместо этого они используют информацию, уже полученную при решении системы = Г^-1) для случая /¡.-версии и редуцированной системы = Г^-1) для случая р-версии.
В [11] была предложена следующая стратегия построения и Пусть — некоторое
множество целых чисел, п — число элементов этого множества, а в — параметр, |0| < 1. Пусть /1, /2,... — некоторый набор положительных чисел. Определим
1
1тгп = ^ 1тах = тах 1и = 1д=(П1Ч ' (6)
геэт ге'ГЯ
/ V Ь " П
п 9
1/п
:%-> и(.к-1)
е = Ф) = ^^гп + (1 - И)'т(„) + (7)
I I (8)
/г
Ь I Ь — 1 ь
где ^ = х\к ^ — х^К Обозначим через в^, вр, £а > 0, £\ > 0, е2 > 0 и е+ > 0 параметры метода, < 1 и \вр\ < 1.
Построение множества X. Сначала вычисляются индикаторы коррекции и величины
(к)
Щ2 ■> которые простым образом выражаются через величины, определяемые при вычислении ин-
дикаторов. В работе [9] было показано, что равенство в!^ = 0 есть необходимое и достаточное условие того, что добавление пробной функции ¡рц не меняет приближенное решение. Следовательно, в этом случае индикатор коррекции будет равен нулю, тогда как приближенное решение 1)
на интервале ег может быть как близко к точному, так и сильно отличаться от него. Поэтому, если ^ £(Ь т0 дополнительно вычисляется г\к\ чтобы различать указанные случаи. Далее
с использованием в (6) I^^ вместо и множества 91+ = |г: ^ £21 вместо вычисляются
1+. , и По формуле (7) вычисляется барьер е^ = и строится множество
ТПЪЛ ТТ№*ЗС ТП
Хк = 1 < г < пк-1 + 1, I^^ ^ или < ^ ео и ^ }■
Построение множества Х^. Для г Е + I ^ Х^ ъ + вычисляют-
ся индикаторы коррекции С использованием в (6) ^ вместо Ц и множества 91" =
= |г: ^ £2 \ вместо 91 вычисляются и Г. По формуле (7) находится барьер
^ } ТПЪЛ ттмхзс с/
£^ = и строится множество Х^ = |г: i Е 91", ^
Если построенное множество X¿7 содержит кластеры (т.е. подряд идущие номера узлов), то исключение на одной итерации всего кластера может породить очень негладкую сетку. Для уменьшения влияния этого эффекта в [11] была предложена процедура, состоящая в прореживании множества Xпо следующему правилу. Кластер с четным числом элементов предварительно преобразуется: из двух его крайних элементов удаляется тот, которому соответствует больший индикатор коррекции Кластер с нечетным числом элементов не меняется. Затем из кластера удаляются
элементы через один, начиная со второго до предпоследнего. Построенное в результате выполнения этой процедуры множество Iи используется при генерации новой сетки.
Построение множества XПостроение этого множества аналогично построению . Вычис-
(к) (к) (к) I (к)I
ляются индикаторы коррекции ц и при 3 = . Если • | ^ е0, то по формуле (8), в кото- 1,(^-1) и (/г —1) (к) п
рои п\ заменено на щ = х^ — и е\ заменено на е*, вычисляется С использованием в (6) вместо Ь и множества 91 = I %: /„■ ^ ^ £2 I вычисляются 1,1 и / . По формуле (7)
\ ; г 1 \ ) тгп' тах д ^ 1 х '
вычисляется барьер £р = е(6р) и строится множество
Хь —
{г: 11к) ^£Р или (о < \n\f\ < £0 и ^ е^ }.
3.2. Адаптивные алгоритмы СР для //-версии. Опишем алгоритм сглаживания решения, предполагая, что решение задачи минимизации на (к — 1)-й итерации первого этапа уже построено. Пусть ломаная*, проходящая через точки Р^ = (х:1, и^), ] = г ± 1,г, есть приближенное решение и(х) на отрезке [жг_1, Жг+1], а Т^ — треугольник, образованный этими точками (см. рис. 1). Имеем
и {Хг - 0) = ---, и {Хг + 0) = ---, /ц = Х^ - Х^-1, /¿¿+1 = - Ж*,
Пг Лг_|_1
= -
1 Хг-1 Щ-1 1 X 2
1 Х1+1 Щ+1
п,
*Для упрощения записи зависимость от к узлов сетки, индикаторов коррекции и приближенного решения не указывается.
так что \Si\ площадь треугольника Т*. Полагая в условии Галеркина (4) v = f>i. получим
iif(xi + 0) — и'(xi — 0) = 77—г / (Lii, — f)(pi dx.
k{Xi) J
■i-i-i
-• • •-
XI Хг-\-\
Рис. 1. Сеточное решение и на [жг_1, ж,+1]
Отсюда следует, что |5г| оценивает снизу норму Ь2 по интервалу (а;*-!, ж^х) от невязки дифференциального уравнения, вычисленной на текущем приближенном решении. В случае относительно большого значения |5г| интервалы е* и вг+1 целесообразно измельчить, а в случае относительно малого значения исключить узел х%.
Индикатор коррекции для удаления узлов определим, полагая = |5г| для 1 ^ ?! ^ п. Заметим,
что индикатор коррекции 1~ равен нулю тогда и только тогда, когда функция и(х) является линейной на отрезке [а;*-!, £¿+1] . В этом случае приближенное решение не имеет излома в точке х^. Индикатор коррекции для измельчения интервалов определим, полагая
/+ = /", /+ = 0.5(7" +/А 2<г<п. /+ =/".
1 1 " г \ г-1 г )■ ' н + 1 п
Для построения множества сначала вычисляется индикатор коррекции ■ Если имеет место одно из условий: ?! = 1 и (-Г-^! ^ во; 2 ^ ?! ^ ^ = пк-1 + 1;
| ^ во, то дополнительно вычисляется г\Ь2 \ Далее с использованием в (6) /г+<,Я'! вместо Ц и множества = { ?!: ^ е-2 \ вместо УХ вычисляются /+ , /+ и По формуле (7) вычисляется
I. ) тл п' тал■ т
барьер е£ = и строится множество
Т~1 = |?!: 1 < ?! < Пк-1 + 1, ^ ££ или ^ £1}.
Построение множества I¿7 выполняется так же, как и в алгоритме ПБФ.
4. Численные эксперименты. Эксперименты проводились в рамках Ир-версии МКЭ. На этапе /г-версии использовался как предложенный алгоритм СР. так и алгоритм адаптации ПБФ, на этапе р-версии использовался алгоритм адаптации ПБФ. Целью экспериментов являлась оценка эффективности предложенного алгоритма на примере модельной сингулярно возмущенной задачи. Было проведено сравнение указанных алгоритмов также и с алгоритмом последовательной бисекции интервалов сетки, используемым в варианте /г-версии МКЭ.
4.1. Тестовая задача. Расчеты проводились для задачи (1), в которой
к(х) = 2е (е:г + е~:г) , д(х) = е 1_:г~, а = -1, Ъ= 1, /х0 = Ц1 = 2, где е > 0 малый параметр. Функция /(ж) задавалась так, чтобы функция
и*(х) = 2 + Мх(1 - |я|)е_с|:г|,
с = е-1'\
являлась решением краевой задачи (1). Нормирующий множитель М был выбран так, чтобы 1 ^ и(х) ^ 3 для х € [—1,1]:
1 г
М =
z( 1 — z)e cz
z =
1 + г + л/l + г2
= 2л/ё.
Рис. 2. Точное: решение и* на [жо, гп+1]
Так как точки х = в которых функция « принимает максимальное и минимальное на [—1,1] значения, принадлежат интервалу (—у/ё, у/ё), то при е< 1 функция и* сильно меняется на малом интервале в окрестности точки х = 0. Функция и* € Н и имеет непрерывную на (—1,1) производную, вид и* приведен на рис. 2 для е = Ю-4. Поведение функции и* в окрестности точки х = О характеризуют следующие данные:
£ 10"3 10"6 10"9
(«*)'( 0) 8.87- 101 2.72 • 103 8.60- 104
z 3.06 • 10"2 9.99 • 10"4 3.16- 10"5
Вторая производная от функции и* (а вместе с ней и функция /) имеет разрыв в точке х = 0. Поэтому эта точка должна быть узлом сетки на любой итерации.
4.2. Базисные функции и квадратурная формула. Используемые в численных экспериментах базисные функции, принимающие на интервале е^ = (ж*-!, х^) не равные тождественно нулю значения, определялись на нем следующим образом:
<Pi-i(x) = ip0(t) Здесь 4'Q(t),4'i(t)
<Pi(x) = 4'l(t), (РгЛх) = ФЛ^: 2 < S < s\ \ |i| < 1, г , hit- Xi = U.5(Xi + Xi — i), hi = 0.5(Жг - Жг-l); 1 < ?' < П + 1. заданные на отрезке [—1,1] мастер-полиномы, такие, что Фа(t) = 0.5(1 - i), 4ч(t) = 0.5(1 + i),
а мастер-полином ips и его производная ф'8 для s ^ 2 определяются через полином Лежандра Ps~i степени s — 1 следующим образом:
t
фs(t) = y"ps_1(c)dc = ^Ty(i2^i)ps,_1(i), Ф'М = Ps-tit), о 2. -1
Значения ips(t) и ip's(t) вычисляются с использованием рекуррентных формул (см. [11]):
g_3 t2 — 1 1Л - I
Фз = t^s-l + —- ^s-г), О 4, ф2^) = ---, Фз(*) = 2 ,
s ^ 2 3i2 — 1
ф'8=ьф'8_1 + —чР'а_2), о 4, = = — -
Реализация метода требует вычисления значений а (и, v) и b(v), где и и v — различные базисные функции. С учетом определения базисных функций эта задача сводится на к-й итерации
(к)
к вычислению для 0 ^ г, s ^ s^ интегралов
ill
krs(i) = J ЬФ'гФ'в dt,, Qrs(i) = J шфгфз dt, fs(i) = J dt, 1 5$ i 5$ П + 1, -1 -1 -1
где ki(t) = k(x), qi(t) = q(x), fi(t) = f(x), а x и t связаны соотношением (9). Для вычисления интегралов применялся метод численного интегрирования с контролем точности — квадратурная формула Гаусса-Кронрода (GK6i) с числом узлов, одинаковым для всех интервалов сетки и каждой итерации [11]. Так как мастер-полиномы не зависят ни от интервала сетки, ни от номера итерации, то целесообразно вычислить значения функций 4j/r4j/s, iprips и ips в узлах квадратурной формулы для требуемых значений индексов г и s и затем использовать эти значения при вычислении соответствующих интегралов.
Значение J(u) на функции и = и* определяется по формуле (см. (15) в [11]):
и
J (и) = —— J \к{и)2 + qu2] dx + кии
— кии
х=Ь
Для вычисления интеграла использовался адаптивный вариант указанного метода отдельно для отрезков [а, 0] и [О, Ц, который строился следующим образом. На первом шаге с помощью формулы Гаусса-Кронрода вычисляется приближенное значение интеграла и оценка погрешности. Если оценка велика, то отрезок интегрирования разбивался пополам и указанная процедура выполнялась на каждой из частей отрезка. Если суммарная погрешность по-прежнему оставалась большой, то делился пополам отрезок с максимальной погрешностью интегрирования. Процесс продолжался до тех пор, пока не достигалась заданная точность или барьер, ограничивающий число делений отрезка.
4.3. Характеристики качества и эффективности метода. В качестве основных показателей, характеризующих эффективность сравниваемых адаптивных алгоритмов, в экспериментах использовались величины:
еи — относительная погрешность приближенного решения на к-й итерации; е'1 — относительная погрешность значения функционала на к-й итерации;
к — величина, оценивающая качество адаптации; ]\[к — размерность подпространства, которому принадлежит приближенное решение; Iк — число интегралов, которое потребовалось вычислить. Приведенные характеристики определяются формулами:
ек = шах ек(г), ек(г) = шах *--——-ек =-——--,
где а к (г) — вспомогательная равномерная сетка на отрезке х\к^], число узлов которой выби-
ралось равным 104. Использование вспомогательных сеток позволяет оценить погрешность в большом числе точек, расположенных между узлами сетки яПк, и тем самым смоделировать оценку для погрешности в норме С [а, Ь].
Величина С,}), оценивает распределение относительной погрешности приближенного решения по области интегрирования и определяется следующим образом:
1 / \ 1/>-к Як = Е1/Е?е(0,1], Е? = Г Т1ек{г), = (Д «*(*)) ,
где 91/; = {г: ^ £}, 1к — число элементов множества 91/;. Выбор (¿к для оценки качества ада-
птации мотивирован тем, что максимально возможное значение (¿к, равное единице, достигается лишь тогда, когда все относительные погрешности вк(г) равны. Другими характеристиками алгоритмов являются: Пк — число внутренних узлов сетки, построенной на к-й итерации; п+ — общее число узлов сетки, добавленных на этапе использования /¡.-версии; п~ — общее число узлов сетки, удаленных на этапе использования /¡.-версии. В приведенных таблицах А^, кр — число итераций, выполненных на этапах применения /¡.-версии, р-версии, а к — общее число итераций. В табл. 1 втах — максимальная степень полиномов в представлении приближенного решения на последней итерации, /¡.Тогп и Ьтах — минимальный и максимальный шаги построенной сетки.
Таблица 1
/г-версия МКЭ /гр-версия МКЭ
Алгоритм бисекции Алгоритм ПБФ Алгоритм СР
£ ю-3 10-® Ю-9 ю-3 ю-® Ю-9 ю-3 ю-® Ю-9
Як ЯкК 2.1 • Ю-3 5.3-Ю-3 1.2 • Ю-5 0.76 6.1 • ю-2 0.73 1.6 • ю-4 0.91 5.1 • Ю-5 0.60 6.7- Ю-2 0.53 1.7- Ю-5 0.83 1.7-10-®
Пк П\ 51199 3199 1638399 51199 3276799 1638399 165 15 9 173 35 49 149 29 89 203 21 9 217 37 49 183 39 89
1к 595198 19353598 29491178 2794 3745 4981 3130 3833 4358
Рь> вкн е\ 5.5 • Ю-7 1.4 • Ю-4 5.1 • Ю-7 5.2 • Ю-4 1.3 • ю-4 5.0- Ю-4 1.5 • Ю-12 1.9 • Ю-1 9.7- Ю-1 1.1 • ю-12 3.7- Ю-1 1.0 • 10° 1.7- Ю-12 7.8-Ю-1 1.0 • 10° 1.4 • Ю-12 1.9- Ю-1 9.7- Ю-1 1.8- Ю-12 3.6-ю-1 1.0 • 10° 1.8 • Ю-12 7.8-Ю-1 1.0 • 10°
ск е-7 е{ 1.9 • Ю-8 4.9 • 10-® 0.6 • Ю-0 5.8-Ю-7 4.5 • Ю-0 1.8 • Ю-8 4.2 • Ю-16 7.3 • Ю-3 4.5 • Ю-2 6.6 • ю-1® 4.7- Ю-4 1.5 • Ю-3 4.3 • Ю-1® 3.5 • Ю-5 4.8 • Ю-5 4.2 • Ю-1® 7.0- Ю-3 4.5- Ю-2 8.7- Ю-1® 4.6 • Ю-4 1.5 • Ю-3 2.2 • Ю-1® 3.5 • Ю-5 4.8 • Ю-5
п+ п~ 48000 1587200 1638400 6 0 10 24 16 76 12 0 28 40 48 98
$тах Ь*тах 1 3.9 • Ю-5 1 1.2 • 10-® 1 6.1 • ю-7 15 2.5 • Ю-2 2.0 • Ю-1 14 1.3 • ю-3 5.2 • Ю-1 14 8.7- Ю-5 8.7- Ю-1 11 2.5- Ю-2 2.0- Ю-1 11 1.2 • Ю-3 6.4 • Ю-1 13 8.7- Ю-5 9.8-Ю-1
fc/l кр 5 0 6 0 2 0 4 24 6 27 9 28 4 27 6 19 9 20
4.4. Результаты численных экспериментов. В табл. 1 приведены результаты решения тестовой задачи для различных значений параметра е. Расчеты проводились при следующих значениях внутренних параметров метода: = 0.25, в^ = ^0.5 и 6Р = 0.25. Во всех расчетах использовались значения: е0 = Ю-15, £1 = Ю-12, £2 = Ю-14, = Ю-7, е = Ю-14. Начальная сетка задавалась равномерной.
В табл. 1 для различных значений параметра е приведены результаты решения модельной задачи с применением: а) стандартного варианта /¡.-версии МКЭ с алгоритмом бисекции всех интервалов сетки; б) варианта /г.р-версии МКЭ с алгоритмом ПБФ; в) варианта /г.р-версии МКЭ, реализующего алгоритм СР. В алгоритме бисекции для очередного е выбор п\ равным полученному значению для предыдущего е позволяет сравнить е" и — относительные погрешности на одной и той же сетке, соответствующие соседним значениям е. Для этого алгоритма в таблице указаны округленные значения ктгП = ктах = /г., где к — шаг равномерной сетки, построенной на последней итерации.
В табл. 2 и 3 для одного значения е представлены распределения узлов конечной сетки и степеней полиномов по заданным частям области интегрирования. В них для каждой части отрезка [—1, 1] приведены: процентное отношение длины этой части к длине всего отрезка (а), п — число попавших на нее узлов сетки и их процентное отношение к общему числу узлов (Ь), ктгП и ктах — минимальное и максимальное расстояние между узлами, а также «Тогп и £тах — минимальная и максимальная степени полиномов, используемых для представления решения на интервалах между узлами.
Табл и ца 2
Алгоритм ПБФ (е = Ю-9)
Интервал а, % п Ь, % ^тгп ^тах $тгп $тах
[-1.0000 -0.0100] 49.50 7 24.14 1.1 • Ю-2 8.7- Ю-1 1 1
[-0.0100 -0.0050] 0.25 1 3.45 5.6-Ю-2 5.6- Ю-3 1 1
[-0.0050 -0.0010] 0.20 2 6.90 1.4 ■ Ю-3 2.8 • Ю-3 1 1
[-0.0050 -0.0003] 0.04 2 6.90 3.5 • Ю-4 6.9 • Ю-4 12 14
[-0.0003 0.0003] 0.03 5 17.24 8.7- Ю-5 1.7- Ю-4 12 14
[ 0.0003 0.0010] 0.04 2 6.90 3.5 • Ю-4 6.9 • Ю-4 12 14
[ 0.0010 0.0050] 0.20 2 6.90 1.4 • Ю-3 2.8 • Ю-3 1 1
[ 0.0050 0.0100] 0.25 1 3.45 5.6-Ю-з 5.6-Ю-з 1 1
[ 0.0100 1.0000] 49.50 7 24.14 1.1 • ю-2 8.7- Ю-1 1 1
Табл и ца 3
Алгоритм СР (е = 10 9)
Интервал а, % п Ь, % ^тгп Ьтах ^тгп Згпах
[-1.0000 -0.0100] 49.50 2 5.13 1.1 • Ю-2 9.7- Ю-1 1 1
[-0.0100 -0.0050] 0.25 2 5.13 1.4 • Ю-3 4.2 • Ю-3 1 1
[-0.0050 -0.0010] 0.20 7 17.95 3.5 • Ю-4 1.4 • Ю-3 1 1
[-0.0010 -0.0003] 0.04 5 12.82 8.7- Ю-5 1.7- Ю-4 3 9
[-0.0003 0.0003] 0.03 7 17.95 8.7- Ю-5 8.7- Ю-5 9 13
[ 0.0003 0.0010] 0.04 5 12.82 8.7- Ю-5 1.7- Ю-4 3 9
[ 0.0010 0.0050] 0.20 7 17.95 3.5 • Ю-4 1.4 • Ю-3 1 1
[ 0.0050 0.0100] 0.25 2 5.13 1.4 • Ю-3 4.2 • Ю-3 1 1
[ 0.0100 1.0000] 49.50 2 5.13 1.1 • ю-2 9.7- Ю-1 1 1
Приведенные в таблицах результаты свидетельствуют, что для модельной задачи алгоритм СР так же, как и алгоритм ПБФ, позволяет получить приближенное решение с высокой точностью на сетках с малым числом узлов в подпространствах, размерность которых слабо зависит от малого параметра е. По сравнению с алгоритмом ПБФ узлы сетки, построенной при использовании алгоритма СР, распределены более плотно в той части области, где решение сильно меняется, и более разреженно там, где решение меняется слабо. При использовании алгоритма СР требуется выполнить меньше итераций второго этапа метода, что приводит к меньшим затратам на вычисление соответствующих интегралов и меньшим временным затратам.
СПИСОК ЛИТЕРАТУРЫ
1. Zienkiewiez О.С., Kelly D.W., deGago J.P., Babuska I. Hierarchical finite element approaches, adaptive refinement and error estimates // The Mathematics of Finite Elements and Applications. London: Academic Press, 1982. P. 313-346.
2. Zienkiewiez O.C., Craig A. Adaptive refinement, error estimates, multigrid solution, and hierarchical finite element method concepts // Accuracy Estimates and Adaptive Refinements in Finite Element Computations. Chichester: John Wiley & Sons, 1986. P. 25-59.
3. Bank R. E., Smith R. K. A posteriory error estimates based on hierarchical bases // SIAM J. Numer. Anal. 1993. 30. N 4. P. 921-935.
4. Николаев E. С. Метод решения краевых задач для систем линейных ОДУ на адаптивно измельчаемых сетках // Вестн. Моск. ун-та. Сер. 15. Вычисл. матем. и киберн. 2000. № 3. С. 11-19. (Nikolaev E.S. Adaptive mesh refinement in boundary value problems for linear ODE systems // Moscow Univ. Comput. Math, and Cybern. 2000. N 3. P. 1-12.)
5. Николаев E. С. Метод решения краевой задачи для ОДУ 2-го порядка на последовательности адаптивно измельчаемых и укрупняемых сеток // Вестн. Моск. ун-та. Сер. 15. Вычисл. матем. и киберн. 2004. № 4. С. 5-16. (Nikolaev Е. S. Method of solving a boundary value problem for a second order ODE on a sequence of adaptively refined and coarsened meshes // Moscow Univ. Comput. Math, and Cybern. 2004. N 4. P. 1-14.)
6. Николаев E.C., Шишкина О.В. Метод решения задачи Дирихле эллиптического уравнения 2-го порядка в полигональной области на последовательности адаптивно измельчаемых сеток // Вестн. Моск. ун-та. Сер. 15. Вычисл. матем. и киберн. 2005. № 4. С. 3-15. (Nikolaev Е. S., S hishki па О. V. A method of solving the Dirichlet problem for the second-order elliptic equation in a polygonal domain on an adaptively refined mesh sequence // Moscow Univ. Comput. Math, and Cybern. 2005. N 4. P. 1-14.)
7. Золотарёва H. Д., Николаев E. С. Метод построения сеток, адаптирующихся к решению краевых задач для обыкновенных дифференциальных уравнений второго и четвертого порядков // Дифференц. уравн. 2009. 45. № 8. С. 1165-1178. (Zolotareva N.D., Nikolaev E.S. Method for constructing meshes adapting to the solution of boundary value problems for ordinary differential equations of the second and fourth orders // Differential Equations. 2009. 45. N 8. P. 1189-1202.)
8. Золотарёва H. Д., Николаев E. С. Адаптивная p-версия метода конечных элементов решения краевых задач для обыкновенных дифференциальных уравнений второго порядка // Дифференц. уравн. 2013. 49. № 7. С. 863-876. (Zolotareva N. D., Nikolaev E.S. Adaptive p-version of the finite element method for solving boundary value problems for ordinary second-order differential equations // Differential Equations. 2013. 49. N 7. P. 835-848.)
9. Золотарёва H. Д., Николаев E. С. О стагнации вр-версии метода конечных элементов // Вестн. Моск. ун-та. Сер. 15. Вычисл. матем. и киберн. 2014. № 3. С. 5-13. (Zolotareva N. D., Nikolaev Е. S. Stagnation in the p-version of the finite element method // Moscow Univ. Comput. Math, and Cybern. 2014. 30. N 3. P. 91-99.)
10. Золотарёва Н.Д., Николаев E. С. Адаптивная hp-версия метода конечных элементов решения краевой задачи для стационарного уравнения реакции-диффузии // ЖВМиМФ. 2015. 55. № 9. С. 15121529. (Zolotareva N.D., Nikolaev E.S. Adaptive Лр-finite element method for solving boundary value problems for the stationary reaction-diffusion equation // Computational Mathematics and Mathematical Physics. 2015. 55. N 9. P. 1484-1500.)
11. Золотарёва H. Д., Николаев E. С. Реализация и исследование эффективности адаптивной hp-версии метода конечных элементов решения краевой задачи для стационарного уравнения реакции-диффузии // Ж. вычисл. матем. и матем. физики. 2016. 56 № 5. С. 69-88. (Zolotareva N. D., Nikolaev Е. S. Implementation and efficiency analysis of an adaptive hp- finite element method for solving boundary value problems for the stationary reaction-diffusion equation // Computational Mathematics and Mathematical Physics. 2016. 56. N 5. P. 764-782.)
Поступила в редакцию 01.02.16
THE COMPARATIVE ANALYSIS OF ADAPTIVE ALGORITHMS IN THE FINITE ELEMENT METHOD FOR SOLVING BOUNDARY VALUE PROBLEM FOR THE STATIONARY REACTION-DIFFUSION EQUATION
Zolotareva N. D., Nikolaev E. S.
The new adaptive algorithm of construction of a grid in the hp-version of a finite element method with piecewise-polynomial basis functions, oriented on solving a boundary value problem having local singularity for the one-dimensional reaction-diffusion equations is offered. The algorithm ensures smoothing of a grid solution by adaptive elimination and adding of nodes of a grid. Comparison with the algorithm offered earlier in which for an adaptive mesh refinement and removal of nodes are based on the estimation of local effect from trial adding new and removals of old basis functions is spent. Results of numerical experiments for an estimation of efficiency of the offered algorithm on the singular perturbed modelling problem with a smooth solution are presented.
Keywords: finite element method, adaptive methods, correction indicator, stationary reaction-diffusion equation, singularly perturbed problems.
