CC BY
45
УДК 519.63; 536.252; 004.75 С.В. СИРИК*
О ВЫБОРЕ СТАБИЛИЗИРУЮЩИХ ПАРАМЕТРОВ В КОНЕЧНОЭЛЕМЕНТНОМ МЕТОДЕ ПЕТРОВА-ГАЛЕРКИНА ПРИ РЕШЕНИИ ЗАДАЧ КОНВЕКЦИИ-ДИФФУЗИИ
Национальный технический университет Украины "КПИ", Киев, Украина_
Анотація. Розглянуто питання вибору та динамічного (адаптивного) налаштування стабілізуючих параметрів вагових функцій у скінченноелементному методі Петрова-Гальоркіна при інтегруванні рівнянь конвективно-дифузійного типу. На ряді обчислювальних прикладів продемонстровано ефективність відповідного методу Петрова-Гальоркіна.
Ключові слова: метод Петрова-Гальоркіна, метод скінченних елементів, задачі конвекції-дифузії, вагові функції, стабілізуючі параметри, стабілізовані методи.
Аннотация. Рассмотрены вопросы выбора и динамической (адаптивной) настройки стабилизирующих параметров весовых функций в конечноэлементном методе Петрова-Галеркина при интегрировании уравнений конвективно-диффузионного типа. На ряде вычислительных примеров продемонстрирована эффективность соответствующего метода Петрова-Галеркина.
Ключевые слова: метод Петрова-Галеркина, метод конечных элементов, задачи конвекции-диффузии, весовые функции, стабилизирующие параметры, стабилизированные методы.
Abstract. The choice and dynamical tuning of stabilization parameters of the weighting functions in finite-element Petrov-Galerkin method for solving convection-diffusion problems are considered. The efficiency of the corresponding Petrov-Galerkin method is illustrated and confirmed by a number of computational examples.
Keywords: Petrov-Galerkin method, finite element method, convection-diffusion problems, weighting (test) functions, stabilization parameters, stabilized methods.
1. Введение
В настоящее время метод Петрова-Галеркина (МИГ) в форме метода конечных элементов (МКЭ) является одним из наиболее универсальных методов построения численных схем для решения задач математической физики [1]. При использовании МИГ для расчета задач с доминированием конвективных процессов ключевым вопросом является корректный выбор весовых функций (ВФ), который предотвращал бы появление в численном решении нефизических осцилляций и неустойчивостей (обеспечил стабилизацию [1] численного решения) при сохранении приемлемой его точности. В работе [2] предложены структуры кусочно-полиномиальных ВФ МПГ для решения двухмерных задач конвекции-диффузии (КД), а в работе [3] - для интегрирования задач с тремя пространственными переменными. Для настройки формы функций используются стабилизирующие параметры (СП), связанные с ребрами сетки разбиения. Анализ аппроксимации, устойчивости и сходимости для данных ВФ [2, 3] (а также их кусочно-степенных обобщений) для нестационарных одномерных случаев проведен в работе [4]. Данные ВФ впоследствии применялись для численного решения различных нестационарных задач КД (в том числе и в тех случаях, когда скорость в конвективном слагаемом с течением времени резко изменяется как по величине, так и по направлению), а также нелинейных уравнений Бюргерса [5, 6] и уравнений магнитной гидродинамики [7]; в цитированных работах также были предложены методы (способы) выбора СП, реализующие идею адаптивности (возможности подстраивать их под эволюционирующее во времени решение для обеспечения устойчивости счета).
В настоящей статье результаты работ [5-7], касающиеся выбора СП, получили дальнейшее развитие, и на примерах (считающихся "тяжелыми" для численного решения
© Сирик С.В., 2014
ISSN 1028-9763. Математичні машини і системи, 2014, № 4
139
по причине наличия в них доминирующих конвективных процессов, пограничных или внутренних слоев, решений солитонного типа и прочих особенностей, доставляющих осложнения при численном счете) показано эффективность ("робастность" [1]) предложенных вариантов выбора СП в плане точности и устойчивости численного решения.
2. Постановка задачи
Будем считать, что в качестве базисных функций {N\ (х)} (каждая функция Nt (х) ассоциирована с соответствующим узлом х;) используются кусочно-линейные финитные функции [1-7]. В качестве ВФ Wi (х) используем функции работ [2-4]. Обозначим через W(i) многоугольник, образованный объединением тех элементов , для которых узел i является одной из вершин. Множество W(i) является носителем для функции Ni, а также для функции Wi, которая определяется как [2-4]:
W (х) = Ni (х) + £ aLtW(it)(x), (1)
kkKt
где Ki - множество номеров узлов k, каждый из которых соединен одним ребром с узлом
i, ai k - СП, соответствующий ребру (i, k) и позволяющий регулировать изгиб ВФ на нем.
Подробная конструкция функций W(i,k), а также всей ВФ Wi, приведена в [2, 3]. Целью
данной работы является выбор совокупности СП {ai,k} и их динамическая настройка, при
котором бы обеспечивалась стабилизация (устранение нефизических осцилляций) численного решения при сохранении приемлемой точности получаемого решения.
3. Динамическая настройка стабилизирующих параметров
При интегрировании полудискретных аппроксимаций [1-7] в виде системы обыкновенных дифференциальных уравнений (СОДУ), построенной МПГ для аппроксимации уравнений КД, на n + 1-м шаге интегрирования предлагается выбирать СП {ai k } ВФ Wt (х) в виде
a
(n +1) =
i,k
(2)
где b(n+1) - некоторое фиксированное неотрицательное число (зависящее от номера шага), [bz, | bz | < 1
функция Fn (z) ° функция a(g) ° coth(g) -1/ g (при g = 0 полагаем, что
н [sign z, | bz | > 1
a(0) = lim(coth(g)-1/g) = 0), число Т® =(u(n)•Kk)/(2k(n+1)), h,k °хк -хг, u{n) - вектор
g®0
скорости (векторная величина, которая умножается скалярно на градиент неизвестной искомой величины в рассматриваемом уравнении типа КД) в i -м узле на предыдущем n -м
шаге интегрирования СОДУ (в качестве u (n) также можно взять взвешенное среднее данных векторов в узлах i и k [7]), k(n+1) - коэффициент, регулирующий величину вводимой
в расчетную схему искусственной вязкости относительно i -го узла на n + 1 -м шаге интегрирования (соответственно, в силу определения a(g), чем меньше значение k(n+J), тем большей будет искусственная вязкость). При изотропном (диагональном) диффузионном тензоре, в большинстве случаев, как показывают расчеты для устранения нефизических осцилляций, выбор данного коэффициента вязкости в виде k (где k - диффузионный ко-
140
ISSN 1028-9763. Математичні машини і системи, 2014, № 4
эффициент при диффузионных членах решаемого уравнения КД) оказывается вполне достаточным. Параметр 0(n+1) - некоторая константа, увеличение или уменьшение которой
приводит соответственно к увеличению или уменьшению вводимой в численную схему искусственной диссипации. В большинстве случаев для нужд практических расчетов вполне достаточно выбирать ее из промежутка [0; 2]. Увеличение этой константы может понадобиться при возникновении в решении больших градиентов и нефизических осцилляций в соответствующих тонких слоях (поскольку увеличение искусственной диссипации приводит к "монотонизации" численных схем и сглаживанию неоднородностей, разрывов
и скачков). Заметим, что при = 0 по определению получаем классический МКЭ Га-
леркина. Положив в (2) b(n+1) = +¥ и b(n+1) = 1, получим соответственно
a
(n+1) i ,k
= 0
(n+1)
i,k
sign (a(g (nj)) ) = q
(n+1) i ,k
sign (її),
(3)
a
(n+1) = i,k
qS^aefi1).
(4)
Отметим, что выражения (2)-(4) для СП пригодны не только при интегрировании полудискретных аппроксимаций, но без изменений могут быть использованы при нахождении решений разностных схем (при переходе от предыдущего n -го временного слоя схемы на текущий n + 1-й). Естественно, при интегрировании систем уравнений настройка СП в ВФ для каждого уравнения системы осуществляется индивидуально. Для определения СП при решении стационарных задач используются те же формулы (2)-(4) (однако, теперь они априори являются независимыми от n ). При 0(flj+^'> = 1 формула (4) дает выражение, которое
для стационарного линейного одномерного уравнения КД с постоянными коэффициентами обеспечивает совпадение численного и точного решений в узлах равномерной сетки (то есть, получаем схему Ильина-Аллена-Саусвелла [1]). Дополнительные замечания, касающиеся выбора СП, упрощения их вычисления, а также стабилизации решений в окрестностях фронтов ударных волн приведены в [5] (см. замечания 1-5 в разделе 2 работы [5]).
4. Численное моделирование
Пример 1. Рассмотрим граничную задачу для стационарного двухмерного уравнения КД v • Vu = кДи (тут v - вектор скорости, к - диффузионный коэффициент) в единичном квадрате W (центр квадрата находится в начале координат). Для упрощения записи будем
обозначать х ° хь у ° Х2 . Поле скоростей имеет вид v(x, y) = (—y; x) (вращение против
—8
часовой стрелки), к = 10 . На границе квадрата W решение и равно нулю, а внутри квад-
рата, на прямой х = 0, —1/2 < у < 0 имеется разрез, на котором и = и(0, у) = w(y), где w( у) = (cos(4py + p) +1)/2. Постановка граничных условий и линия разреза изображены на рис. 1 [8]. На рис. 2 показан график функции w( y) при —1/2 < y < 0. На рис. 3 показано решение данной задачи [8] - "воронка", которая образуется при вращении графика функции w(y) против часовой стрелки относительно начала координат. Данная задача использовалась в [8] для тестирования точности и работоспособности стабилизированных методов SUPG, GLS, USFEM и RFB (однако, в [8] было взято к = 10—6). Поскольку процесс конвекции в задаче является доминирующим над процессом диффузии ( к - малое), то при вращении графика w(y) его высота должна оставаться практически такой же, какой она была на разрезе (то есть график не "проседает" под действием диффузии, распределенной в данном случае изотропно) [8]. Это, в свою очередь, для оценки точности и качества чис-
ISSN 1028-9763. Математичні машини і системи, 2014, № 4
141
ленного метода позволяет использовать максимальную разницу между высотой w(y) на линии разреза и высотой графика вычисленного решения на вертикальной линии узлов, находящихся слева от линии разреза на ближайшем расстоянии (для расчетов будем использовать сетку с прямоугольными треугольниками, ориентированными вдоль координатных осей). Обозначим данную разницу в текущей задаче через errmax . Отметим, что в
силу особенностей данной задачи ее можно рассматривать как "стресс-тест" для проверки того, вводит ли численный метод избыточную диффузию поперек потока [1, 6, 8]. В случае утвердительного ответа при полном повороте вокруг начала координат график под воздействием диффузии будет "разглажен", а величина errmax, соответственно, будет большой [8]. Таким образом, данная задача предъявляет высокие требования к численным методам (в особенности, к их диссипативности) и считается весьма тяжелой для численного решения [8].
и=0
11=0 У 1/2
s' ~~
\
у
и=Л\’(у)
и=0
и=0
1/2
Рис. 1. Постановка граничных условий (пример 1)
Рис. 2. График w = w(y) .
Рис. 3. Пространственный профиль решения (пример 1)
Для метода SUPG используем поэлементную стабилизацию [1, 3] с СП {at} вида
ax= hxa(Pex)/ (21| v ||^ (5)
где ht - диаметр элемента Tt, Pet = || v || ht /(2к), || v || - евклидова норма вектора v , вычисленного в центре элемента Tt [8] по поводу данного выбора параметров. Данный выбор параметров в SUPG является оптимальным для одномерных стационарных задач КД с постоянными коэффициентами (дает решение, совпадающее с точным решением в узлах) и считается также вполне подходящим для большинства задач (в том числе, многомерных) [1, 8, 9] (особенно таких, в которых процесс конвекции не сильно доминирует над остальными процессами). Используем сначала равномерную сетку с треугольными элементами и числом узлов 21х 21. Для МПГ с ВФ (1) (используем здесь, для определенности, кусочноквадратичные ВФ) и при выборе параметров в виде (4) для всех Qik = 1/10 получаем
errmax » 0,097, а график решения практически совпадает с истинным решением (рис. 3); при возрастании или убывании Qik от данного значения погрешность растет. Так, для %ік = 0 (МКЭ Галеркина) получаем errmax » 0,145 (а в численном решении имеются осцилляции, что противоречит физическому смыслу задачи), при Qik = 1/2 будет errmax » 0,109 , а при Qik = 1 получим errmax » 0,119. Метод SUPG с выбором {at} в виде (5) дает errmax » 0,129. В работе [8] для решения данной задачи также применялся метод
SUPG с выбором {at} в виде at
ht cos f + sin f
21| v || (1 + 3cos f sin f)
a(Pet), где f - угол (локально-
142
ISSN 1028-9763. Математичні машини і системи, 2014, № 4
го) направления потока (относительно декартовой системы координат, см. подробно в [8]). В работе [8] показано, что для стационарных двухмерных задач такой выбор стабилизирующих параметров в методе SUPG обеспечивает наилучшую точность в случае доминирования конвекции (по сравнению не только с другими вариантами выбора параметров, но, как показали расчеты [8], также и с другими стабилизированными методами, в особенности, методами GLS, USFEM и RFB); метод SUPG с таким выбором {at} дает
errmax » 0,117.
Рассмотрим теперь результаты расчетов для равномерной сетки 31х 31 узлов. Снова же, при Qik = 1/10 получаем минимальную погрешность errmax » 0,043. При Qik = 0
получим errmax » 0,072 , при qi k = 1/2 — errmax » 0,046 , при qi k = 1 - errmax » 0,048. Для метода SUPG с (5) получим errmax » 0,05, а для вышеописанного оптимального выбора параметров {at} - errmax » 0,046 . Для МПГ с ВФ (1), параметрами (4) и 0i k = 1/10, а также
для метода SUPG c вышеописанным оптимальным выбором параметров при использовании разбиения 21х 21 узлов число обусловленности (в спектральной норме [10]) »126 (для разбиения 31х 31 узлов приближенно равно 206); для метода SUPG c параметрами (5) и использовании разбиения 21х 21 число обусловленности »210 (для разбиения 31х31 узлов - приближенно равно 352). Таким образом, расчеты свидетельствуют, что предлагаемая в данной работе версия МПГ способна в данной задаче обеспечить большую точность решения по сравнению с другими стабилизированными методами и, кроме того, по сравнению с методом SUPG, формирует систему уравнений с меньшим числом обусловленности и не вводит избыточную диффузию поперек потока (что следует из определения метода SUPG и рассчитанных значений errmax ).
Пример 2. Рассмотрим задачу для нестационарного двухмерного уравнения КД:
Эu du du Э f Эu) Э f Эu ^
+ Vi + v2 = Эt Эxl Эx2 Эxl к— 1 ^1J + ЭX2 к 1 ^2 7
(6)
в прямоугольной области, 0 < xt < Ц (i = 1, 2) при постоянном коэффициенте к и векторе
скорости V = (v1; v2)T = const. Начальное условие задается в виде u(0,x1,x2) = eX1 . Граничные условия являются продолжением по непрерывности начальных условий на границу области: u(t,0,x2) = 1, u(t,L1,x2) = eL1, u(t,x1,0) = u(t,x1,L2) = eX1. Аналитическое решение данной задачи для (6) выражается в виде [2]
u (t, xb x2 ) = e
¥ ¥ TT
at+bx1 + cx2^PST' Hij
A ^“1 Wii - a
і=1 j=1 ij
—at ~Wiit
e~"j I sin pix1sin pjx2 + ex1 L1 L2
где
rj л 2• .(n i\2 2 )1 + eL1(1-b)(-1)i+1 1 + e-L2C (-1)j+1
Hj = 4p j(k(l-b) +KC +a) r2„ „2 2.2 • ,2 2. 2 2
Т2ҐЛ ,2 , 2-2 т-2 2 , 2-2
Lj(1 -b) + Pi L2c +P j
Wjj =K((Pi / L1)2 + (P / L2)2).
Определим следующие значения параметров задачи: L = L2 = 1, V1 = V2 = 50,
к = 10-1. Поскольку || V || / к » 707, то в данной задаче конвективные процессы значительно преобладают над процессом диффузии.
ISSN 1028-9763. Математичні машини і системи, 2014, № 4
143
Для интегрирования возникающих СОДУ во всех приведенных примерах нестационарных задач использовались стандартный явный метод Рунге-Кутты 4-го порядка [11] и явный адаптивный метод 3-го порядка из работы [12]. Шаг по времени в текущем примере
_3
положим t = 10 (начальный для метода из [12]. При этом в настройках выставлялись
значения допустимой абсолютной и относительной ошибок [12] интегрирования, также
_3
равные 10 ). Отметим, что в приведенных примерах оба метода интегрирования приводят
к практически идентичным результатам. Отметим также, что рассмотренные и проанализированные в работе (стабилизированные) методы конечных элементов, а также предлагаемый вариант МПГ, осуществляют дискретизацию задач по пространственным переменным, в результате чего получается СОДУ. Для интегрирования получившихся СОДУ в большинстве случаев оказывается вполне достаточным использования методов низкого порядка [1, 8, 9, 13]. Однако, поскольку погрешности и неустойчивости при применении конечноэлементных методов к решению задач КД обусловлены преимущественно дискретизациями пространственных членов и, следовательно, имеют скорее "пространственное" происхождение, чем "временное" [1-9, 13], то в данном разделе при рассмотрении примеров нестационарных задач для интегрирования СОДУ были использованы методы высоких порядков (3-го и 4-го). Это было сделано для того, чтобы минимизировать влияние дискретизации по временной переменной на окончательный результат и в явном виде выделить преимущество МПГ с ВФ (1) над другими стабилизированными методами пространственной дискретизации, рассмотренными в примерах.
Рассмотрим сначала результаты расчетов на равномерной сетке 15 х15 узлов с треугольными элементами. Гипотенузы треугольников ориентированы поперек потока (известно [8], что уменьшение угла между направлением потока и гипотенузами/диагоналями элементов приводит, как правило, к уменьшению погрешности, поэтому мы будем тестировать методы на наиболее "неудобном" с данной точки зрения расположении сетки). График численного решения МКЭ Галеркина в зависимости от переменной x ° x при t = 0,4 и фиксированном Х2 = 0,5 приведен на рис. 4 (здесь и в дальнейшем на рисунках данного примера тонкой линией будет обозначаться истинное решение, а жирной - вычисленное). Видно, что численное решение является колебательным и не имеет ничего общего с истинным решением. Обозначим через err^lX величину максимального уклонения на графи-
ке истинного решения от вычисленного, а через err
glob
величину максимального (пото-
чечного) уклонения на всей области. Тогда для МКЭ Галеркина имеем err^lX »1,47,
а
err^0 » 3,34. Для метода SUPG с выбором параметров (5), а также для методов GLS и
USFEM получаем err°x » 0,72, а errmg°b »1,55 (рис. 5). При этом для интегрирования
СОДУ адаптивным методом [12] понадобился M = 341 шаг, а показатель 1 жесткости системы (отношение максимального модуля действительных частей собственных чисел к минимальному) равен 3,76. Видно, что данные методы сильно сглаживают решение (стремятся его "выровнять" в окрестности пограничного слоя), что свидетельствует об их чрезмерной диссипативности в данной (нестационарной) задаче. Отметим, что использование метода SUPG с оптимальным для двухмерных стационарных задач выбором параметров [8] (см. пример 1) не приводит к существенному улучшению решения. Для него
err^lx » 0,64 , err^lf » 1,54, 1» 1,88, M = 232, однако решение имеет осцилляции в окрестности пограничного слоя (рис. 6).
144
ISSN 1028-9763. Математичні машини і системи, 2014, № 4
Рис. 4. Решение МКЭ Рис. 5. Решение SUPG Рис. 6. Решение SUPG c
Галеркина (пример 2) c выбором СП в виде (5) оптимальным выбором СП
и Qik = 1,5 (15ХІ5 узлов) и Qik = 1,75 (15ХІ5 узлов) и Qik = 1,75 (25X25 узлов)
Для МПГ с ВФ (1) (для определенности используем кусочно-квадратичные функции), параметрами (4) и 0i к = 1 получим errox » 0,61, err» 1,53, 1» 1,64, M = 233 .
Однако решение будет иметь небольшие осцилляции возле приграничного слоя (график решения практически совпадает с графиком на рис. 6, поэтому здесь не приводится). При
0i к = 1,5 получим errlOC » 0,66 (рис. 7), M = 252 и 1» 2,45 . Видно, что в окрестности пограничного слоя имеется тенденция к осцилляции (образуется небольшая "ступенька"), хотя и значительно меньшая, чем для метода SUPG (рис. 6). Увеличение 0i к позволяет удалить эту "ступеньку". Так, график решения для 0i к = 1,75 приведен на рис. 8 (при этом
err°x » 0,67 , errm0 » 1,52, 1» 2,98 и M = 290). Отметим, что использование параметров (3) в данном примере дает слегка худшие результаты. Так, для 0i к = 1,75 будет
err^ » 0,66 , M = 363 , однако в окрестности пограничного слоя наблюдается "ступенька" (аналогично ситуации на рис. 7), убрать которую можно еще увеличив 0i к (однако при
этом будет расти err^X ). Измельчение сетки ведет к уменьшению погрешности (хотя соотношения между погрешностями (в плане "больше" или "меньше") протестированных методов с различными вариантами СП, а также характер поведения решения в целом сохраняются). Так, для разбиения 25 X 25 узлов для МКЭ Галеркина получаем err^X » 0,64, err^X»3,47, 1» 2,02 и M = 551. Для метода SUPG с (5) будет err^X »0,5,
err£ » 1,49, 1» 3,69, M = 678 . Для рассматриваемого МПГ с ВФ (1), параметрами (4) и 0i к = 1 получим errjOC » 0,34, єгг^Х » 1,47 , 1» 1,65 , M = 502; при 0i к = 1,5 будет
ISSN 1028-9763. Математичні машини і системи, 2014, № 4
145
errmx » 0,42, err^ »1,48, 1» 2,36 и M = 583; при 0a = 1,75 будет err^ » 0,45, errmax » 1,48 , 1» 2,96 и M = 630 (рис. 9).
Пример 3. Рассмотрим начально-краевую задачу при х є [а; b] для одномерного уравнения Бюргерса [5, 6]:
Ъы
dt
л Ъы
+ 1ы— Ъх
_э_
Ъх
к-
Ъы
Ъх
+ f
(7)
при u(t, х)
1 = 1
и
f = 0
х /(t +1)
1 + д/(t +1) / exp(1 / 8к) exp(х2 /(4k(t +1)))
известным аналитическим решением Данная задача представляет собой мо-
с
дель ударной волны [14] и применялась в [14-17] для тестирования предложенных там численных методов, основанных на методе RKPM и являющихся его дальнейшим усовершенствованием [14], а также методах, основанных на применении В-сплайнов, эрмитовых сплайнов и методов коллокаций [14-17]. Начальное (при t = 0) и граничные условия получаются из выписанного аналитического решения путем его непрерывного продолжения на гиперплоскости с t = 0 и х = а , х = b соответственно. Для оценки уклонения численного решения и от аналитического и здесь и в дальнейших примерах используем величину
errmax °maxy |~~(х(,t)-и(х(,t)|. Положим, к = 10-4, а = 0, b = 1. При использовании (3)
положим [5, 6] 0(П+1) = 0 • m, где m = max | a(g'f‘1)) |, константа 0 = 1, а при использова-
(«)-
i,k
нии (4) положим 0(”k+1^ = 1. В табл. 1 приведено значение величины errmax (при t = 2) в зависимости от выбора СП a(?k+1^ и количества N равноотстоящих узлов сетки.
Таблица 1. Значение errmax (пример 3)
„(«+!) ai,k Количество узлов пространственной сетки, N
200 250 300 350
(3) 0,33203 0,32471 0,31836 0,31259
(3)* 0,12784 0,11491 0,08867 0,06558
(4) 0,19181 0,18664 0,18053 0,15949
(4)* 0,28647 0,28646 0,28644 0,28643
* при расчете использовалось сосредоточение !5, 18].
Здесь и в дальнейших примерах на рисунках жирной линией обозначен график вычисленного решения в зависимости от х , тонкой - график известного аналитического. На рис. 10 представлен результат расчета задачи МПГ (здесь и в дальнейших примерах в МГП используем кусочно-квадратичные ВФ (1)) при использовании формулы (3) и сосредоточения [5, 18], N = 350. На рис. 11 область укручения решения показана крупным планом. На рис. 12 показана область укручения решения при использовании формулы (4), N = 350, а на рис. 13 - при использовании МКЭ Галеркина с сосредоточением. Видно, что в последнем случае решение является осциллирующим (в силу преобладания конвекционных процессов и нелинейности); отметим, что в случае, когда сосредоточение не используется, осцилляции будут гораздо большими, поскольку сосредоточение, как известно [5, 18], увеличивает устойчивость схем и осуществляет их "монотонизацию". При использовании (3) без сосредоточения, а также при использовании (4) с сосредоточением погрешность боль-
146
ISSN 1028-9763. Математичні машини і системи, 2014, № 4
ше (по сравнению с другими вариантами, перечисленными в табл. 1) и с возрастанием N она медленно убывает. При использовании МПГ с постоянными СП (не зависящими от индексов узлов сетки и момента времени t) решение получается неустойчивым и колебательным в окрестностях тонких слоев, а в случае, когда применяется сосредоточение — чрезмерно диссипативным и сглаженным. Так, на рис. 14 изображен результат решения при прежних условиях, с использованием сосредоточения и СП, тождественно равными единице. Видно, что решение чрезмерно сглажено в окрестности слоя. На рис. 15 показан соответствующий результат, когда сосредоточение не используется. Как видим, в таком случае появляются осцилляции. Увеличение СП не приводит к улучшению решения и исчезновению осцилляций (а при использовании сосредоточения решение становится еще более сглаженным, чем на рис. 14). Уменьшение СП также не приводит к улучшению решения (а при стремлении их к нулю получим МКЭ Г алеркина, который порождает осцилляции независимо от того, использовалось сосредоточение или нет подобно ситуации на рис. 13). При интегрировании СОДУ здесь и в дальнейших примерах использовались методы, описанные в примере 2 (в данной задаче с шагом по времени t = 10- ).
Рис. 10. Решение МПГ с СП (3) и
использованием сосредоточения
Рис. 12. Решение МПГ с СП (4) (область укручения)
Рис. 11. Решение МПГ с СП (3) и сосредоточением (область укручения)
Рис. 13.Решение МКЭ Галеркина с использованием сосредоточения
Рассмотрим теперь данную задачу при условиях к = 0,005 , a = 0, b = 1,2 с шагом по времени t = 0,01 и по пространству h = 0,0012 (то есть N = 1001). При таких значениях параметров задача рассматривалась в [14] (однако там был взят еще более мелкий шаг h, h = 0,001). Тогда при использовании (3) ( 0 = 1) с сосредоточением для моментов времени
t = 1,7, t = 2,5, t = 3 , t = 3,5 имеем err„
67 10
-5
err
62 • 10
-5
err
89 • 10
-5
err
36 • 10 5 соответственно. Для метода работы [14] с использованием пространств
кусочно-линейных функций соответствующие значения errmax равны errmax » 13,47 • 10
-4
errmax » 15,55 -10
-4
err
15,53 • 10-4, errn
15,22 • 10-4 [14].
ISSN 1028-9763. Математичні машини і системи, 2014, № 4
147
Рассмотрим еще задачу при условиях к = 0,005, а = 0, b = 1 с шагом по времени t = 0,01 и по пространству h = 0,005 . При таких значениях параметров задача рассматривалась в [14, 15]. Тогда, аналогично, при использовании (3) с сосредоточением для моментов времени t = 1,7, t = 2,5, t = 3,25 имеем errmax » 48 • 10-5, errmax » 18,5 -10-5 , errmax » 63 • 10-5 соответственно. Для метода QRKM (Quadratic Reproducing Kernel function
_3
Method), предложенного в работе [14], соответственно имеем errmax »0,092 10 ,
errmax » 0,115 10_3, errmax » 8 10_3 (табл. 11 на стр. 433 в [14]). Для методов QBCM (Quadratic B-spline Collocation Method) и CBCM (Cubic B-spline Collocation Method), предложенных в работе [15] и основанных на использовании В-сплайнов, имеем соответственно
значения [14, 15] errmax » 0,312 10_3, errmax » 0,189 10_3, errmax » 8,98 10_3 и
errmax » 27,577 • 10_3, errmax » 25,152 10_3, errmax » 21,049 10_3 . Таким образом, приведенные расчеты показывают, что предлагаемая версия МПГ с адаптивным выбором СП в виде (3) и использованием сосредоточения дает лучшие результаты (в том числе даже и на более грубых сетках), чем методы, предложенные в работах [14, 15].
Рис. 14. Решение МПГ с СП, равными единице, и сосредоточением
равными единице
Для сравнения c методами, предложенными в работах [16, 17], рассмотрим данную задачу при условиях к = 0,005 , а = 0, b = 1, t = 0,001 (при таких условиях она рассматривалась в указанных работах). Тогда при использовании (3) (0 = 1/2) с сосредоточением при N = 61, N = 81, N = 101 и N = 121 для фиксированного момента t = 3,6 имеем
» 0,134 10_3, errm^ » 0,088• 10_3, errm^ » 0,059 10“"
errmax » 0,238 10
err
соответ-
ственно. При использовании (4) (при прежних условиях) соответственно имеем
errmax » 0,354 -10
_3
errmax » 0,146 10
_3
err
0,156 10
_3
err
0,046 10 3. Для ме-
-3
тода работы [16] будем иметь errmax » 0,47 10 , err,
0,27 • 10
-3
err
0,17 • 10
-3
err
-3
0,12 10 соответственно (см. [16]), а для метода работы [17] — соответственно
0,45 40
-3
err
0,31 10
-3
err
: 0,23 40
-3
err
0,16 10_3 [16,
получаем errmax
17]. Таким образом, данные расчеты показывают, что предлагаемая версия МПГ с адаптивным выбором СП в виде (3) и использованием сосредоточения, а также с выбором параметров в виде (4) (без использования сосредоточения) дает в данном примере более точные результаты, чем методы работ [16, 17]. Отметим, что для интегрирования полудискретных аппроксимаций в [16] использовались методы Рунге-Кутты SSP-RK43 и SSP-RK54 (3-го, 4-го и 5-го порядков) (детально в [16]).
148
ISSN 1028-9763. Математичні машини і системи, 2014, № 4
Пример 4. Рассмотрим начально-краевую задачу при х є [a; b] для одномерного уравнения Бюргерса (7) при 1 = 1 и f = 0 с известным аналитическим решением
u(t, х)
2 лк - exp(-p2 kt) • sin(px) 2 + exp(-p Kt) • cos(px)
. Данная задача применялась в [16, 19] для проверки и тес-
тирования предложенных там численных методов решения уравнения Бюргерса. Положим,
a = 0, b = 1, h = 1/40, t = 10-4 . Тогда при использовании (3) с сосредоточением (подроб-
_3
ный выбор СП описан в примере 3) в фиксированный момент времени t = 10 для
Ї-5 ___ о о ю і і а—6
к = 5/10, к = 2/10 и к = 1/10 имеем errmax »1,2741 10_
■ max ■
errmax » 2,2121-10
и
err„
0,5987 10 6 соответственно. При использовании формул (4) соответственно полу-
чим errmax » 4,0752 -10
max
-5
errmax » 6,6653 10
-6
err
1,6787 10 6. Соответствующие
значения величины errmax для
errmax »1,22 -10-5 и errmax » 3,08 10
errmax » 3,00 -10 -6 и errmax » = 2,00 • 10
для метода работы [16] равны errmax »7,44 10
-5
-6
2,00 • 10 6 (табл. 3 работы [19] на стр. 575). Таким образом,
приведенные расчеты показывают, что предлагаемая версия МПГ с адаптивным выбором СП в виде (3) и использованием сосредоточения дает в данном примере более точные результаты, чем методы, предложенные в работах [16, 19].
Пример 5. Рассмотрим начально-краевую задачу при х є [a; b] для одномерного уравнения Бюргерса (7) при 1 = 1 и f = 0 с известным аналитическим решением
. ч a + m + (m-a)eh (х -mt -B)a w ™
u (t, x) =------------, где h =----------. Данная задача применялась в [16, 20, 21]
а для метода работы [19] - errmax » 2,00 10
-5
1 + e
h
к
для проверки и тестирования предложенных там численных методов решения уравнения Бюргерса. Ее решение представляет собой движущуюся "ступеньку" (решение типа
"кинк"). Положим, а = 0,4, b = 0,125, m = 0,6, к = 10 2, t = 10 2, h = 1/36 [16, 20]. Тогда при использовании (3) с сосредоточением (подробный выбор СП описан в примере 3) в фиксированный момент времени t = 1/2 при различных 0 = 0, 0 = 1/4, 0 = 1/2, 0 = 3/4 и 0 = 1 имеем соответственно errmax » 0,03693, errmax » 0,03045, errmax » 0,03239,
0,04670. Аналогично, при использовании (4) и 0(П+1'), приобретаю-
-2
errmax 0,04, errmax
щем значения 0, 1/4,
1/2, 3/4 и 1 соответственно, имеем errmax »0,005.
errmax » 0,00409, errmax » 0,00604, errmax » 0,00928, errmax » 0,01246. Поскольку в данной
задаче нет сильного доминирования конвекции (число Пекле имеет порядок 10 ) и резких неоднородностей в решении, то, соответственно, нет нужды во введении сильной сглаживающей диссипации, поэтому наиболее эффективными оказались "мягкие" варианты выбора СП (в особенности, (4)) и при небольших значениях параметров 0 и 0(П+1'). Для метода работы [16] при тех же условиях расчета имеем errmax » 0,00597, а для конечноэлементного метода работы [20] - errmax » 0,045 [16, 20].
-3
Положим, теперь h = 1/18 и t = 10 (такие значения параметров использовались в [21]). Тогда при использовании (4) и 0(П+1\ приобретающем значения 0, 1/4, 1/2, 3/4 и
1 соответственно,
имеем
errmax » 0,02494,
errmax » 0,01347,
errmax » 0,02165 ;
errmax » 0,03432, errmax » 0,0481. При указанных значениях параметров для конечноэле-
ISSN 1028-9763. Математичні машини і системи, 2014, № 4
149
ментного метода SGA (Standard Galerkin Approach) имеем errmax » 0,096, для конечноэлементного метода PAG (Product Approximation Galerkin) errmax »0,082, для метода CD (Compact finite Difference) errmax » 0,151 (описание данных методов и указанные значения errmax в [20, 21]). Таким образом, расчеты показывают, что и в данном примере предлагаемая версия МПГ с адаптивным выбором СП (3)-(4) способна дать более точные результаты, чем методы, предложенные в работах [16, 20, 21].
Пример 6. Рассмотрим начально-краевую задачу для одномерного уравнения Бюр-герса (7) при 1 = 1 и f = 0 на отрезке [0; 1] с известным аналитическим решением
, ч kJu2 -1 кsh(x-кt) ^
u(t,x) = к +-----------------------, где m - некоторая константа. Данное решение
ch( x -к t) + m ch( x -к t) + m
имеет солитонный тип [22] и применялось в работах [16, 22] для проверки и тестирования предложенных там численных методов. Положим, к = 75 10 , m = -3, h = 10 , t = 10
(при таких условиях задача использовалась в [16, 22]). Тогда при использовании (3) (c 0 = 1, пример 3) с сосредоточением в фиксированный момент времени t = 1/20 имеем
errmax »10-13. При использовании (4) ( 0(Ид-+1) = 1) имеем errmax » 2 • 10-9 . Для метода рабо______________________________о
ты [16] получаем errmax » 2,1 • 10 (табл. 10.1 в [16]), а для 12-точечной схемы работы [22] - errmax » 1,058 • 10-4 (табл. 6 в [22]).
Пример 7. Рассмотрим начально-краевую задачу для уравнения (7) при 1 = 1 и f = 0 на отрезке [0; 1] с нулевыми граничными условиями первого рода и начальным условием U0 (x) = sin(px). Отметим, что задача в такой постановке дает описание нелинейных стоячих медленных волн в корональных магнитных петлях [23] в солнечной плазме. Аналитическое решение начально-краевой найдено в [20, 23]. Возьмем сначала значение
к = 10-4. В силу необходимости выполнения нулевого граничного условия вблизи точки x = 1 возникает тонкий приграничный слой. Сравним погрешность предлагаемой версии МПГ с погрешностями методов, использованных в работе [20] при t = 1/2 на наборе точек x = 0,5, x = 0,556, x = 0,611, x = 0,667, x = 0,722, x = 0,778, x = 0,833, x = 0,889, x = 0,944 (именно такой тестовый набор точек был взят в [20]). Взяв h = 1/216 и t = 1/8 [20] для конечноэлементного метода работы [20], получим errmax »0,097, для конечноэлементного метода работы [24] с кубическими В-сплайнами получим errmax » 0,078 (табл. 3 в [20]), а для предлагаемой версии МПГ с (3) (при 0 = 1) и использования сосредоточения получим
-3
errmax »0,022 (при 0 = 1,3 получим errmax »0,02). Положим, теперь h = 1/18, t = 10 . Тогда для метода SGA получим errmax » 0,124, для метода PAG - errmax » 0,105, а для метода CD - errmax » 0,087 (по поводу данных методов см. табл. 3 в [20]). Для МПГ с (3) (при
0 = 1) и сосредоточением получим errmax » 0,0617. Для МПГ с сосредоточением и постоянными СП, равными единице, получим errmax » 0,063 (а решение получается чрезмерно сглаженным возле приграничного слоя), а без использования сосредоточения получим errmax » 0,119. Для МКЭ Галеркина решение получается колебательным вблизи приграничного слоя и имеет большую погрешность (с использованием сосредоточения будем иметь errmax » 0,839, а без его использования - errmax » 0,975 ).
Пример 8. Рассмотрим начально-краевую задачу для уравнения Бюргерса (7) при
1 = -1 и f = 0 на отрезке [0;1] с известным аналитическим решением
150
ISSN 1028-9763. Математичні машини і системи, 2014, № 4
u(t, x) :
sh( x /(2 k))
ch( x /(2k)) + exp(-t /(4 k))
. Данная задача использовалась в [25] для тестирования
__5 __2
предложенных там неявных численных схем. Положим, к = 75 10 , h = 1/20, t = 10
(при таких условиях проводились расчеты в [25]). Тогда при использовании (3) (0 = 1) с
_2
сосредоточением в фиксированный момент времени t = 1/20 имеем errmax » 5 • 10 . При
использовании (4) ( 0(Ид-+1) = 1) имеем errmax » 6,5 • 10_2 . Для 8-точечной неявной схемы ра-
боты [25] имеем errm
_2
7 • 10 , а для симметричной схемы Кранка-Николсона -
err
5 • 10 1 (табл. 4 в [25]).
Пример 9. Рассмотрим начально-краевую задачу для двухмерного уравнения Бюр-герса [6, 26-28]:
ди » ди » ди
——+ K1U—---+ 12U
д
dt
dx
dx2 dx
K-
ди
dx
+ -
д
1J
dx0
du
dx0
K-
V ил2
+ f
(где 1 =1i (t, x), K = K(t, x), f = f(t, x)) при І1 = 12 = 1, к = 1/20, f = 0,
(xb x2) є [0;1] x [0;1] с известным аналитическим решением
u(t,x1,x2) = (1 + exp((x1 + x2 _t)/(2k)))_1 [26-28]. Данное решение служит моделью ударной волны, волновой фронт которой движется вдоль линии x1 + x2 _ t (интерпретацию решения и его гидродинамическую суть см. в [26, 27]). Используем равномерное разбиение 21 x 21 узлов и положим t = 10 3. Тогда, при использовании (4) (с 0(”k+1) = 0,2) для моментов времени t = 1/2, t = 3/4, t = 1, t = 5/4, имеем соответственно errmax »0,00032, errmax » 0,00097, errmax » 0,0007, errmax » 0,00056. Отметим, что увеличение или уменьшение 0("k+1), как показывают расчеты, ведет к (медленному) росту ошибки (по сравнению
с ошибкой для 0(”к+1^ = 0,2 ); это, как уже отмечалось ранее, связано с тем, что коэффициент вязкости к не является малой величиной, следовательно, конвекция в данной задаче не является доминирующей, поэтому небольшие значения коэффициентов 0(n+1) приводят
к более точным решениям (однако при нулевых значениях 0(П+1') погрешность увеличивается). Для meshless-метода ELMFS (Eulerian-Lagrangian Method of Fundamental Solutions) работы [26] при тех же условиях интегрирования задачи имеем соответственно errmax » 0,00061, errmax » 0,00105, errmax » 0,00299, errmax » 0,00303 (табл. 1 в работе [26]
_2
на стр. 399). Рассмотрим теперь задачу при разбиении 20x20 узлов и положим t = 10 (при таких условиях задача рассматривалась в [27]). Тогда на промежутке 0 < t < 5/4 для МПГ с (4) имеем errmax » 0,007 . Для методов MFS-DRM (Method of Fundamental Solution -Dual Reciprocity Method) и Kansa работы [27] при тех же условиях решения имеем соответственно errmax » 0,0703 и errmax » 0,0719 (стр. 218 в [27]). Пусть теперь к = 1, t = 10_4 и
используется разбиение 11 x 11 узлов. Тогда при t = 1/2 для МПГ с (4) (0(Ид+1') = 1) имеем
_7
errmax » 8,72 10 , а для метода LMAPS (Local Method of Approximate Particular Solutions)
работы [28] - errmax » 6,9 • 10_6 (табл. 1 в [28]; отметим, однако, что авторы [28] не указали в своей работе ни метода интегрирования своих полудискретных аппроксимаций, ни кон-
ISSN 1028-9763. Математичні машини і системи, 2014, № 4
151
кретного значения шага по времени, с которым они решали данный пример). Таким образом, предлагаемая версия МПГ с выбором СП в виде (4) в данной задаче дает более высокую точность, чем методы работ [26-28]. Отметим, что для многомерных задач использование сосредоточения является "рискованным" (в силу того, что может ввести диффузию поперек потока [1], которая приведет к возникновению больших погрешностей), а вариант (3) выбора СП без сосредоточения показывает плохие результаты даже в одномерном случае (пример 3), поэтому при решении многомерных задач вариант (4) выбора СП является предпочтительным.
Пример 10. Рассмотрим задачу для системы Бюргерса [26, 28]:
Эи , Эи , Эщ
—+ ----+12Щ2
Э ( Эи
\
dt
Эщ
dt
2 + I1U1
dx1
Эи
Эх
2
к
Эх1 ^ Эх1
+ -
Э
Эх2
к
Эи
Эх
2 J
Эх
2 Эи2
+ I2U2 —“
Э
/
Эх2 Эх
к
Эи
2
Эх
+ -
Э
1J
Эх
к
Эи
\
Эх
2J
при 11 =12 = 1, к = 1/100, (х1, х2) є[0;1] х [0;1] с аналитическим решением [26]
u1(t, х1, х2) = 3/4 - (1 + exp((-4х1 + 4х2 -1) / (32к)))-1 / 4, u2(t,х1,х2) = 3/4 + (1 + exp((-4х1 + 4х2 -1)/(32к)))-1 /4.
Данное решение служит моделью ударной волны (гидродинамическую интерпретацию решения см. в [26]) и использовалось в работе [26] для проверки работоспособности предложенного там метода ELMFS (пример 9). Используем равномерное разбиение 11 х 11
узлов и положим t = 0,005. Тогда, при использовании (4) (с 0г(П+1) = 0,2, см. комментарии
к примеру 9) для моментов времени t = 10 , t = 1/2, t = 2 для U1 имеем соответственно
errmax »3,07 10 , errmax »7,32 10 , errmax »7,54 10 , для величины u2 погрешности такие же (с точностью до малых более высокого порядка). Для метода ELMFS работы [26] при тех же условиях интегрирования задачи соответствующие погрешности для U1 имеют
-4 -3 -2
вид errmax » 5,39 10 , errmax » 7,66 10 , errmax »1,03 10 (для величины и2 они будут
такими же, см. табл. 3 и 4 в работе [26] на стр. 404-405). Для метода LMAPS работы [28]
-3
(см. комментарии к примеру 9) при тех же условиях и t = 1/2 имеем errmax » 7,4 10
(табл. 2 в [28]). Таким образом, предлагаемая версия МПГ с выбором стабилизирующих параметров в виде (4) в данной задаче дает более высокую точность, чем методы работ [26, 28].
Пример 11. Рассмотрим важный случай системы уравнений магнитной гидродинамики (в предположении адиабатичности и бесконечной проводимости среды [7, 29]), в которой все величины зависят от одной пространственной координаты, например, координаты z [29]. Тогда имеем вектор скорости и = (их(t,z),Uy(t,z),Uz(t,z)), магнитную индукцию B = (Вх (t, z), B (t, z), Bz (t, z)), плотность p = p(t, z) и давление p = p(t, z). К
этой ситуации сводится математическое описание распространения волн при возмущении какой-либо из МГД величин на границе рассматриваемой области, а также движение плазмы в конфигурациях, которые могут быть описаны (по крайней мере, в некотором приближении) единственной пространственной координатой (например, явление спикул [29]).
Введем сетку zt, i = 1..n на отрезке z є [L1;L2], z1 = zn = L2 . Начальные усло-
вия, соответствующие стационарному равновесному состоянию плазмы, имеют вид
152
ISSN 1028-9763. Математичні машини і системи, 2014, № 4
u = u00 = 0, B = Bq = 0, p = Pq = Cp0 = Cp (где C — некоторая константа),
p = pо (z) = Po(0)e_gz/C , где сила тяжести (направлена противоположно оси z ) g = const. Соответственно, po = Cpo. По компоненте скорости uz на левом конце области дадим следующее синусоидальное возмущение:
duz (t)
sin(2pt /(3/20)), 0 £ t £ 3/20, 0, t > 3/20.
Без ограничения общности, для демонстрации работоспособности обсуждаемой версии МПГ положим при расчетах C = 1/2, pq(0) = 1, g = 1 и коэффициент вязкости _______2
и = 10 [7]. Тогда при использовании МПГ с выбором СП в виде (3) с сосредоточением,
_2 _2
0 = 1 (пример 3), h = 10 (используем плавающую сетку работы [30]) и t = 10 графики
плотности p и компоненты uz скорости в момент времени t = 0,15 в зависимости от координаты z изображены на рис. 16 и 17 соответственно. На рис. 18 и 19 показаны графики соответствующих величин в момент времени t = 0,25 .
Рис. 16. График плотности p при t = 0,15 Рис. 17. График скорости uz при t = 0,15
Рис. 18. График плотности p при t = 0,25 Рис. 19. График скорости uz при t = 0,25
Видно, что возмущение duz (t) сохраняет свою форму, но с течением времени затухает по амплитуде, что объясняется действием диссипации (вязкость с коэффициентом и) и нелинейности конвективного члена, причем, чем большим будет значение коэффициента вязкости u , тем сильнее будет затухать амплитуда скорости с течением времени. Возмущение скорости приводит к возмущению плотности (и, в свою очередь, давления) на фоне
равновесного установившегося распределения плотности po(z) = po(0)e_gz/C . Как видно из рис. 16 и 18, происходит "скатывание" возмущения со стороны левого конца области. С течением времени данное возмущение затухает (поскольку затухает возмущение переменной скорости) и плазма, при отсутствии других внешних возмущений, постепенно переходит в равновесное состояние. Таким образом, полученные результаты расчетов полностью согласуются с физической картиной поведения рассматриваемого процесса.
ISSN 1028-9763. Математичні машини і системи, 2014, № 4
153
5. Заключение
В настоящей статье получили дальнейшее развитие результаты работ [5-7], касающиеся выбора и динамической (адаптивной) настройки стабилизирующих параметров в конечноэлементном методе Петрова-Г алеркина с предложенными в работах [2-4] весовыми (поверочными) функциями. На многочисленных примерах (разд. 4) продемонстрирована высокая эффективность метода Петрова-Галеркина с весовыми функциями (1) и предложенными вариантами выбора стабилизирующих параметров (разд. 3) при решении различных задач конвекции-диффузии (стационарных и нестационарных линейных задач, а также нелинейных уравнений и систем типа Бюргерса). Во всех приведенных примерах рассматриваемый метод Петрова-Галеркина оказался способным дать более точное решение, чем методы работ [8, 13-17, 19-22, 24-28], в которых данные примеры (те или иные) также применялись для проверки точности и работоспособности предложенных там методов. Расчеты показали, что в целом, в случае одномерных задач с наличием тонких слоев и ударных волн (особенно задач для уравнений Бюргерса), среди вариантов (3) и (4) выбора стабилизирующих параметров более эффективным (в плане точности и подавления осцилляций), как правило, оказывается "жесткий" вариант выбора (3) с одновременным использованием сосредоточения. В случае же, когда конвекция не является доминирующей над процессом диффузии, а также для многомерных задач, "мягкий" вариант (4) выбора стабилизирующих параметров в приведенных примерах оказывается более эффективным по сравнению с вариантом (3).
Благодарности. Автор выражает глубокую благодарность к.т.н., с.н.с. ИКИ НАНУ Сальникову Н.Н. и д.т.н., проф. Молчанову А.А. за плодотворное обсуждение работы и ценные замечания к ней.
СПИСОК ЛИТЕРАТУРЫ
1. Roos H.-G. Robust numerical methods for singularly perturbed differential equations / H.-G. Roos, M. Stynes, L. Tobiska. - Berlin, Heidelberg: Springer-Verlag, 2008. - 604 p.
2. Сальников Н.Н. О построении конечномерной математической модели процесса конвекции-диффузии с использованием метода Петрова-Галеркина / Н.Н. Сальников, С.В. Сирик, И.А. Терещенко // Проблемы управления и информатики. - 2010. - № 3. - С. 94 - 109.
3. Сальников Н.Н. Построение весовых функций метода Петрова-Галеркина для уравнений кон-векции-диффузии-реакции в трехмерном случае / Н.Н. Сальников, С.В. Сирик // Кибернетика и системный анализ. - 2014. - № 5. - С. 173 - 183.
4. Сирик С.В. Выбор весовых функций в методе Петрова-Галеркина для интегрирования линейных одномерных уравнений конвекции-диффузии / С.В. Сирик, Н.Н. Сальников, В.К. Белошапкин // Управляющие системы и машины. - 2014. - № 1. - С. 38 - 47.
5. Сирик С.В. Численное интегрирование уравнения Бюргерса методом Петрова-Галеркина с адаптивными весовыми функциями / С.В. Сирик, Н.Н. Сальников // Проблемы управления и информатики. - 2012. - № 1. - C. 94 - 110.
6. Молчанов А.А. Выбор весовых функций в методе Петрова-Галеркина для интегрирования двухмерных нелинейных уравнений типа Бюргерса / А.А. Молчанов, С.В. Сирик, Н.Н. Сальников // Математичні машини і системи. - 2012. - № 2. - С. 136 - 144.
7. Сальников Н.Н. О построении конечномерных математических моделей для двухмерных процессов магнитной гидродинамики с использованием метода Петрова-Галеркина / Н.Н. Сальников,
С.В. Сирик, В.К. Белошапкин // Управляющие системы и машины. - 2014. - № 4. - С. 23 - 32.
8. Harari I. Streamline design of stability parameters for advection-diffusion problems / I. Harari, L.P. Franca, S.P. Oliveira // Journal of Computational Physics. - 2001. - Vol. 171 (1). - P. 115 - 131.
9. John V. Finite element methods for time-dependent convection-diffusion-reaction equations with small diffusion / V. John, E. Schmeyer // Comput. Methods Appl. Mech. Engrg. - 2008. - Vol. 198. - P. 475 -494.
10. Хорн Р. Матричный анализ / Р. Хорн, Ч. Джонсон; пер. с англ. - М.: Мир, 1989. - 655 с.
11. Калиткин Н.Н. Численные методы / Калиткин Н.Н. - М.: Наука, 1978. - 512 с.
154
ISSN 1028-9763. Математичні машини і системи, 2014, № 4
12. Скворцов Л.М. Простые явные методы численного решения жестких обыкновенных дифференциальных уравнений / Л.М. Скворцов // Вычислительные методы и программирование. - 2008. -Т. 9. - С. 154 - 162.
13. Brooks A.N. Streamline upwind Petrov-Galerkin formulations for convection dominated flows with particular emphasis on incompressible Navier-Stokes equations / A.N. Brooks, T.J.R. Hughes // Computer Methods in Applied Mechanics and Engineering. - 1982. - Vol. 32 (1 - 3). - P. 199 - 259.
14. Xie S.-S. Numerical solution of one-dimensional Burgers' equation using reproducing kernel function /
S.-S. Xie, S. Heo, S. Kim [et al.] // J. Comput. Appl. Math. - 2008. - Vol. 214. - P. 417 - 434.
15. Dag I. B-spline collocation methods for numerical solutions of the Burgers equation / I. Dag, D. Irk, A. Sahin // Math. Problems in Eng. - 2005. - Vol. 5. - P. 521 - 538.
16. Mittal R.C. Numerical solutions of nonlinear Burgers equation with modified cubic B-splines collocation method / R.C. Mittal, R.K. Jain // Applied Mathematics and Computation. - 2012. - Vol. 218. -P.7839 - 7855.
17. Korkmaz A. Polynomial based differential quadrature method for numerical solution of nonlinear Burgers equation / A. Korkmaz, I. Dag // Journal of the Franklin Institute. - 2011. - Vol. 348 (10). -P.2863 - 2875.
18. Сирик С.В. О применении сосредоточения в методе конечных элементов Петрова-Галеркина при решении задач конвекции-диффузии / С.В. Сирик, Н.Н. Сальников // Доповіді НАН України. -2014. - № 5. - С. 39 - 44.
19. Ganaie I.A. Numerical solution of Burgers equation by cubic Hermite collocation method / I.A. Ga-naie, V.K. Kukreja // Applied Mathematics and Computation. - 2014. - Vol. 237. - P. 571 - 581.
20. Dogan A. A Galerkin finite element approach to Burgers' equation / A. Dogan // Applied Mathematics and Computation. - 2004. - Vol. 157. - P. 331 - 346.
21. Christie I. Product approximation for nonlinear problems in the finite element method / I. Christie,
D.F. Griffiths, A.R. Mittchell [et al.] // IMA. J. Numer. Anal. - 1981. - Vol. 1. - P. 253 - 266.
22. Hesameddini E. Soliton and numerical solutions of the Burgers Equation and comparing them / E. He-sameddini, R. Gholampour // Int. J. Math. Anal. - 2010. - Vol. 4 (52). - P. 2547 - 2564.
23. Ruderman M.S. Nonlinear damped standing slow waves in hot coronal magnetic loops / M.S. Ruder-man // Astronomy & Astrophysics. - 2013. - Vol. 553 (A23). - P. 1 - 6.
24. Ali A.H.A. A collocation solution for Burgers equation using cubic B-spline finite elements / A.H.A. Ali, G.A. Gardner, L.R.T. Gardner // Comput. Meth. Appl. Mech. Eng. - 1992. - Vol. 100. - P. 325 - 337.
25. Tabatabaei A.H. Some implicit methods for the numerical solution of Burgers equation / A.H. Tabata-baei, E. Shakour, M. Dehghan // Applied Mathematics and Computation. - 2007. - Vol. 191. - P. 560 -570.
26. Young D.L. The Eulerian-Lagrangian method of fundamental solutions for two-dimensional unsteady Burgers' equations / D.L. Young, C.M. Fan, S.P. Hu [et al.] // Engineering Analysis with Boundary Elements. - 2008. - Vol. 32. - P. 395 - 412.
27. Li J. Numerical comparisons of two meshless methods using radial basis functions / J. Li, Y.C. Hon, C.S. Chen // Engineering Analysis with Boundary Elements. - 2002. - Vol. 26. - P. 205 - 225.
28. Zhang X. Local method of approximate particular solutions for two-dimensional unsteady Burgers' equations / X. Zhang, H. Tian, W. Chen // Computers and Mathematics with Applications. - 2014. -Vol. 66. - P. 2425 - 2432.
29. Ладиков-Роев Ю.П. Математические модели сплошных сред / Ю.П. Ладиков-Роев, О.К. Черем-ных. - К.: Наукова думка, 2010. - 551 с.
30. Молчанов О.А. Проекційні методи інтегрування рівнянь магнітної гідродинаміки / О.А. Молчанов, М.М. Сальніков, С.В. Сірик // ІІІ наукова конференція магістрантів та аспірантів «Прикладна математика та комп'ютинг»: тези доповідей, (Київ, 13-15 квітня 2011). - К.: Просвіта, 2011. - С. 272 - 278.
Стаття надійшла до редакції 21.10.2014
ISSN 1028-9763. Математичні машини і системи, 2014, № 4
155
