Компактные схемы для систем уравнений второго порядка с конвективными членами Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 3, № 1, 1998

КОМПАКТНЫЕ СХЕМЫ ДЛЯ СИСТЕМ УРАВНЕНИЙ ВТОРОГО ПОРЯДКА С КОНВЕКТИВНЫМИ ЧЛЕНАМИ *

В. И. ПААСОНЕН Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: paas@net.ict.nsc.ru

In present paper the properties of special approximations for dissipative terms are investigated, and on this basis the compact difference fourth order accuracy schemes for many-dimensional systems of differential equations with convective and dissipative terms without mixed derivatives are constructed. The factors of the equations are assumed by variables, and they can be depended also on a solution. Alongside with elliptic systems of equations the evolutionary systems are considered also. The order of approximation on time is second for the parabolic equations and it is fourth for the hyperbolic equations.

Введение

Компактными принято называть разностные схемы, которые имеют повышенный порядок аппроксимации, но пишутся на шаблоне, несущественно отличающемся от традиционного для данного уравнения. Для конвективно-диффузионных уравнений это обычно схемы третьего-четвертого порядка аппроксимации, шаблон которых представляет собой m-мерный параллелепипед с размерами ребер в два пространственных шага по каждому из m координатных направлений, называемый иначе m-мерным ящиком. В отличие от схем повышенной точности на многоточечных шаблонах для компактных схем аппроксимация является улучшенной не на любых гладких функциях, а именно на решениях дифференциального уравнения, так как в этом случае она достигается путем использования следствий уравнения, полученных его дифференцированием. Начало этим исследованиям положил Микеладзе задолго до появления самого термина "компактная схема" в работах 1934-41 гг. (см. [3] и приведенную там библиографию), где были созданы первые схемы четвертого порядка точности для уравнений с постоянными коэффициентами на шаблоне типа "ящик".

Самарский в одной из своих работ [5] распространил этот результат на уравнение теплопроводности с самосопряженным оператором, использовав в качестве исходного элемен-

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №97-01-00819.

© В. И. Паасонен, 1998.

та традиционный по форме разностный аналог для "сложной" второй производной

(здесь и далее Д ± — разделенные односторонние разности), имеющий специальное разложение по степеням шага сетки к:

Выполнение этого свойства для каждого из координатных направлений позволяет в случае переменных коэффициентов с помощью уравнения выразить "сложные четвертые" производные по одной переменной через смешанные порядка 2 + 2, обеспечивая при создании компактной схемы аналогию со случаем постоянных коэффициентов. Заметим, что представленное выше разложение не является тривиальным. Вид разложения определяется зависимостью коэффициентов разностного оператора от ряда значений функции а дифференциального оператора Ь, называемой, согласно терминологии Самарского, шаблонным функционалом. Так, например, простой выбор

не дает нужного разложения. В работе [5] дан следующий пример более тонкой зависимости довольно необычной формы

обеспечивающей нужное разложение.

Для уравнения колебаний, вероятно, впервые компактная схема четвертого порядка точности (по пространству и времени) построена и исследована в работе [2]. В [4] даны обобщения компактных схем для уравнений трех основных типов с переменными коэффициентами без смешанных и младших производных. Наконец, вопросы построения и исследования компактных схем для уравнений в неортогональных координатах рассматривались, например, в работах [2, 6]. Отметим также монографии [1, 7], посвященные компактным схемам и ценные также приведенной в них обширной библиографией.

В течение нескольких десятилетий упомянутые компактные схемы и многочисленные их обобщения и модификации назывались "схемами повышенной точности, построенными без расширения шаблона". Краткое наименование они приобрели благодаря нашим англоязычным коллегам, предложившим к хорошо разработанной теме удачный и выразительный термин. Вплоть до настоящего времени и у нас, и на Западе к повышенной точности проявляется большой интерес, причем в связи с изменившимися возможностями вычислительной техники многое открывается несколько по-новому — с акцентом на технологический аспект. Поэтому становится понятной попытка автора записать компактные схемы для уравнений в как можно более общей форме. Результаты такого исследования приводятся в настоящей работе, где строятся компактные схемы для многомерных систем

