CC BY
62
Вычислительные технологии
Том 7, № 1, 2002
О НЕЯВНЫХ АБСОЛЮТНО УСТОЙЧИВЫХ СХЕМАХ РУНГЕ-КУТТЫ ЧЕТВЕРТОГО
ПОРЯДКА*
В. И. Пинчуков
Институт вычислительных технологий СО РАН, Новосибирск, Россия
e-mail: pinchvi@ict.nsc.ru
Implicit fourth order schemes of Runge — Kutta type for ordinary difference equations and for multidimensional transfer equation are considered. The absolute stability of schemes is proved under the assumption about frozen coefficients. The results of test calculations are presented.
В настоящее время методы повышенной точности доказали свою эффективность в вычислительной аэродинамике. Помимо классических явных схем третьего порядка типа Рунге — Кутты — Русанова [1], Бернстейна — Мирина [2], Кутлера — Ломекса [3] сконструированы современные адаптивные ENO-схемы третьего и более высоких порядков (см. [4, 5] и обзор [6]). В [7] предложены компактные неявные условно устойчивые схемы третьего порядка, а также схемы более высоких порядков. Абсолютно устойчивые схемы Рунге — Кутты третьего порядка развиты в [6, 8-11]. Здесь предлагаются неявные абсолютно устойчивые схемы четвертого порядка.
1. Схемы для обыкновенных дифференциальных уравнений
Рассмотрение начнем со случая обыкновенных уравнений. Пусть решается задача Коши для уравнения
dv
— = ^(М), v(0) = vo.
Предлагаемые схемы могут быть получены модификацией известных явных схем четвертого порядка. В качестве базовой рассмотрим следующую схему:
v(1) - vn = <^(vn,t)
Т = 2 '
v(2) - vn = y(v(1),t + т/2)
Т = 2 '
* Работа выполнена при частичной финансовой поддержке Российского фонда фундаментальных исследований, грант №02-01-00922, и Civilian Research & Development Foundation, проект №10086. ©В. И. Пинчуков, 2002.
v(3)_vn
<^(v(2),i + т/2),
vn+i _ vn
<^(vn,i) ^(v(1),t + т/2) <^(v(2) ,t + т/2) ^(v(3),t + т)
+
+
+
т 6 3 3 6
Второе уравнение заменим на его разность с первым уравнением. Для получения неявной схемы будем использовать на втором и четвертом шагах схемы стабилизирующие операторы
v
(1) - vn ^(vn,t)
т
[1 + А(Т^П)2 ]
(2) - v(1) ^(v(1),t + т/2) - ^(vn,t)
v - v т
v(3) vn т
= ^(v(2),t + т/2),
(1) (2) (3)
[1 + £ (т^П)4]
vn+1 - vn ^(vn,t) ^(v(1),t + т/2) ^(v(2),t + т/2) ^(v(3),t + т)
6
+
3
+
3
+
6
(4)
где £ > 0, Л > 0 — числовые параметры, которые следует определить в анализе. Для исследования устойчивости и точности этой схемы рассмотрим случай <^(v,t) = - Av, A = const > 0. Исключая промежуточные временные слои, получаем
vn+1 = vnP^A),
где
Р(х)
1
(2 - х)х2
1---1---+
2 12 24(1 + Лх2) J
х
(1 + £x4)'
(5)
Чтобы иметь четвертый порядок, необходимо выполнение условия р(х) — ехр(—х) = 0(х5), х — 0, что эквивалентно соотношению р(х) — (1 — х + х2/2 — х3/6 + х4/24) = 0(х5). Подставляя р(х), нетрудно получить
2 3 4
/•v* ry '-> /V»"*
, ч K\J K\J K\J
р(х) - 1 + х - Т + ¥ - 24
(2 - х)Л
x x2 x3 24(1 + Лх2) + V - 2+"6 - 24 J £
х5
(1 + £х4)
= 0(х5).
Таким образом, четвертый порядок имеет место при любых значениях параметров £, Л. Абсолютная устойчивость схемы заключается в выполнении условий —1 < р(х) < 1 при х > 0. Условие р(х) < 1 будет выполнено, если выражение в квадратных скобках в соотношении (5) неотрицательно, что эквивалентно неравенству
(1 + Лх2) ( 1 - - + уз ) +
(2 - х)х2
24
0.
В результате тождественных преобразований оно может быть приведено к виду
1 - х+ х2
1 , / Л2 Л
6+ Л - 3+4Л 12
+ х 3
2 Лх2 4Л / 12
0.
2
2
Если положить Л = 1/5, то дискриминант квадратного трехчлена, образованного первыми тремя слагаемыми в данном неравенстве, отрицателен, а значит, этот трехчлен положителен, и неравенство выполнено. Естественно, значение Л =1/5 носит характер достаточного и может быть завышенным. Как показали расчеты на ЭВМ левой части неравенства (4), минимальное значение параметра Л, гарантирующее ее неотрицательность, равно 0.19122.
Для получения абсолютно устойчивой схемы необходимо также обеспечить выполнение условия —1 < р(х) при х > 0, которое эквивалентно неравенству
/ х2 х3 \ (2 _ х)х3
(х — 1 + У (1 + Лх2) + р < 2(1 + £х4)(1 + Лх2).
Подбор на ЭВМ минимального значения параметра £, гарантирующего выполнение этого неравенства, дает £ = 1.4130 ■ 10-3 для Л =1/5 и £ = 1.3932 ■ 10-3 для Л = 0.19122.
Данная схема имеет четвертый порядок аппроксимации при любых значениях параметров £, Л. Поэтому можно использовать округленные значения параметров, например £ = 1.5 ■ 10-3 и Л = 0.2, которые, как и точные значения, гарантируют абсолютную устойчивость схемы.
2. Схемы для уравнения переноса
В данном разделе строятся схемы, которые могут рассматриваться как обобщение схемы (1)-(4) для решения обыкновенных уравнений. При численном исследовании разрывных решений необходимо введение диффузии в какой-либо форме. Однако предварительный анализ показал, что центрирование по времени диффузионных членов, идентичное центрированию конвективных членов, не порождает абсолютно устойчивую схему. Построение схемы для уравнения с диффузией оказалось существенно более сложной задачей. Поэтому пока мы ограничимся рассмотрением m-мерного уравнения переноса без диффузии
<л m г\
du du
-r—+ / aiо— = 0, ai = const, i = 1,..., m, 1 < m < 3. (6)
at ' dxl i=i l
Введем разностный оператор конвекции
Ru = £ aiA^ 1 - ^ ) u = £ aiЦ + O(h4), (7)
h = maxm= 1 (Axi), и вспомогательные операторы, используемые при построении схемы:
mm
R = £afA+Af, R = £a4(A+Af)2. l= l=
Здесь и далее приняты обозначения 8± u = ±[u(x 1,...,xi ± Axi, ...,xm) — u(x 1,...,xm)],
A± u = 8± u/Axi. Рассмотрим четырехшаговую схему
u(1) _ un
U U + R un = 0, (8)
т/2
?/2> - w(1)
(1 - AT2R2)-+ Ri(u(1) - = 0, (9)
т/2
И(3) -
t
+ Riu(2) = 0, (10)
(1 + 5т4 ^ + д,( Т +- + - + -)=00 (11)
В этих уравнениях Л — параметры, которые определим далее. Исключим промежуточные функции . В результате можно получить
„ п+1_„п
(1 + т4Д4^(1 - ЛТ2Я2)--—+
т
+ (1 - Лт2Д2^Д1 - 2д2 + 12 -п + ^д3 - ^^ -п = 0. (12)
Для исследования порядка по времени умножим это уравнение на оператор (1-Лт2Д2)-1 и перегруппируем слагаемые:
-п+1 - -п / „ т „о т2 о т3
+ ( - 2 + у- 24) +
/ т 2 т 3\ —п+1 _ —п
+ [(1 - Лт2Я2)-1 - 1] (Д312 - д4-п + т4Д4^- т - = 0.
Первые два слагаемые этого уравнения соответствуют схеме четвертого порядка по времени, так как второе слагаемое содержит разложение по времени в ряд Тэйлора с нужным количеством членов. Третье и четвертое слагаемые этого уравнения могут быть оценены как 0(т4), поэтому сконструированная схема имеет четвертый порядок по времени.
Для исследования устойчивости применим метод Фурье. Представив решение в виде
-п = рп ехр ^^ щц^ , 0 < а < 2п, %1 = 1, ... , I, ] = V-1, получаем множитель перехода
для схемы (12) р =1 - ^'Р^/А - С/А, где
р 2 р 4 / р 2 \ р 2
С = - —; А = X(1+ №); В = £ 1 - —--^; X =1 + ЛР2;
2 24 V 12/12
Pi = ^ аг sin а N +
т
г=1
6 / Ахг
P2
Е
l=1
агт Ахг
Ф;
P4
Е
агт Ахг
(13)
l=1
т. е. P1 и P4 являются спектральными образами операторов тР1 и t4R4, P2 есть спектральный образ оператора — т2R2, q^ = 2-2 cos . Условие устойчивости | р |< 1 для этой схемы
записывается в виде (1 — C/A)2 + (P1B/A)2 < 1 или P2B2 < C(2A — C), т.е.
P12
P2 P2
Z (1 — ^) — ^ 1 12 ; 12
<
P 2 P 4
Z—1--P1-
2 24
P 2 P 4
2(1 + £P4)Z — Z^ + _4
(14)
2
4
2
Отметим соотношения между величинами Р\, Р2, Р4, которые потребуются при исследовании устойчивости [9, 10]:
Р2 < вР2, в = т(1 + 7/243), Р22 < вР4- (15)
Теорема. Схема (8)-(11) устойчива при любом шаге по времени т, если
Л = в/6, £ = в 3/72.
Действительно, подставив эти выражения в неравенство (14), заменив Р4 на его минимальное ввиду второго из неравенств (15) значение Р22/в и перегруппировав слагаемые, можно получить
0 < Р?
72 Р2 (вР2 - р2) +6 Р2 - Ш+144(2+4 Р2 -т) +
+ ^Р^ в2Р| - )+^Р22С в2Р| - вР2Р2 - ^
432 V 2 4 у 2592 2 V 2 2 1 2
Выражения в каждой из круглых скобок являются невозрастающими функциями аргумента Р^ > 0. Следовательно, эти выражения принимают минимальное значение при максимально возможном значении аргумента Р 2, равном в силу первого из неравенств (15) вР2. Прямой подстановкой можно убедиться, что минимальные значения каждого из этих выражений неотрицательны, и нетрудно видеть, что правая часть данного неравенства поэтому неотрицательна, т. е. оно верно. Теорема доказана. В
Таким образом, соотношения Л = в/6, £ = в3/72 обеспечивают устойчивость схемы (8)-(11). В одномерном случае эти схемные параметры принимают значения Л = с/6, £ = с3/72, в двумерном случае — Л = с/3, ^ = с3/9, в трехмерном — Л = с/2, £ = с33/8, с = 250/243. Следует отметить, эти значения обеспечивают достаточные, но не необходимые условия устойчивости. Возможно, существуют и другие наборы параметров исходной схемы, удовлетворяющие условию устойчивости в приближении метода Фурье.
Сделаем некоторые замечания по данной схеме.
1. Числовой коэффициент 250/243 определяется использованием формул (7) четвертого порядка аппроксимации для конвективных слагаемых уравнения (6). Возможны другие варианты аппроксимации оператора конвективного переноса (естественно, лишь антисимметричные, имеющие чисто мнимый спектральный образ). В этом случае приведенное выше доказательство устойчивости остается справедливым, лишь изменится данный коэффициент.
2. В многомерном случае реализация схемы достаточно проста лишь при использовании факторизации стабилизирующих операторов. Однако для факторизованной схемы приведенное выше доказательство устойчивости теряет справедливость.
3. Факторизация стабилизирующих операторов не портит аппроксимационных свойств схемы как по времени, так и по пространству.
4. Схема имеет свойство полной аппроксимации, т. е. стационарное решение не зависит от временного шага, что позволяет использовать схему в рамках метода установления и выбирать оптимальный шаг по времени в соответствии с критерием наиболее высокой скорости сходимости.
3. Схемы для уравнения переноса с диффузией
Необходимым условием конкурентоспособности схем повышенной точности для интегрирования гиперболических уравнений первого порядка является способность качественного моделирования разрывных решений. Для подавления появляющихся при этом осцилля-ций необходимо вводить искусственные диффузионные слагаемые. Однако оказалось, что естественные наиболее устойчивые, казалось бы, способы центрирования их по времени не приводят к успеху. Анализ устойчивости различных вариантов аппроксимации диффузионных членов позволил построить нижеприведенную схему. Итак, рассмотрим т-мерное уравнение
О т ~ т
ОП тг-^ ОП тг-^ О2 и
дь
дх1 ^ ' дх2'
а1 = сопэ^ = сопэ^ I = 1,
1=1 ' 1=1 Введем разностный диффузионный оператор
, т.
т т п2
Уи = £ *АГА+и = £ * ^ + 0(П2). 1=1 1=1 1
Рассмотрим четырехшаговую схему
и(1) - ип
т/2
+ Ехип = Уип
(16)
(1 - Ат 2К2)и (2)т/2и(1) - (2 У - т ^^ + Я1(и(1) - ип) = 0,
и(2) - ип
т/2
т/2
(17)
и
(3)- ип
+ ^и(2) = Уип
(18)
(1 + £т4Я4 - птУ + т37УР2)
ип+1 - ип
+ 2'ПУ (и(1) - ип) +
и"
+*1{ -6+ 3
и(1) и(2)
+
3
+
и
(3)
6
Уи(2).
(19)
Здесь А — параметры, заданные в предыдущем разделе, 7 и ц определим ниже. Отметим, что при любых значениях этих трех параметров диффузионные члены аппроксимируются со вторым порядком. Для исследования устойчивости применим метод Фурье.
Представив решение в виде ип = рпехр ^ ^ а1 , 0 < а1 < 2п, ц =1,..., 11, ] =
получаем множитель перехода для схемы (16) - (19) р =1 - ]Р1Б/А - С/А, где
W
3
X W\
12 24 )
1
W
Б = X 1--- Р2---- ----7 (^Х + Р2), X = X + W/2 + рШР2
12 24
Р2
3
fWZ Р4 24
С = X(W + ^ ] - Р2 ( ^ + ) + Р2 ( — - ) (2Х - W);
12
1 Ж
12 24
А = X(1+ ^ + П^ + 7^2), Ж = ^ ^гдг/Аж?, ® = 2 - 2со8(а)
1=1
величины Р-, Р2, Р4, ! определены формулами (13).
Условие устойчивости |р| < 1 для этой схемы записывается в виде (1 — С/А)2 + (Р1В/А)2 < 1, или Р2В2 < С(2А — С), т.е.
Р2
т)
Xi1 — — р2 ! — Ж + 1
12
+(т—о т
<
<
2Х(1 + СР4 + пТ + 7^Р2) —X (V + Р^ + Р2
Р,2 --+ —
12 24
— Р2
1—Т)
12;
X
X
х + Р2) — Ч !!
+
Р2
24
+(1—I) (!
12
(20)
При Ж = 0 это неравенство превращается в уже доказанное неравенство (14). Можно несколько упростить анализ устойчивости, заменив Р4 на меньшее ввиду второго из неравенств (15) выражение Р22/в. Рассмотрев асимптотический случай большой диффузии Ш ^ ж, можно показать, что существуют по меньшей мере следующие два набора свободных параметров, гарантирующие выполнение неравенства (20) в этом случае, а именно ^ = в/24, 7 = в/36, п = 1/3 и ^ = в/72, 7 = в/18, П = 1/2. Однако проверка на ЭВМ этого неравенства во всей области изменения переменных Рр > 0, Р2 > Рр/в, ! > 0 показала, что для каждого из приведенных выше наборов параметров существуют значения переменных Р^2, Р2, !, при которых неравенство не выполняется, т.е. устойчивость схемы не обеспечивается. Далее с помощью ЭВМ был проведен поиск значений параметров ^ и 7, гарантирующих выполнение неравенства (20). В результате перебора нескольких десятков вариантов вместо указанных выше двух наборов получены следующие два, обеспечивающие абсолютную устойчивость схемы: ^ = в/12, 7 = в/22, п =1/3 и ^ = в/15, 7 = в/18, П = 1/2. По-видимому, представляет интерес задача минимизации значений параметров стабилизирующих операторов с определением единственного оптимального их набора. Ее решение здесь не производилось.
2
4. Численное решение многомерного уравнения переноса
Итак, абсолютно устойчивые схемы для линейных уравнений построены. Однако при численном решении как линейных, так и нелинейных задач возникают дополнительные проблемы. Основной из них является необходимость подавления осцилляций вблизи разрывов решения. В принципе для этой цели разработаны эффективные алгоритмы [6], но они обоснованы лишь для явных схем или неявных схем с простейшими стабилизирующими операторами. Здесь для подавления осцилляций будет использоваться схемная диффузия в виде параметрических комбинаций слагаемых низкого и высокого порядков [11, 12]. Параметр, переключающий режимы диффузии, будет вычисляться локально по соотношению градиентов решения.
Рассмотрим двумерное линейное уравнение переноса
+ 5 (а/)/дж + 5 (ау / )/ду = 0, ах = ах(ж,у), ау = ау (ж, у).
Для его решения используется аналог схемы (16)-(20). Однако реализация схемы более удобна при максимальной однородности шагов схемы. Этой однородности у схемы (16) - (20) нет, так как второй шаг соответствует разности двух уравнений. Поэтому в нижеследующей схеме произведена обратная операция, т. е. второе уравнение соответствует сумме уравнений (16) и (17).
Используя нелинейную диффузию в консервативной форме, запишем
7,(!) _ ип
- + К^ = V пи(2), г/2
т т ?/(2) _ ип
(1 - Ат2 Я2х - - к1 + 2т 3^ТКХ)(1 - Ат 2Д2у - - Уу + 2т3^ТКУ) т/2 +
+2ЛтДп(м(1) - ип) + К^(1) = У"ига,
и(3) - иП + Д1«(2) = т
(1+ет4 К2 - 2 к(2) + т з27к(2)д2Х))(1+^2) - 2 ^(2) + т^2^) +
(2)(и(1) - ип) + кА/П + ^ + ^ + ^ = V(2)и(2).
\ 6 3 3 6 /
Здесь е = в3/72, А = в/6, ^ = в/15, 7 = в/18 — параметры, заданные в предыдущих разделах; п положен равным 1/2; К1 — разностный оператор конвекции, отличающийся от оператора (7) лишь внесением переменных коэффициентов с целью сохранения консервативности под оператор разностного дифференцирования
К^а = Д£(1 - ¿-¿+/6)а^ий + Ду(1 - /6)ауиц;
Vк — квазилинейный оператор, связанный с диффузионной правой частью уравнений соотношением V(к)и(к) = Д-^+1/2 + Д-^(+у/2; ^(к) — искусственные адаптивные диффузионные потоки, удовлетворяющие трем условиям [11]. Во-первых, они имеют высокий порядок (пятый в данном случае) на "гладких" решениях и первый возле разрывов; во-вторых, они порождают неотрицательные стоки разностной "энергии"; в-третьих, они могут быть тождественно преобразованы к двухточечному выражению с неотрицательными коэффициентами. В [11] получены формулы
^+1/2 = 2(2 - - + 32и»',
где 0 < дц < 1 — индикаторы локальной гладкости — вычисляются следующим образом:
ди = тах[0, тт(аз+, 5-, 1) - с]/(1 - с), = [(<Ти«)(5+иа) + еР)/[(Гигг)2 + 10-9], = [(¿+иЙ )(ГиЙ) + еи]/[(<ТиЙ )2 + 10-9 ];
константы а и с позволяют регулировать свойства диффузии; параметр ер предназначен для предотвращения понижения порядка диффузии на гладких экстремумах или участках почти постоянного решения. По результатам пробных расчетов были выбраны следующие значения эмпирических величин: с = 0.4, а = 1.3, еи = 2(итах - ит1П).
Сечение x = 0: сетка 201 х 201 (а), сетка 401 х 401 (б).
С помощью описанной схемы численно решалось уравнение переноса с коэффициентами ах = —у, ау = х, которые соответствуют вращению произвольного начального распределения вокруг точки х = 0, у = 0 с единичной угловой скоростью. Расчетная область является квадратом —1/2 < х < 1/2, —1/2 < у < 1/2 и содержит 201 х 201 и 401 х 401 узлов. Начальные данные представляют цилиндр единичной высоты радиуса 1/8 с центром в точке х = 0, у = 1/4 на нулевом фоне. На рисунке приведены расчетные результаты для сечения х = 0. Точками изображено начальное распределение, прямоугольниками — результаты расчета на момент времени t = п, когда исходная картина должна развернуться на 180° относительно центра расчетной области. На рис. а приводятся результаты для сетки 201 х 201, сделано 100 шагов по времени, на рис. б — для сетки 401 х 401, сделано 200 временных шагов, рисуется каждый второй узел. В обоих случаях максимальное локальное число Куранта равно 3.16.
Итак, представленная в работе схема четвертого порядка действительно позволяет решать уравнение переноса при числах Куранта, больших единицы, она сохраняет устойчивость и аппроксимационные свойства при факторизации и поэтому может быть эффективной при решении многомерных задач, гиперболических и параболических.
Список литературы
[1] Русанов В. В. Разностные схемы третьего порядка точности для сквозного расчета разрывных решений // Докл. АН СССР. 1968. Т. 180, №6. С. 1303-1305.
[2] BuRsTEIN S. Z., MlRIN A. A. Third order difference methods for hyperbolic equations // J. Comp. Phys. 1970. Vol. 5, No. 3. P. 547-571.
[3] Kutler P., Lomax H., Warming R. Computation of space shuttle flow fields using noncentered finite-difference schemes // AIAA paper. 1972. No. 72-193. 25 p.
[4] Shu C.-W. Total-variation-diminishing time discretizations // SIAM J. Sci. Stat. Comp. 1988. Vol. 9, No. 6. P. 1073-1084.
[5] Osher S., Shu C.-W. Efficient implementation of essentially nonoscillatory shock-capturring schemes // J. Comp. Phys. 1988. Vol. 77, No. 2. P. 439-471.
[6] Пинчуков В. И., Ч.-В. Шу. Численные методы высоких порядков для задач аэрогидродинамики. Новосибирск: Изд-во Сиб. отд-ния РАН, 2000. 232 с.
[7] Толстых А. И. Компактные разностные схемы и их приложения к проблемам аэрогидродинамики. М.: Наука, 1990. 230 с.
[8] Пинчуков В. И. Компактная схема шестого порядка для решения уравнений Эйлера в криволинейных координатах // Журн. вычисл. математики и мат. физики. 1998. Т. 38, №10. С. 1717-1720.
[9] Пинчуков В. И. О неявных схемах типа Рунге — Кутты третьего порядка по времени // Вычисл. технологии. 1999. T. 4, №2. C. 59-73.
[10] Пинчуков В. И. Абсолютно устойчивые схемы Рунге — Кутты третьего порядка // ЖВМ и МФ. 1999. Т. 38, №11. С. 1855-1868.
[11] Пинчуков В. И. О численном исследовании трансзвуковых турбулентных течений возле крыла неявными схемами высоких порядков // Вычисл. технологии. 2001. T. 6, №2. C. 110-121.
[12] Пинчуков В. И. Нелинейные сеточные фильтры и их использование в схемах высоких порядков для задач аэродинамики // Мат. моделирование. 1999. Т. 10, №11. С. 111-115.
Поступила в редакцию 11 июля 2001 г.
