CC BY
55
Вычислительные технологии
Том 4, № 2, 1999
O НЕЯВНЫХ СХЕМАХ ТИПА РУНГЕ-КУТТЫ ТРЕТЬЕГО ПОРЯДКА АППРОКСИМАЦИИ *
В. И. Пинчуков Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: Pinchvi@net.ict.nsc.ru
Implicit schemes of the third-order approximation with respect to time are suggested. Ordinary equations of the first order are considered as well as the multidimensional convection equations. The schemes are proved to be absolutely stable in the frozen coefficients approximation. The results of test calculations are presented.
В вычислительной практике широко распространены схемы первого и второго порядков аппроксимации по времени. Схемы более высокой точности по временной переменной также известны, хотя и не столь популярны. Это явные схемы третьего порядка типа Рунге — Кутты — Русанова [1], Бернстейна — Мирина [2], Кутлера — Ломекса [3]. Вследствие жесткого ограничения временного шага по условиям устойчивости они, как и конкурирующие с ними явные схемы второго порядка, могут использоваться лишь с небольшим временным шагом. При этом временная составляющая погрешности обоих типов схем невелика, и преимущество по точности схем третьего порядка не имеет большого значения, в то же время они значительно сложнее и дороже по количеству арифметических операций.
Поэтому представляет интерес построение неявных схем третьего порядка, которые могут использоваться при большом временном шаге. Временная погрешность при этом существенна, поэтому они могут иметь преимущество по сравнению со схемами первого — второго порядков. Создание неявных схем повышенной точности связано с преодолением определенных технических сложностей. Чтобы пояснить суть проблемы, укажем, что в классических неявных схемах стабилизирующий оператор при воздействии на временное приращение искомых переменных образует величину порядка 0(т), и, как правило, схема имеет такой же порядок 0(т). Если стабилизирующие слагаемые являются приближением второй производной по времени от искомых переменных, Rstab(fn+l — fn) l(df/dt)n+ — (df/dt)n]/2 « (T/2)d2f/dt2, то порядок аппроксимации по времени можно повысить на единицу и добиться второго порядка. При построении схем третьего порядка подобные подходы не работают.
Если при создании неявных схем повышенной точности отталкиваться от каких-либо явных схем, например, типа Рунге — Кутты, то приходится конструировать стабилизирующие операторы, с одной стороны гарантирующие устойчивость, с другой — удовлетворяющие ограничениям по точности, которые являются более жесткими, чем для схем
* Работа выполнена при финансовой поддержке Американского фонда гражданских исследований и развития для независимых государств, образованных на территории бывшего Советского Союза (CRDF), грант RM1-212, и Российского фонда фундаментальных исследований, грант №97-01-00711.
© В. И. Пинчуков, 1999.
второго порядка. Таким образом, построение неявных схем третьего порядка сопряжено с преодолением трудностей. В частности, в [4] построены компактные неявные схемы третьего порядка, однако эти схемы являются условно устойчивыми. В работе приведен пример построения неявных абсолютно устойчивых схем.
Упомянутые выше явные схемы Кутлера — Ломекса и Русанова могут рассматриваться как обобщение известных схем Рунге — Кутты третьего порядка для обыкновенных дифференциальных уравнений [5]. Предлагаемые схемы допускают также построение аналога для обыкновенных уравнений, с чего и будет начато их изложение. Отметим, что цель создания неявных методов может трактоваться как расширение их области устойчивости. Определенные возможности для этого имеются и в классе явных методов, примером чему могут служить явные методы с заданной областью устойчивости [5].
1. Двухшаговые схемы
Пусть решается задача Коши для уравнения
du/dt = <^(u,t), u(0) = (1)
Рассмотрим схему
[3/2 — т^П/2 + А(т V,I)2](«n+2/3 — un) = T?(un, t), (2)
[1 + С(тVS2](«"+1 - un) - Ст2(тVn)2(un+2/3 — un)3/2 =
= т [^(un,t)/4 + ^(un+2/3,t + т 2/3)3/4], (3)
где £ > 0, А > 0 — числовые параметры, которые следует определить в анализе. Оказалось, что необходимым и достаточным условиям точности и устойчивости схемы можно удовлетворить, если положить
С = 2А. (4)
Чтобы не загромождать изложения, все выкладки далее проведем с учетом этого соотношения.
Для исследования устойчивости и точности этой схемы рассмотрим случай <^(u,t) =
—Au, А = const > 0. Умножая уравнение (2) на Ат2^Ц3 + т^Ц3/4 = Ат2А23 — тА3/4, уравнение (3) на 3/2 — т^Ц/2 + Ат2= 3/2 + тА/2 + Ат2A2 и складывая их, получаем соотношение un+1 = ипр(тА), где
р(х) = 1 — х[1 — х/6 + х2А8/3]/[1 + х/3 + Ах22/3]/[1 + 2Ах2]. (5)
Абсолютная устойчивость схемы заключается в выполнении условия | р(х) |< 1 при
х > 0, что эквивалентно двум неравенствам:
р(х) < 1, р(х) > —1, которые с учетом ограничения x > 0 переписываются в виде зависимостей
1 — х/6 + х2А8/3 > 0,
х(1 — х/6 + х2А8/3) < 2(1 + х/3 + Ах22/3)(1 + 2Ах2), при этом первое эквивалентно ограничению А > 1/384, второе — неравенству
2 — х/3 + (А16/3 + 1/6)х2 — (А4/3)х3 + А2х48/3 > 0,
которое можно преобразовать к виду
1 + (1 - х/6)2 + х2{Л16/3 + 5/36 - (Л4/3)х + х2Л28/3} > 0.
Условие неотрицательности квадратного трехчлена в фигурной скобке эквивалентно следующему соотношению на его коэффициенты:
(Л16/3 + 5/36)А28/3 > (А2/3)2,
которое выполняется при Л > 1/192. Итак, можно положить Л = 1/192 (условие Л > 1/384 при этом также выполняется), £ = 1/96, при этом будем иметь абсолютно устойчивую схему.
Чтобы иметь третий порядок, необходимо выполнение условия р(х) — ехр(—х) = 0(х4), х ^ 0, что эквивалентно соотношению р(х) — (1 — х + х2/2 — х3/6) = 0(х4). Подставляя р(х), нетрудно получить
р(х) — (1 — х + х2/2 — х3/6) = х4[Л2х32/9 + (—Л2/3 + 1/9)Лх2+ + (Л/9 + Л24/3)х — 2Л/3 + 1 /18]/(1 + х/3 + Лх22/3)/(1 + 2Лх2) = 0(х4).
Таким образом, данная схема имеет третий порядок аппроксимации, причем при любых значениях параметров £, Л. Это означает, что можно использовать упрощенные значения величины <рЦ. В частности, при решении систем обыкновенных уравнений на последнем шаге Якобиан можно заменить на мажорирующую диагональную матрицу, обращение которой элементарно.
Аналогичное временное центрирование правой части и стабилизирующих слагаемых, совместимое с нужным нам третьим порядком по времени, можно использовать при построении схемы для т-мерного уравнения переноса
т
д//д* + ^ агд//5хг = 0. (6)
1=1
Введем разностный оператор конвективного переноса
т т т
#1/ = X! а1д°(1 — /6)/ = X! а«д//дх1 + Дх4) (7)
1=1 1=1 1=1
и вспомогательные операторы, используемые при построении схемы,
т
#2 = X]а12А+Д-, Я4 = Е1=1 а4(А+Д-)2. (8)
1=1
Здесь и далее £±/ = ±(/(х1,..., х1 ±Ах1,..., хт) —/(х1,..., хт)), Д±/ = ^±//Дх1. Рассмотрим двухшаговую схему
(1 + т Я1/3 + £г4Д4 — Лт2 Д2)(/* — /п)3/(2т) + Я/ = — т3 р^Г, (9)
(1 + ^т4Я4 — пт2 Я2)(Г+1 — / п)/т + птЯ2 (/* — /п)3/2+
+ВД п/4 + / *3/4) = -т 3рД4/п,
(10)
где £, Л, п — параметры, которые необходимо определить при анализе устойчивости
схемы. Чтобы не загромождать изложения, воспользуемся некоторыми итогами такого анализа и ограничимся рассмотрением двухпараметрического семейства Л = п, £ = <Р = ^/2. Для анализа аппроксимационных свойств схемы умножим уравнение (9) на оператор —Лт2Я2 — Rit/2, (10) — на оператор 1 + rRi/3 + r4R4^/2 — Лт2R2 и сложим их, исключив тем самым промежуточную функцию f *. В результате получим уравнение
(1 + т 4Я4^/2 — Лт 2R2)(1 + TRi/3 + т 4R4^/2 — Лт 2R2)(fra+1 — Л/т+
+ [1 — т Ri / 6 + т4^4^/2 — 2Лт 2R2](Ri — т VRf = 0, которое может быть преобразовано к следующему виду:
fn+1 = fn — (1 + т 4R4^/2 — Лт2 R2)-1 (1 + т#1 /3 + т 4R4^/2 — Лт 2R2)-1x
X [1 — т#1/6 + т4R4^/2 — 2Лт2R2](R1 — т3^4)т^.
Разлагая правую часть в ряд по малому параметру т, можно получить уравнение
fn+1 = fn — (1 — тRl/2 + т ^^Rf + 0(т4).
Учитывая, что R1f= — df /dt+O(h4), h = тах™1(Аж1) и соответственно R^f = (—1)kdkf/dtk+ O(h4), имеем
fn+1 = (f + тд^д* + т2д2f/dt2/2 + т3d3f/dt3/6) Ur +0(т4 + тЛ4).
Поскольку это выражение совпадает с разложением в ряд Тэйлора до слагаемых третьего порядка, из него следует, что схема имеет третий порядок аппроксимации по времени. Она имеет четвертый порядок по пространственным переменным, так как разностный оператор конвективного переноса (7) имеет четвертый порядок.
Для исследования устойчивости применим метод Фурье. Представив решение в виде fn = pn exp(Em=1 ), 0 < а < 2n, j = \/ —1, получим множитель перехода для схемы
р = 1 — (D + ЛР2 — jP1/6)(jP1 + ^P4)/(D + ,7Р1/3)/Д
где
m m
P1 = ^2 a sin ai(1 + ®/6)т/Ажг, P2 = 5^(агт/Ажг)2®,
1=1 1=1
m
P4 = ^(агт/Ажг )4qf, qi = 2 — 2cos D =1 + P4^/2 + Р2Л. (11)
1=1
Отметим соотношения между величинами P1, P2, P4, которые потребуются при исследовании устойчивости:
P2 < вРг, в = m(1 + 7/243), P22 < mP4. (12)
Докажем первое из неравенств (12) в одномерном случае. Перепишем его, опуская индекс /, поскольку он принимает единственное значение, равное единице:
(a sin а(1 + д/6)т/Аж)2 < в(ат/Аж)2д.
Подставляя (sin а)2 = q(1 — q/4), сокращая одинаковые сомножители слева и справа, получим выражение
(1 — q/4)(1 + q/6)2 < в.
Поскольку функция в левой части этого неравенства имеет экстремум в точке q = 2/3, где она принимет максимальное значение на интервале изменения переменной q : 0 < q < 4, то, подставляя q = 2/3, получаем (1 —1/6)(1 + 1/9)2 = 250/243 < в. Таким образом, при m = 1 неравенство (12) доказано. При m> 1 квадрат суммы в левой части неравенства (12) порождает смешанные произведения слагаемых, оценивая которые по формуле 2ab < a2 + b2, сводим нужное нам неравенство к сумме одномерных неравенств с одним отличием — слева оценка смешанных слагаемых породила множитель m, который учтен аналогичным множителем в определении параметра в.
Аналогично, второе из неравенств (12) достаточно рассмотреть в одномерном случае; параметр m учитывает отличие многомерного варианта от одномерного. Опуская в одномерном случае индекс /, имеем
[^т/Аж)^]2 < ^т/Аж)^2.
Ясно, что левая и правая часть этого неравенства тождественно равны. Таким образом, неравенства (12) доказаны.
Исследование устойчивости рассматриваемой схемы сводится к доказательству условия | р |2< 1. Оно будет выполнено, если после приведения к общему знаменателю слагаемых коэффициента перехода квадрат модуля знаменателя будет больше квадрата модуля числителя, т. е.
p2D2/9 + D4 > [D2 — (D + ЛРг)^ — P2/6]2 + p2(D2/3 + ЛР2 — ^Pi/6)2.
Если перенести все слагаемые этого неравенства в левую часть, то получим квадратный трехчлен относительно величины Pj2. Поскольку коэффициент при второй степени аргумента отрицателен, то минимальное значение трехчлена достигается на границах отрезка определения аргумента, т. е. ввиду его неотрицательности и с учетом первого из соотношений (12) — в точках Pj2 = 0 и Pj2 = вР2. Таким образом, надо проверить два неравенства:
D4 — [D2 — (D + ЛР2)^Р4]2 > 0, D4 — [D2 — (D + ЛР2)^Р4 — вР2/6]2+ +вP2[D2/9 — (D2/3 + ЛР2 — ^Р4/6)2] > 0.
Учитывая, что D4 — [D2 — (D + ЛР2)^Р4]2 = [(1 + Р2Л)(2 + 2Р2Л) + ^P4](D + ЛР2)^Р4, первое неравенство выполнено всегда. Аналогично, подставляя определение D и раскрывая разности квадратов, второе неравенство перепишем в следующем виде:
2[^Р4(1 + ^Р4/2 + 2Р2Л) + вР2/6][1 + ^Р4/2 + 2ЛР2 + (ЛР2)2 — вР2/12]-
— (вР2/3)(1 + 2ЛР2 + ^Р4/3)(1 + 4Р2Л) > 0.
Собирая члены, содержащие одинаковые произведения переменных Р2, Р4, можно получить неравенство
{Р43^3/2 + Р|Р|(^Л)2} + {р42р2(4^2Л — в^2/12) + Р4Р23^Л34} +
+{Р422^2 + Р4Р22 (10^Л2 — в^Л7/9)} + {Р4^2(^Л8 — в^/9) — Р23вЛ27/3}+
+{Р4(^2) — Р*2 (вЛ4/3 + в2/36)} > 0. (13)
В левой его части сгруппированы в фигурные скобки слагаемые, которые, ввиду неравенств (12), удобно рассматривать совместно. Отметим: во втором из неравенств (12) вместо параметра т будем использовать параметр в, что лишь усиливает это неравенство и позволяет упростить дальнейшие выкладки. Заменим в третьей фигурной скобке неравенства (13) Р42 на Р4Р22/в, тем самым уменьшив, возможно, в силу неравенств (12), значение выражения в скобке. Потребовав, чтобы получившееся выражение было неотрицательным при любых неотрицательных значениях Р2, Р4, получим неравенство 2^2/в + 10^Л2 — в^Л7/9 > 0. В следующей скобке заменим Р23 на вР2Р4, также уменьшив при этом, в силу (12), значение выражения. Аналогичным образом сделаем однородными выражения в остальных скобках. В результате получим систему неравенств для коэффициентов схемы:
^3/(2в) + (М)2 > 0, 4^2Л — в^2/12 + в^Л34 > 0,
2^2/в + 10^Л2 — в^Л7/9 > 0, ^Л8 — в^/9 — в2Л27/3 > 0,
^2/в — вЛ4/3 — в2/36 > 0.
Нетрудно проверить, что данным неравенствам удовлетворяют параметры Л = в/48, ^ = в3/36. Таким образом, в одномерном случае это параметры Л = с/48, ^ = с3/36, в двумерном — Л = с/24, ^ = с32/9, в трехмерном — Л = с/16, ^ = с33/4, с = 250/243. Следует отметить, что наложенные нами условия на коэффициенты полинома (13) имеют характер достаточных, но не необходимых для его положительности. Возможно, существуют и другие наборы параметров исходной схемы, удовлетворяющие условию устойчивости в приближении метода Фурье.
2. Трехшаговые схемы
Выше были рассмотрены неявные схемы третьего порядка аппроксимации по времени для решения обыкновенных уравнений или уравнения переноса. Далее будут построены другие схемы третьего порядка, по некоторым качествам обладающие преимуществами по сравнению с приведенными выше.
Как и ранее, начнем со случая обыкновенных уравнений.
Итак, рассмотрим снова задачу Коши для уравнения (1). Рассмотрим трехшаговую схему
(и* — ига)3/2 = т^(ига,£), (14)
(и** — ип)3/2 = т [<^(м*, £ + т 2/3) + ^(ига, £)]/2,
(15)
(1 + £т4^)(иП+1 — «га) = т HuV)/4 + <^(u*V + т 2/3)3/4], (16)
где £ > 0 — числовой параметр, который следует определить в анализе устойчивости.
Для исследования устойчивости и точности этой схемы рассмотрим случай <^(u,t) = —Au, А = const > 0. Введя параметр x = тА, умножая уравнение (14) на выражение —тА/3 = — ж/3 и складывая с (15), имеем
(u** — un) = —xun2/3 + 2un(x/3)2.
Уравнение (16) в результате подстановки <^(u,t) = — Au принимает вид
(1 + £x4)(un+1 — un) = —x(u**3/4 + un/4).
Умножая предпоследнее уравнение на выражение —x3/4 и складывая с последним, получаем
(1 + £x4)(un+1 — un) = (—x + x2/2 — x3/6)un
Таким образом, коэффициент перехода для схемы имеет вид
p(x) = 1 — x(1 — x/2 + x2/6)/(1 + £x4). (17)
Абсолютная устойчивость схемы заключается в выполнении условия | p(x) |< 1 при x > 0, что эквивалентно двум неравенствам:
p(x) < 1, p(x) > —1, которые с учетом ограничения x > 0 переписываются в виде соотношений
1 — x/2 + x2/6 > 0,
x(1 — x/2 + x2/6) < 2(1 + £x4).
Первое из них выполнено всегда, второе эквивалентно неравенству
0 < (1 — x/4)2 + x2(3/4 — x/6)2/3 + (£ — 1/108)x4,
которое выполнено, если £ > 1/108. Итак, положив £ = 1/108, будем иметь абсолютно устойчивую схему.
Чтобы иметь третий порядок, необходимо выполнение условия p(x) — exp(—x) = O(x4), x ^ 0, что эквивалентно соотношению p(x) — (1 — x + x2/2 — x3/6) = O(x4). Учитывая, что в (17) в числителе выражение совпадает с разложением экспоненты до нужной степени, а в знаменателе первое отличное от единицы слагаемое имеет порядок (x4), аппроксимаци-онное условие выполнено, и схема имеет третий порядок по времени.
Отметим, что данная схема имеет третий порядок аппроксимации при любых значениях параметра £, причем можно использовать упрощенные значения величины <^. Это позволяет облегчить реализацию схемы при решении систем большого количества обыкновенных уравнений.
Аналогичная схема со стабилизирующим оператором, содержащим четные степени шага по времени, может быть сконструирована для m-мерного уравнения переноса (6).
Рассмотрим трехшаговую схему
(f * — fn)3/(2т) + R1 fn = 0, (18)
(f ** — fn)3/(2т) + R(fn + f *)/2 = 0, (19)
(1 + пт^4)(Г+1 — Г)/т + R1(fn/4 + f **3/4) = 0. (20)
Здесь R1, R4 — разностные операторы, определенные формулами (7), (8), п — параметр, который определим позже.
Исключим разностную функцию промежуточного слоя f *, умножив уравнение (18) на оператор — R^/3 и сложив его с уравнением (19). Получим уравнение
(f ** — Л3/(2т) + Rfn — R2fn/3 = 0.
Исключим разностную функцию промежуточного слоя f **, умножив последнее уравнение на оператор — R^/2 и сложив его с уравнением (20), получим уравнение
(1 + пт^4)( Л1 — Л/т + R1(1 — т R1/2 + т 2R2/6)fn = 0. (21)
Учитывая соотношение R1 = d/dt + O(h4) и отбрасывая стабилизирующий оператор, который имеет порядок малости выше, чем порядок схемы, который мы желаем показать, последнее уравнение преобразуем к следующему виду:
Л1 = fn + Mf/dt + т252f/dt2/2 + т3d3f/dt3/6) |t=rar +0(т4).
Поскольку это выражение совпадает с разложением в ряд Тэйлора до слагаемых четвертого порядка, то из него следует, что схема имеет третий порядок аппроксимации по времени. Четвертый порядок по пространственным переменным определяется использованием формулы (7) для аппроксимации оператора конвективного переноса.
Для исследования устойчивости снова применим метод Фурье. Представив решение в виде fn = pn exp(£m=1 a^j), 0 < a < 2n, j = v/—T, получим множитель перехода для схемы
p = 1 — jP1B/A — СУ-^
где
A =(1 + nPi), B =1 — P2/6, C = P2/2,
а P1, P4 даются формулами (11). Несложно проверить, что условие устойчивости | p |2< 1 записывается в виде неравенства
P2B2 < C(2A — C) или, с учетом определений параметров A, B, C, в виде
P2(1 — P2/6)2 < P2(1 + ПР4 — р2/4).
Сократив в последнем выражении неотрицательный множитель Pj2 и перенеся все в левую часть, получаем квадратный трехчлен относительно переменной Pj2. Поскольку коэффициент при старшей степени неотрицателен, максимальное значение этот трехчлен может принимать лишь на границах интервала определения переменной Pj2. Левая граница интервала — ноль, очевидно, обеспечивает выполнение нужного неравенства. Из соотношений (12) можно получить ограничение
Pi2 < Y(P4)1 /2, Y = m(m(1 + 7/243))1/2, (22)
используя которое в качестве правого интервала переменной Р-)2, получаем неравенство
(1 - 7(Р4)1/2/6)2 < (1 + пР4 - 7(Р4)1/2/4). (23)
Положим п = 72/36. Раскрывая квадрат в левой части неравенства (23) и сокращая одинаковые члены, получаем соотношение
-7(Р4)1/2/3 < -7(Р4)1/2/4,
которое, очевидно, верно. Таким образом, выбор п = 72/36 = т3(1+7/243)/36 обеспечивает устойчивость схемы (18)-(20) при любом значении т.
Сделаем некоторые замечания по данной схеме.
1. Числовой коэффициент 1 + 7/243 определяется использованием формул (7) четвертого порядка аппроксимации для конвективных слагаемых уравнения (6). Возможны другие варианты аппроксимации оператора конвективного переноса. В этом случае приведенное выше доказательство устойчивости сохраняет силу, изменится лишь данный коэффициент.
2. В многомерном случае приходится использовать факторизацию стабилизирующего оператора в уравнении (20). Не повторяя всех выкладок, укажем, что, например, для двумерных задач доказательство устойчивости сводится к доказательству неравенства (23), в котором параметр п заменен выражением п[1 + Р4хР4у/(Р4х + Р4у)], где Р4 = Р4х + Р4у соответствует расщеплению стабилизирующего оператора по пространственным переменным. Ясно, что такое неравенство тем более верно, т. е. факторизация данной схемы не нарушает ее устойчивости.
3. Факторизация не ухудшает аппроксимационных свойств схемы как по времени, так и по пространству, так как стабилизирующий оператор имеет в любом случае порядок 0(т 4).
4. Схема имеет свойство полной аппроксимации, т. е. стационарное решение не зависит от временного шага, что позволяет использовать эту схему в рамках метода установления и выбирать оптимальный шаг по времени в соответствии с критерием наиболее высокой скорости сходимости.
Естественно, схемы повышенной точности для решения гиперболических уравнений первого порядка необходимо тестировать на разрывных решениях. При этом для подавления осцилляций решения возле разрывов следует вводить искусственные диффузионные слагаемые. Отметим, что используемый здесь способ временного центрирования, характерный для схем типа Рунге — Кутты, может быть применен для построения схем повышенной точности по времени [6] для решения параболических уравнений без конвективных слагаемых. Рассмотрим уравнение теплопроводности
т т
д//д*= ЕЕ 2//5жйдх, аы = , 1 < т < 3.
1=1 к=1
Для того чтобы это уравнение было параболическим, потребуем выполнения условий: в одномерном случае а11 > 0, в двумерном — а12 < а11а22, в трехмерном — < а11а22,
а13 < а11а33, «32 < а33а22.
Введем диффузионные операторы = ^2т=1 агг А+ А--, и = Udiag+Z] 1=^к=г А0Ак. Рассмотрим двухшаговую схему
(1 - АтЦ^Х/* - /п)3/(2т) = и/п,
(24)
(1 - тАиЛчї)(/”+1 - /")/т + тАиЛч,(/* - /")3/2 = и(/“/4 + /*3/4). (25)
Здесь А — параметр, который определяется при анализе устойчивости схемы. Оказывается, что выбор А =т/4 гарантирует [6] устойчивость данной схемы, т. е. в одномерном случае следует использовать А = 1/4, в дву- и трехмерном случаях — А =1/2 и А =3/4.
Данная схема имеет второй порядок как по времени, так и по пространственным переменным. Достоинством схемы является то, что факторизация стабилизирующего оператора по пространственным переменным не ухудшает ее аппроксимационные свойства и оставляет устойчивой. Второй порядок имеет место несмотря на несогласованность стабилизирующего оператора с правой частью, а именно, этот оператор не содержат смешанных производных.
Однако, как показывает анализ Фурье, объединение приведенных выше аппроксимаций для диффузионных членов и аппроксимаций (9)-(10) или (18)-(20) не позволяет получить устойчивую схему. Попытки расчетов оказались безуспешными также в случае, когда диффузия вводилась максимально устойчивым, согласно обычным представлениям, образом, а именно, в уравнения (18) - (20) в правую часть добавлялись явные диффузионные слагаемые, а в уравнение (20) соответствующая добавка с коэффициентом 1 введена в стабилизирующий оператор (обычно в одномерных задачах достаточно использовать коэффициент 1/2).
В дальнейших поисках было найдено обобщение схемы (18) - (20), в котором аппроксимационные свойства по времени снижены, однако имеет место абсолютная устойчивость. Итак, рассмотрим т-мерное уравнение переноса с эллиптическим оператором в правой части (ориентируясь на включение в схему (18)-(20) искусственной вязкости, будем полагать, что этот оператор не содержит смешанных производных):
т
д//ді + ^ агд//дх = Р/. (26)
1=1
Введя соответствующий диффузионный разностный оператор и, построим трехшаговую схему:
(/* - /")3/(2т) + Я/" = и/", (27)
(1 - £т2Я - ти/2)(/** - /")3/(2т) + Яі(/" + /*)/2 = и/", (28)
(1 + пт4Я - £т2Я - ти/2)(/"+1 - /")/т + Яі(/"/4 + /**3/4) = и/". (29)
Здесь Я1, Я2, Я4 — разностные операторы, определяемые формулами (7), (8), п = т3(1 + 7/243)/36, £ — параметр, который определим в анализе устойчивости.
Исключим разностную функцию промежуточного слоя /*, умножив уравнение (27) на оператор -Я1т/3 и сложив его с уравнением (28). Получим уравнение
(1 - £т2Я2 - ти/2)(/** - /")3/(2т) + Я/" - тЯ2/"/3 = 0.
Исключим разностную функцию промежуточного слоя /**, умножив последнее уравнение на оператор -Я1т/2 и сложив его с (29), умноженным на оператор 1 - £т2Я2 - ти/2; в результате получим уравнение
(1 - ти/2 - £т2Й2 + пт4ВД(1 - £т2Я2 - ти/2)(/"+1 - /")/т+
+ (Я1 - и)(1 - £т2Я2 - ти/2 - т^/2 + т2Р?/6)/га = 0. (30)
Используя для исследования устойчивости метод Фурье, получаем множитель перехода для схемы
Р = 1 - ,7Р1В/А - СУА
где
А = (1 + £Р2 + ПР4 + У/2)(1 + + ^/2), В =1 + £Р2 - Р2/6,
С = Р2/2 + V(1 + £Р2 + У/2 - Р2/6),
Р\, Р4 даются формулами (11), V > 0 — спектральный образ оператора -тЦ. Несложно проверить, что условие устойчивости | р |2< 1 записывается в виде неравенства
С(2А - С) > Р^В2 или, с учетом определений параметров А, В, С,
[Р2/2 + V(1 + £Р2 + У/2 - Р2/6)][2(1 + £Р2 + пР4 + У/2)(1 + £Р2 + У/2)-
-Р2/2 - V(1 + £Р2 + У/2 - Р2/6)] - Р2[1 + £Р2 - Р2/6]2 > 0. (31)
При введении величины О = (1 + £Р2 + пР4)(1 + £Р2 + V/2) - Р^/12 - Рр/4, то (31) можно переписать в следующем виде:
^(1 + £Р2 + V/2 - Р2/6)]О + Р2[О - (1 + £Р2 - Р?)2] > 0.
Таким образом, если величина О и фигурирующие в последнем неравенстве выражения в квадратных скобках неотрицательны, то оно, очевидно, верно. Если выбрать £ = в/6, где в = т250/243 — параметр, связывающий, ввиду неравенств (12), величины Р2, Р1, то выражение в первой квадратной скобке (31) неотрицательно. Выражение во второй квадратной скобке заведомо меньше О, и, следовательно, если оно неотрицательно, то и О неотрицательно, и нужное нам неравенство верно. Итак, докажем неотрицательность выражения во второй квадратной скобке, причем запишем его, используя тождество (1 + £Р2 - Р2/6)2 = (1 + £Р2)(1 + £Р2 - Р2/3) + (Р2/6)2. В результате получаем следующее неравенство:
(72/36Р4 + Р2/3)(1 + £Р2) + (1 + 72/36Р4 + ^^/2 + P12V/12 - Р*/4 - (Р2/6)2 > 0.
Отбросим множитель (1+£Р2) возле первой скобки, тем самым лишь уменьшая, возможно, левую часть искомого неравенства. Поскольку слагаемые в ней в совокупности с двумя вычитаемыми в конце соотношения, с учетом (22), образуют неотрицательное выражение и остальные слагаемые неотрицательны, то неравенство (31) справедливо.
Итак, абсолютно устойчивая схема с диффузионными слагаемыми построена. Однако включение диффузии достигнуто за счет понижения аппроксимационных свойств схемы. Варианты построения диффузионных операторов и способы избежать понижения аппрок-симационных свойств схемы обсуждаются далее при численном исследовании одного точного решения уравнения Бюргерса. Рассмотрим нелинейное уравнение
ди/д£ + 5(м2/2)/5ж = 0.
(32)
Построим трехшаговую схему:
(и* - ига)3/(2т) + Лі(мга)2/2 = игамга,
(33)
(1 - £т2Д* - ти*/2)(и** - )3/(2т) + Ді[(ига)2 + (и*)2]/4 = и*ип,
(34)
Здесь п = (1 + 7/243)/36, £ = (1 + 7/243)/6, Л1 = Д0(1 - 5 5+/6), Д2, Д4, и — нелиней-
ные разностные операторы, верхний индекс которых отмечает временной слой, на котором вычисляются фигурирующие в них коэффициенты. Например, Д* = (1/2)4Д-(и* + и*+1)4Д+Д-Д+. Операторы и с различными индексами предназначены для подавления осцилляций решения возле скачков и строятся на основе подхода [7], в частности:
Другие верхние индексы у данного оператора повлекут использование таких же индексов в правой части приведенных формул. Как обычно, штшо^ж,у) = sign(ж)x тах[0,min(sign(ж)y, | у |)]. Отметим, что для “гладкого” решения вне экстремумов имеет место 7г = 1, и, нетрудно видеть, диффузионные слагаемые исчезают. Возле разрывов коэффициенты 7г принимают значения, меньшие единицы. Малая константа е введена для того, чтобы не включать диффузионные слагаемые низкого порядка вблизи “гладких” экстремумов. Этот прием предложен в [8] для построения ТУБ-схем, т. е. схем с ограниченной полной вариацией (вообще говоря, возрастающей, в отличие от ТУБ-схем).
В приведенных выше формулах коэффициенты являются индикаторами включения искусственной диффузии. Поэтому их удобно использовать при построении стабилизирующего оператора Я2, который, как показывает анализ схемы (18) - (20), может отсутствовать при нулевой искусственной диффузии и должен включаться, когда она появляется.
В итоге использовалась формула
здесь с > 0 — эмпирическая константа, влияющая на размеры области, где включен стабилизирующий оператор Д2. Отметим, что на гладком решении е = 0, и схема имеет третий порядок аппроксимации по времени. Максимальное значение параметра е равно, как и следует в соответствии с анализом устойчивости схемы (27)-(29), (250/243)/6. В пробных расчетах найдено, что величина с = 0.125 обеспечивает приемлемые результаты.
3. Результаты расчетов
Рассмотрим результаты решения уравнения (32) с начальным распределением, представляющим собой прямоугольник единичной высоты на нулевом фоне. При этом передний скачок распространяется со скоростью 0.5, второй скачок, имеющийся в начальный момент времени, превращается в линейное распределение, левая граница которого, смыкающаяся
7* = [тіпто^| 5+и* |, ^5 и*2А + є) + тіпто^| 5 и* |, 5і-15+м*2А + є)]/ /(| 5+и* | + | 5-и* | +10-8), є = 10(итах* - итіп*)//2, 1 < г < I.
тах *
Д* = 0.25Д-(и* + <+і)2ет/2Д+, Єі+і/2 = тіп(2 - 7* - ^*+і,с)(250/243)/(6с),
Рис. 1. Уравнение Бюргерса; число Куранта равно 0.5 (а), 1 (б), 2 (в).
Рис. 2. Уравнение переноса; число Куранта равно 0.5 (а), 1 (б), 2 (в).
с нулевым фоном, неподвижна, а правая, на которой значение и равно единице, движется с единичной скоростью вправо, постепенно догоняя разрыв.
На рис. квадратиками изображено численное, точками — точное решения в момент времени t = 0.65 для чисел Куранта 0.5, 1, 2. Используется равномерная сетка, содержащая 108 узлов по переменной х. Произведено 140, 70, 35 шагов схемы (33)-(35) численного интегрирования уравнения (32). При дальнейшем увеличении временного шага решение значительно ухудшается.
Для тестирования численных методов популярной является задача линейного переноса возмущений различных форм. Эта задача решалась численным интегрированием одномерного уравнения переноса с помощью схемы (33) - (35), в которой в качестве потоков вместо и2/2 используется и, в стабилизирующих операторах и заменяется на 1, однако диссипативные операторы и адаптивные алгоритмы сохранены для возможности моделирования возмущений с разрывами функции или ее производных. Начальное распределение включает последовательно треугольное, параболическое и прямоугольное возмущения единичной высоты.
На рис. 2 квадратиками изображено численное, точками — точное решение в момент времени t = 0.33 для чисел Куранта 0.5, 1, 2. Используется равномерная сетка, включающая 217 узлов. Произведено 144, 72, 36 шагов схемы. При дальнейшем увеличении временного шага решение, как и для уравнения Бюргерса, значительно ухудшается.
Список литературы
[1] РусАнов В. В. Разностные схемы третьего порядка точности для сквозного расчета разрывных решений. Докл. АН СССР, 180, №6, 1968, 1303-1305.
[2] Burstein S.Z., Mirin A.A. Third order difference methods for hyperbolic equations. J. Comput. Phys, 5, No. 3, 1970, 547-571.
[3] Kutler P., Lomax H., Warming R. Computation of space shuttle flow fields using noncentered finite-difference schemes. AIAA Paper, No. 72-193, 1972.
[4] Толстых А.И. Компактные разностные схемы и их приложения к проблемам аэрогидродинамики. Наука, М., 1990.
[5] Новиков Е.А. Явные методы для жестких систем. Наука, Новосибирск, 1997.
[6] Пинчуков В.И. Численные методы аэрогидромеханики высоких порядков аппроксимации. Новосиб. гос. ун-т, 1997.
[7] Пинчуков В.И. Энтропийные и квазиэнтропийные условия для схем высоких порядков аппроксимации. Вычислит. технологии, 2, №6, 1997, 71-87.
[8] Shu C.W. TVB Uniformly High-Order Schemes for Conservation Laws. Math. Comput., 49, No. 179, 1987, 105-121.
Поступила в редакцию 18 июня 1998 г., в переработанном виде 4 декабря 1998 г.