или

1

уравнений трех основных типов, содержащих диссипативные и конвективные члены в произвольном виде. Исследование представляет собой развитие результатов работы [4] и, в некотором смысле, подводит итог многолетним усилиям многих авторов, обобщая результаты отдельных построений для частных вариантов уравнений и систем с постоянными и переменными коэффициентами на широкий класс систем уравнений с переменными коэффициентами и нелинейностью.

1. О структуре элементарных разностных операторов

Прежде всего рассмотрим вопрос о том, каким образом следует определить шаблонный функционал в разностном операторе

с(х+1)д+ - с(х - 2)д-

Ли(х) = ^-^---^-^-и(х), (1.1)

Ъ

чтобы в окрестности точки х имело место разложение специального вида

Ъ2 д { д \

Ли = Ьи + —Ь(а-1Ь)и + 0(Ъ4), Ь = д ад . (1.2)

12 дх \ дх/

При этом в отличие от аналогичного разложения, приведенного во введении, здесь предполагается, что () — не скалярная функция, а матрица с достаточно гладкими элементами, положительно определенная и имеющая обратную матрицу, также достаточно гладкую в области изменения независимой переменной х.

Предположим, что шаблонный функционал симметричен относительно средины ячейки и аппроксимирует значение в ней функции а. Тогда его разложение в окрестности

± Ъ Ъ

точки х ± 2 содержит лишь четные степени Ъ, и следовательно его можно записать в виде

с(х ± 2) = а(х ± Ъ) + Ъ2 ± £) + °(Ъ4)-

Путем выбора функционала о можно, очевидно, влиять на коэффициенты оператора (1.1), и таким образом задача сводится к определению условий на о, при выполнении которых разложение (1.2) имело бы место. Тот же функционал в окрестности центральной точки х шаблона может быть представлен в виде суммы симметричной и антисимметричной составляющих

с(х ± Ъ)= А ± Ъ В + 0(Ъ4),

где (у функций и производных в правых частях подразумевается аргумент х)

Ъ2 Ъ2

А = а + От (За" + 2о), В = а' + Ъ- (а'" + 2о'). 24 24 '

Тогда для оператора (1.1) получается разложение

а'' ди

+ 0(Ъ4),

Ъ2 д

Ли = Ьи + 12 дх

д3и ' д2и / а!'\ ди дх3 дх2 \ 2 ) дх

которое, как нетрудно убедиться непосредственной проверкой, совпадает с (1.2) тогда и только тогда, когда выполняется равенство

а" а" а = у + а(а-1)'а' + О(к2) или а =у - а'а-1а' + О(к2).

(1.3)

Отсюда следует, что (1.3) является необходимым и достаточным критерием выполнения для оператора (1.1) разложения (1.2). Очевидно, критерий (1.3) эквивалентен тому, что с точностью до членов О (к4) шаблонный функционал представляется в виде

с х ±

к2 а'' ' -1 '

а +12\Т - аа а

х±

(1.4)

Разумеется, при выводе учтен запрет на перестановку сомножителей, и следовательно формулы (1.3), (1.4) справедливы как для скалярного, так и для матричного оператора Л. Разложение (1.4) справедливо, например, для шаблонного функционала

с ^х ± = 6 |а-1 ^х ± к) + 4а-1(х ± ^ + а-1(х)

полученного для скалярного случая в работе [5].

Очевидно, существует бесчисленное множество трехточечных функционалов, удовлетворяющих условию (1.4). Это, например, шаблонные функционалы

с Гх ± к | = со Гх ± к | — 112[а(х ± к) — а(х)]а-1 Гх ± ^ | [а(х ± к) — а(х)],

х±

к

со

где

к\ 1

х ± 2 ) — 2^[а(х ± к) — а(х)][а(х ± к) + а(х)]-1[а(х ± к) — а(х)],

