ДОКЛАДЫ АКАДЕМИИ НАУК, 2012, том 446, № 5, с. 504-509
ИНФОРМАТИКА
УДК 519.6
МОНОТОННАЯ БИКОМПАКТНАЯ СХЕМА ДЛЯ КВАЗИЛИНЕЙНЫХ УРАВНЕНИЙ ГИПЕРБОЛИЧЕСКОГО ТИПА
© 2012 г. Б. В. Рогов
Представлено академиком Б.Н. Четверушкиным 05.03.2012 г. Поступило 05.03.2012 г.
Уравнения и системы уравнений гиперболического типа составляют значительную часть математических моделей, используемых для решения разнообразных прикладных задач [1]. Из них ряд задач связан с нахождением обобщенных решений, которые описываются кусочно-гладкими функциями. Для нахождения численных решений нелинейных задач широкое распространение получили разностные схемы сквозного счета повышенного порядка аппроксимации. Важнейшими свойствами, которыми должны обладать такие схемы, являются монотонность [1] и консервативность [2], поэтому построению монотонных и консервативных разностных схем посвящена обширная литература [1, 2].
Ранее в работе [3] для квазилинейных уравнений и систем уравнений гиперболического типа предложены бикомпактные разностные схемы, имеющие на двухточечном шаблоне четвертый порядок аппроксимации по пространственной координате. Эти схемы являются абсолютно устойчивыми, консервативными и монотонными в области локальных чисел Куранта к > 1 [4]. Кроме того,
4
они являются экономичными и могут быть решены методом бегущего счета.
В настоящей работе для квазилинейных уравнений гиперболического типа предложена новая монотонная бикомпактная схема повышенного порядка аппроксимации, которая так же, как и схемы [3], является абсолютно устойчивой, консервативной и экономичной. Однако в сравнении с высокоточной схемой [3] имеются следующие преимущества. Во-первых, она несколько расширяет область монотонности схемы [3] в сторону меньших значений числа Куранта. Во-вторых, предлагаемая схема является более надежной при решении жестких задач, чем схема [3]. Если в схе-
Институт прикладной математики им. М.В. Келдыша Российской Академии наук, Москва Московский физико-технический институт (государственный университет), Долгопрудный Московской обл.
ме [3] используется ^-устойчивая схема интегрирования по времени, то в новой схеме используется Х-устойчивая схема. Следует также отметить, что в монотонной компактной схеме данной работы не использованы искусственная вязкость (см., например, [5]) и какие-либо монотонизато-ры в виде ограничителей численных потоков (см., например, [6]). Приведены результаты расчетов, демонстрирующие точность предложенной схемы, ее консервативность и монотонность при решении тестовых задач для квазилинейного уравнения Хопфа.
1. БИКОМПАКТНАЯ СХЕМА ПЕРВОГО ПОРЯДКА АППРОКСИМАЦИИ ПО ВРЕМЕНИ
Сначала рассмотрим смешанную задачу Коши [7] для одного квазилинейного гиперболического уравнения, записанного в дивергентном виде:
Ut + f (u) = 0, a (u) = df(u) > 0,
x du (1)
u(x, 0) = u0(x), x > 0; u(0, t) = ц(0, t > 0. Разностную схему для начально-краевой задачи построим с помощью метода прямых и интегро-интерполяционного метода [8]. Для этого введем неравномерную сетку {x, j > 0} на интервале [0, да). На временном слое t = const проинтегрируем уравнение (1) на отрезке x, < x < xj +1, в результате получим дифференциальное по t и разностное по x уравнение
(uj+i)t =- f [+1 - fj ] fj = fj fli
(2)
"j+i
hj+i = Xj+i - Xj, j > 0, где u;+1 — интегральное среднее функции u(x, t) по ячейке [ху-, Xj+1]:
uj+i =
-M j}
udx.
(3)
В дальнейшем индекс при пространственном шаге h опущен.
Вычислим интегральное среднее от функции f(u) на отрезке Xj < x < Xj + 1
fj+i = h J f (u)dx
(4)
двумя различными способами с точностью до 0(к4): по формуле Симпсона [7]
7+1 = 6 (( + 4+1/2 + /) + 0(к\ (5)
7+1/2 = 7 (и]+1/2 ) и по формуле Эйлера—Маклорена [7]
7,1 = 7 + Л)-^'!- + 0(*4)- (6)
Из формулы (5) исключим значения искомой функции в полуцелых узлах с помощью формулы Симпсона для интегрального среднего (3)
и+1 =1 ((+1 + 4и;+1/2 + иу) + 0(к4), (7)
6
а из формулы (6) — значения производной/х в целых узлах сетки с помощью уравнения (1). Затем приравняем преобразованные выражения (5) и (6) друг другу. В результате получим дифференциальное по I и разностное по х уравнение
к +1- и), =
4
fj+i - 2f (2 uJ+i - 4 uJ+i - 4 uJ) + fj
J > 0, (8)
аппроксимациями, получим семейство бикомпактных разностных схем для численного решения уравнения (1). Например, если производные по времени на расчетном слое I = 1п + 1 аппроксимировать разностями назад, то получим разностную схему
+ r
Uj+i + r[f (uj+i) - f (uj)] = Uj+i, j > 0, 1 (uj+i - uj ) +
f (uj+i) - 2f (2uj+i - 4uj+i - 4uj) + f (uj)
i / n n\
=4 (uj+i- uj )•
(9)
(10)
выполняющееся с точностью до 0(к4) на точном решении уравнения (1).
Нетрудно видеть, что уравнение (8) есть разностный аналог дифференциального следствия уравнения (1)
(и, + /х) х = 0 которое так же, как и уравнение (1), имеет дивергентный вид.
Приведенные выше выкладки для скалярного уравнения (1) можно полностью повторить и в случае системы уравнений
и, + Г (и) х = 0
где и(х, I) — искомая вектор-функция с т компонентами, 1(и) — заданная вектор-функция размерности т. В результате можно получить системы обыкновенных дифференциальных уравнений (ОДУ), аналогичные (2), (8).
По терминологии работ [9, 10] дифференциально-разностная схема (2), (8) является бикомпактной (двухточечной компактной). Она может быть решена методом бегущего счета подобно схеме [3]. Кроме того, данная схема сохраняет порядок точности при переходе от равномерной сетки к неравномерной.
Введем неравномерную сетку {1п, п > 0} на интервале [0, да). Заменяя производные по времени в системе ОДУ (2), (8) какими-либо разностными
Здесь верхний индекс п + 1 при значениях сеточных функций на рассчитываемом слое I = 1п + 1
опущен; г = т, т = 1п + 1 — 1п — временной шаг. Схема к
(9), (10) есть неявная схема Эйлера, которая имеет первый порядок аппроксимации по времени и является Х-устойчивой [11], а следовательно, и абсолютно устойчивой. Если исключить из уравнения (10) значение +1 с помощью уравнения (9), то получаем схему бегущего счета для расчета искомой функции и в целых узлах.
2. МОНОТОННОСТЬ
Исследуем монотонность схемы (9), (10) сначала в случае, когда функцияf(u) является линейной:
f (u) = au, a = const > 0. (11)
Уравнения (9), (10) с учетом (11) примут следующий вид:
ku,-+i + u,+i - ku,- = uj+i, j > 0,
(12)
(3 K +1) uJ+i - 3kuj+i +
+ (3к-1)и;. = 1 (и;+1 - и;), ] > 0, (13)
где к = аг — число Куранта.
Анализ монотонности схемы (12), (13) проведем, следуя методике классической работы [12]. Полагая, что любую ограниченную монотонную сеточную функцию можно представить как линейную комбинацию простейших монотонных ступенчатых функций, рассмотрим свойство монотонности (по Годунову [12]) схемы (12), (13) на этих функциях.
Рассмотрим задачу (1) с монотонно убывающей начальной функцией, имеющей вид ступеньки,
u 0(x) = \
0 < x < i, [0, x > i и граничным условием
= 2, t > 0.
(14)
x
x
Для упрощения последующих выкладок положим, что пространственная сетка — равномерная.
Дифференциальная задача (1), (14), (15) с разрывным начальным условием аппроксимируется разностной схемой (12), (13) со следующими начальными и граничными условиями:
и
[2, 0 < ] < л,
0, У > 7о,-[2, 1 < У < Уо,
1, У = У о +1, .0, у > Уо +1,
и0 = 2, п > 0.
(16)
иу =
2, 0 < у < У0,
2Р1, У = Уо + 1,
Р2и)-1, У > Уо + 1,
2
Р1
К
2,1 , 1 ' к + -кн—
2 12
Р2
2 1,1
к —кн—
2 12
2,1 , 1 ' к + -кн—
2 12
(17)
2р2,
рз
р2иу--1, У > Уо + 2, 3_24
(18)
к 2 + 1 к + 1.
к2 + 1 к + -1'
2 12
которая также монотонно убывает при любом к > 0.
Условием монотонности решения на слое п = 1 является неравенство
1.-1.1 и У _1 > и У > и У,
которое с учетом формул (17), (18) приводит к ограничению на число Куранта
к>
(19)
Здесь х к = 1; значения сеточной функция и у
определяются по формуле (3) с и = и0(х). При этом интеграл в (3) рассчитывается либо аналитически, либо по квадратурной формуле точности 0(к5) или более высокой на разностных интервалах, где функция ы°(х) непрерывна. В разностных ячейках, где эта функция является разрывной, интеграл рассчитывается по квадратурной формуле трапеций.
На одном временном шаге разностные уравнения (12), (13) переведут монотонно убывающее
начальное условие и0 в сеточную функцию
Аналогичные выкладки для (12), (13) показывают, что эта разностная схема также переводит монотонно возрастающее начальное условие в форме элементарной ступенчатой функции в монотонно возрастающую сеточную функцию при выполнении условия (19). Численный анализ показал, что и при дальнейшей эволюции разностных решений, рассчитываемых по схеме (12), (13), их монотонность сохраняется.
Исследование монотонности схемы (9), (10) в случае произвольной нелинейной функции /(и) затруднительно, поэтому ограничимся анализом свойств монотонности схемы (9), (10) для уравне-
2
ния (1) с /(и) = —, т.е. для уравнения Хопфа [13]:
иI + а(и)их = 0, а(и) = и > 0.
(20)
В уравнении (20) величина а(и) является локальной скоростью переноса. Этой скорости соответствует локальное число Куранта к(и) = га(и). Заменяя нелинейное уравнение (20) линейным в приближении "замороженных коэффициентов", можно воспользоваться формулой (19) для оценки свойств монотонности схемы (9), (10) при решении уравнения (20). Для уточнения оценки (19) с к(и) = ги для нелинейного уравнения (20) были проведены численные эксперименты.
На рис. 1 показаны результаты расчета смешанной задачи Коши для уравнения (20) со следующими начальными и граничными условиями:
которая монотонно убывает при любом к > 0. В то же время й0 перейдет в сеточную функцию
2, 1 < У < Уо,
2рз, У = Уо + 1, У = Уо + ^
'(*) =
1, 0 < х < 1,
2 - х, 1 < х < 2 -5, |а (I) = 1, (21) 5, х > 2 - 5.
У > 1,
Точное решение этой задачи таково, что при I = 1 формируется разрыв, который затем движется с постоянной скоростью Б = 0.5(1 + 8) в положительном направлении оси х [13]. Расчеты проведены на сетке с отношением шагов г = 1. В каждый момент времени локальное число Куранта к(и) = ги изменяется от 8 до 1. Значения решений в целых узлах показаны на
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.