к

. а(х ± к) + 4а I х ± — I + а(х)

си *± 2)=-^—,

а также их любые линейные комбинации с суммой коэффициентов, равной единице. Здесь также нельзя менять порядок сомножителей в том случае, когда а — матрица.

Из (1.4) следует, что никакая линейная комбинация значений функции а не приводит к (1.2), так как ее разложение не может содержать нелинейных членов, присутствующих в (1.4).

Из (1.4) вытекает также невозможность в общем случае обойтись без точки с полуцелым индексом, если, конечно, не привлекать при этом более удаленных точек (см. [4]), так как вторую производную функции а(х) невозможно аппроксимировать лишь по двум точкам. Если же функция а и ее производные не зависят от решения и их зависимость от х задана аналитически, то в качестве шаблонного функционала можно использовать непосредственно выражение (1.4).

Из (1.4) можно также извлечь важные частные случаи для уравнений с постоянными коэффициентами в цилиндрической (а(х) = х) и сферической (а(х) = х2) системах координат. Вот соответствующие этим случаям шаблонные функционалы с произволом до величин О (к4):

с

2

к\ / к\ к2 1 / ± к\ / ± к\2 к2

С11 х± ^ = 1х± V — СЧх± ^ = 1х± 2] —

х± 2,

принадлежащие параметрическим семействам, иным способом полученным в работе [4].

Следует заметить, что в действительности для построения компактной схемы достаточно представимости оператора Л в более свободной, чем (1.2), форме

Ли = Ьи + к2М (Ьи) + О(к4), (1.5)

где М — произвольный дифференциальный оператор второго порядка. Достаточность (1.5) становится очевидной, если учесть, что именно внутреннее вхождение Ь во второе слагаемое позволяет выразить члены погрешности с четвертыми производными по одноименным переменным через четвертые смешанные.

Приведем пример оператора, разложение которого имеет вид (1.5), но не совпадает с (1.2). Пусть Л1, Л2, Л3 — одномерные разностные операторы вида (1.1) с коэффициентами соответственно

1 к к 2 к а( + к) + а( )

+ к)= а(х + а-'(х + ^ а-'(х)а(х + к

Тогда нетрудно показать, что линейная комбинация

л = 1Л1 — ^Л2 — -Л3

3412

удовлетворяет условию (1.5), в котором в качестве оператора М выступает

_ 5 ( д(а-1) \ 52 М = — а \ +

дх \ дх ) дх2

Форма разложения (1.5) в сравнении с (1.2) предоставляет дополнительную степень свободы при построении компактных схем, которая может быть использована для придания им тех или иных свойств. Однако этот предмет требует специального исследования.

2. Исходная компактная схема

Рассматривается эллиптическая система уравнений:

Ки = Ди — /, (2.1)

а также соответствующие ей параболическая и гиперболическая системы уравнений, получающиеся из (2.1) включением в ее левую часть первых или вторых производных по времени. Здесь и и / — вектор-функции, а дифференциально-матричные операторы Д и К представляют собой диссипативные и конвективные члены1

хЗдесь и далее знак суммирования с опущенными пределами означает суммирование от 1 до т.

д ( д \ Дм = V" Ь^м, и = — а—— ) ,

дж^ дxi)

с\ т. г\ г\

V о дм ^ о дР

Км = = > т—, Bi = —.

дxi джi дм

Переменные коэффициенты ai,Ьi,Bi — суть матрицы, а р — вектор-функции. И те и другие могут зависеть также от решения м.

Рассматриваемые системы уравнений формально можно записать в форме эллиптической системы без младших производных

Дм = Д, (2.2)

где Д есть сумма правой части f, конвективных членов и производных по времени ¿, если они присутствуют в исходой системе.

Следуя работе [4], запишем компактную схему для уравнения (2.2). Для этого предположим, что Л¿ — какие угодно трехточечные одномерные аппроксимации операторов для которых выполнен критерий (1.4) по соответствующей координате, обеспечивающий разложение

Л2 д / д \

Л = и + Т2¿¿(«Г%) + 0(Л4) и = дж- ^ , (2.3)

где Л — шаг сетки в г-м координатном направлении. Тогда нижеследующая компактная схема аппроксимирует уравнение (2.2) с погрешностью 0(Л4) (Л = шах(Л^)) и благодаря разложению (2.3) строится аналогично схеме Микеладзе:

Лм = СД, (2.4)

где

Л = £ ЬЛ (к + М а-Ч"1 £ Ь3 Л^ , С = Е ЬДДа-Ч^Е),

а Е — тождественный оператор.

В линейном случае можно показать, что для размерностей т = 1, 2, 3 схема (2.4) при некоторых не очень обременительных ограничениях сверху и снизу для отношений шагов hi/hj является схемой с положительными коэффициентами со всеми вытекающими отсюда благоприятными последствиями (выполнение принципа максимума и сходимость в равномерной норме). Для больших размерностей анализ коэффициентов схемы становится затруднительным.

По построению наш объект (2.4) еще не является разностной схемой, так как содержит в Д производные по пространству и, может быть, по времени. Для построения компактных схем их следует корректно, без потери точности, выразить через разностные отношения, не выходя за пределы нашего "ящика".

3. Случай дивергентной формы записи конвективных членов

Рассмотрим систему уравнений (2.1) в том случае, когда конвективные члены записаны в дивергентной форме. Как отмечалось выше, систему (2.1) можно представить в виде (2.2) с правой частью, включающей конвективные члены

Ди = ^ = / + Ки = / + £ (3.1)

и формально выписать для нее схему (2.4). Заменяя в правой части (2.4) конвективные члены центрально-разностными аппроксимациями и учитывая возникающие при этом погрешности, можно с точностью до членов порядка О (к4) записать

Ли = С/ + Ар + £ М-+ &>), (3.2)

где

Т — Г-1

л \ Л л л Т — Г г

Ар = > Агрг, Аг =--- — центральные разности,

2кг

Тг — оператор сдвига: Тги>(хг) = и>(хг + кг),

(а"1 £ ^ '

Яг р = Ьг ( аг

3=г

( др • \ д 3р' 6р = — 2 . (3.3)

В выражении (3.3) собраны производные третьего порядка, которые, очевидно, без преобразования не аппроксимируются трехточечными разностными отношениями. Но первое слагаемое (3.3) удается выразить через второе:

ЬгЬг (а-Ч-1 др) = др: + Ьг (^г + Я^г), (3.4)

где

Я = А

Яг дхг

д(-1Ь-1) д

д(Ь-1) д2

дх дх

а третья производная представляется в виде суммы

д3р д ^ , ы

Я2 =__,

дх дх2

дх3 дх где

в = др, д1 = А

ди дх

Так как в силу уравнения (3.1)

(В:-1Ьги) + ЯЯ (3.5)

А (Вг-1)аг_д_

дх дх

ььи = / + Ц + Е ® —ь'ь и)

3=г

то на гладких его решениях выражение (3.5) можно переписать в виде

д

дж3

= Д1м + ^ + дг3р - Л2м,

(3.6)

где Л1 определено выше, а

Л=I- (с- £Ь, ,

д

• = & (с-Е >•

о3р = Ь

дж дж,

дж.

В:а-1Ь-1.

Заметим, что выражения Л-, 3-, в конечном итоге входят в (3.2) с множителями Л2, поэтому их достаточно аппроксимировать со вторым порядком. Это можно сделать не выходя за пределы т-мерного ящика, так как каждое из них состоит из слагаемых простых типов, которые аппроксимируются естественным образом:

д /^д^

Д-:

Т-1/2а)

д д^

Д:(аД, ад),

д / , дд ди> ^— а

дж^ \ дж- дж-

Д-

д2ад дж2

дад дж-

д

-— « -— (аи,^ Д:(ал,

здесь оператор Т-1/2 осуществляет сдвиг аргумента на полшага по г-й координате, Д- (без знака) — оператор центральной разности, а Д±- — операторы направленных разделенных разностей.

Если теперь пересчитать (3.3) с помощью (3.4), (3.6) и подставить результат в (3.2), а затем выполнить замену всех производных разностными отношениями по приведенным выше образцам, то получим компактную схему для эллиптического уравнения (3.1) на основе дивергентной записи конвективных членов:

Л м = Др + С f,

(3.7)

где

Л2

Л = Л + Х)Т2;Л:, Л « Л -

К

Д = Д + Ет2 О?", О?" « Ь: (30 + 3 + 3?) - 3

12

К2

С = с _1С с ~ •

С = С — / у 12 •¿, ~

с

с

2

4. Случай недивергентной формы записи конвективных членов

Недивергентные схемы также могут представлять интерес в некоторых задачях, например, в качестве предварительного шага в схемах предиктор — корректор. Применяя описанный выше способ исключения погрешности в случае системы с недивергентной записью конвективных членов

Ки = Ди - /, Ки = У Bi(4.1)

вместо (3.2) получим равенство

к?

Ли = с/ + £ В;дги + £ 12 (ад°и + е

; и +

(4.2)

с несколько иными выражениями для компенсирующих слагаемых

= ¿л

3= 5и \

£ги = ( а"1 б"1 В;- 2В

5 3и 'дж3'

Первое слагаемое (4.3) в силу уравнения (4.1) выражается в виде

(4.3)

ЬА (а-16"1ВгдЖ^ = 6г(#/ + д1и + д?и),

(4.4)

где

А

дж;

% \ 5и СгВ; + а^—а; —

^ 1Вз дЖ^- Ьз ¿зи

3=» ^ 3

д _ _ _ _ _ _ / = дЖ (с/), с = 1В<а_ 1 1, ф = а_ 1 1 В;а_ 1,

а второе в виде

где

2 В;

д3и 5ЖЗ

-2В; (5?/ + Я1и + Д?и) ,

(4.5)

Д1и

А В, + а А ди

5Ж; 7 5Ж;

Д?и = £

ди

(В3 дЖ- - Ь3 ¿3и

д

/ = дж- (А/), А = а_1ь_1.

Подставляяя далее (4.4), (4.5) в (4.3), а затем результат в (4.2) и заменяя производные разностными отношениями, как описано выше, получим компактную схему для эллиптического уравнения (4.1) в случае недивергентной записи младших производных:

Л и = Ди + С/,

(4.6)

где

Л = Л + £_| (Ь-3 : - 2В- Л-) , О? - « + о1 + Я1 Л- « + Л2,

Дм = £ В-Д-м, С = С + £ К! (Ь-? - 2В-.?2), • ~

а Л, С являются операторами компактной схемы (2.4) для соответствующего уравнения без конвективных членов.

5. Компактные схемы для эволюционных систем

При построении компактной схемы для параболического уравнения воспользуемся стандартным приемом, формально заменив правую часть f на f + м^:

Км = Дм - Д, Д = f + м4, (5.1)

причем ограничимся лишь случаем дивергентной записи конвективных членов. На разностном уровне такую операцию нужно провести с компактной схемой (3.7). Но перед этим представляется целесообразным операторы Л, Д, С схемы (3.7) записать в виде сумм главных частей (членов нулевого порядка) и компенсирующих членов второго порядка малости:

Л = ^ Ь-Л: + ^ Т2Л¿, Л: = Ь:Л-а_1Ь_^ Ь,Л, + Л-,

Д = £ Д- + £ _20¿'

С = Е + £ К2«¿' "С: = ЬiЛi(а_!Ь_!E) -

Пусть введена дискретизация по времени = тп, где т — временной шаг. Тогда для достижения точности 0(т2 + К4) главные члены необходимо аппроксимировать симметрично по времени, а поправочные допустимо аппроксимировать как угодно, например, на нижнем временном слое. В результате получается компактная схема для параболического уравнения (5.1):

(Е + £ Т|С:) ^^ + £ - £ Ь-Л-м^ = ф., (5.2)

где

/ ¿.+ 1 + /"га _2 \ _2

ф. = 2 7 + Е ^С-Г] + Е ^ (Л-мп - 4Л

объединяет результаты операций над правой частью и функциями на нижнем слое. Схема (5.2) по форме совершенно аналогична традиционным неявным схемам, для нее применимы обычные методы линеаризации в нелинейном случае и методы приближенной факторизации, использование которых сохраняет порядок аппроксимации и приводит к ряду одномерных разностных уравнений, реализуемых трехточечными матричными прогонками.

Гиперболическую систему второго порядка можно представить как эллиптическую систему с измененной правой частью (здесь предполагается дивергентное представление конвективных членов)

Ки = Ди - Д, Д = / + иЙ. (5.3)

Ввиду того, что (3.7) — схема четвертого порядка аппроксимации по к для эллиптического уравнения, и в силу следствия уравнения (5.3), полученного дифференцированием его дважды по времени, с точностью 0(т4 + к4) получим равенство, справедливое на достаточно гладких решениях системы (5.3):

, / т2 д2 \

Лип = Д+ С / + Аип - - — (и«)га ) , (5.4)

где А — обычный разностный аналог второй производной по времени

„ „ ип+1 - 2ип + ип_1

Аип =-2-.

т 2

Аппроксимируя в (5.4) иц с привлечением уравнения (5.3) со вторым порядком

иы = Ди - К и - / И ^ (Ь^и - ДгРг) - /,

а внешний оператор двойного дифференцирования по времени — оператором А и пренебрегая членами четвертого порядка малости, получим трехслойную компактную схему для гиперболической системы (5.3), имеющую погрешность 0(к4 + т4):

т2

С(Аип) - —А^1 + Дрп = Фга, (5.5)

где 2

Фп = Ли™ - ^(С + 12А) /п, V™ = ^ (ЬДгип - дгРп).

После приближенной факторизации оператора в (5.5) получится компактная схема, эквивалентная (5.5) по порядку аппроксимации и реализуемая трехточечными прогонками. При недивергентной аппроксимации конвективных членов построение проводится аналогично исходя из схемы (4.6).

Список литературы

[1] Вллиуллин А. Н. Схемы повышенной точности для задач математической физики. Изд-во Новосибирского гос. ун-та, 1973.

[2] ВАЛИУЛЛИН А.Н., ПААСОНЕН В. И. Схемы повышенной точности для параболических и эллиптических уравнений со смешанной производной. Численные методы механики сплошной среды, 15, №2, 1984, 36-41.

[2] ВАЛИУЛЛИН А. Н., ПААСОНЕН В. И. Экономичные разностные схемы повышенного порядка точности для многомерного уравнения колебаний. Там же, 1, №1, 1970, 34-47.

[3] МИКЕЛАДЗЕ Ш. Е. О численном интегрировании уравнений эллиптического и параболического типов. Изв. АН СССР. Сер. матем., 5, №1, 1941, 57-74.

[4] ПААСОНЕН В. И. Обобщение методов повышенной точности для нелинейных уравнений 2-го порядка в ортогональных системах координат. Численные методы механики сплошной среды, 8, №2, 1977, 94-99.

[5] Самарский А. А. Схемы повышенного порядка точности для многомерного уравнения теплопроводности. Журн. вычисл. матем. и матем. физики, 3, №5, 1963, 812-840.

[6] Самарский А. А. Теория разностных схем. Наука, М., 1983.

[7] Толстых А. И. Компактные разностные схемы и их применение в задачах аэрогидродинамики. Наука, М., 1990.

Поступила в редакцию 6 января 1998 г.