Введение

          В настоящей работе приведены результаты сравнительного анализа новой и некоторых известных конечно-разностных схем для решения задачи Коши для линейного двумерного уравнения переноса с финитными начальными условиями. Работа является непосредственным продолжением публикаций [1, 2, 3] для одномерного уравнения.

          Для тестирования использован построенный класс точных решений двумерной задачи. Класс определен одной произвольной функцией двух действительных переменных и четырьмя произвольными дифференцируемыми функциями одной действительной переменной.

Точные и численные решения получены в пространственно – временной области достаточно больших размеров. Её размеры выбраны из предположения, что время, за которое финитное начальное распределение пробегает 17 (9 по одной координате) своих длин, является достаточным для проявления свойств вычислительного метода. Для определения ошибок методов использованы конечномерные аналоги норм пространств , , . Предложенная система представления информации об ошибках и визуализации численного решения позволяет сравнить любую новую схему с известными схемами и в единой манере определить ее характеристики.

Авторы благодарны А.П.Фаворскому за полезные обсуждения и консультации.

§ 1. Постановка задачи и класс точных решений задачи Коши.

Пусть в задаче Коши для линейного двумерного уравнения переноса

, , ,                                    (1)

, ,                                                              

с финитным начальным условием ,

где , , , скорости ,  заданы следующим образом

,

,

где , , , .

          Уравнение (1) можно записать в виде:

,

где  – полная производная  по времени вдоль характеристик, задаваемых уравнениями

,  .

          Хорошо видно, что при переносе в несжимаемой среде, в которой , начальный профиль решения сохраняется.

          В данном случае задача (1) имеет следующее решение.

.

          Данный класс точных решений определен одной произвольной функцией двух переменных , , , задающей начальное условие , и четырьмя произвольными дифференцируемыми функциями одной переменной , , , , , , определяющими поле скоростей. Функция  задает начальный профиль. Функции ,  определяют траекторию движения центра решения  в плоскости . Функции ,  определяют законы вращения начального профиля  относительно центра.

          В работе использованы три типа функций .

Ниже на рисунках представлены виды точных решений. Слева показаны линии уровня, справа – видовая проекция.

1.     Конус

,

2.     Конус без сектора

3.     Буква «М»

, , , , ; , .

Рассмотрены следующие закономерности движения центра:

1)     центр неподвижен: , ,

2)     центр движется по прямой: , , ;

3)     центр движется по окружности радиуса , : , ,  -  время одного оборота вокруг центра.

Предусмотрен ряд вариантов движения относительно центра:

1)     фигура вращается как целое: , ,

2)     периферия фигуры вращается быстрее центра ( монотонно возрастает),

3)     периферия фигуры вращается медленнее центра ( монотонно убывает).

При этом предполагается, что угловая скорость:

1)     не зависит от времени: ,  - время обращения фигуры вокруг центра,

2)     возрастает со временем: ,  - время первого «оборота»,

3)     уменьшается со временем: , .

          Отличие численного решения, полученного с помощью какого-либо численного метода, от точного решения оценивалось с помощью конечномерных аналогов норм в пространствах  на  (интегрально по ):

и на  при  (локально по ):

§ 2. Нелинейная монотонизация схемы К. И. Бабенко.

Для численного решения задачи введена равномерная пространственно-временная сетка   

 , где , ,  – шаги разностной сетки по , ,  соответственно. Для сеточных величин использованы принятые в [4] обозначения.

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

          Рассмотрим нелинейную монотонизацию схемы К. И. Бабенко, уделив внимание знакопостоянству и перемене знака скоростей в узлах  и . Предположим, что область сохранения знака скорости содержит хотя бы три соседние точки. Пусть . Хорошо известно [6, с.185], что для получения корректной постановки задачи для гиперболических уравнений необходимо на границах области задать столько граничных условий, сколько характеристик выходят из границы в рассматриваемую пространственную область. В соответствии с этим при переменной скорости возникает ряд различных случаев.

1. Пусть , , . В этом случае характеристики наклонены вправо. Перенос решения происходит с левой границы. Тогда монотонизованная схема имеет вид:

  (2)

где , ,                                     (3)

значения функции  вычисляются следующим образом:

                                  (4)

Искомое значение  является корнем уравнения (2). Это уравнение нелинейно, так как  зависит от .

2. Пусть , , . В этом случае характеристики наклонены влево. Перенос решения происходит с правой границы. Тогда разностная схема имеет вид:

  (5)

где , ,                                         (6)

                                 (7)

Искомое значение  является корнем нелинейного уравнения (5).

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

                   (8)

             (9)

где ,

В этом случае неизвестные ,  являются решениями системы нелинейных уравнений (8, 9). Величины ,  предварительно вычисляются из уравнения (2), коэффициент  в (8) – с помощью выражений (3, 4), а ,  – из уравнения (5), коэффициент  в (9) – с помощью выражений (6, 7).

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

,                                (10)           .          (11)

Коэффициент  в (10) вычисляется из выражений (3, 4), а коэффициент  в (11) определяется соотношениями (6, 7). Неизвестная величина  определяется из первого уравнения, а  - из второго.

§ 3. Решение нелинейных уравнений.

          В случаях 1, 2, 4 одна неизвестная  () определяется из одного нелинейного уравнения. Анализ уравнений показывает, что они во всех случаях имеют единственное решение [7]. Это решение определяется следующим образом:

Случай 1. Перепишем уравнение (2) в следующем виде:

Пусть

Тогда , а  вычисляются, исходя из выражений (4), по следующим формулам

                (12)

Случай 2. Перепишем уравнение (5) в следующем виде:

Пусть .

Тогда , а  вычисляются подобно (12), исходя из (7):

(13)

Случай 3. Уравнения (8), (9) можно записать в виде

где , , ,

,

.

Далее необходимо определить неизвестные ,  либо явно [7], либо итерационным способом.

Случай 4. Запишем уравнения (10, 11) в виде

,

,

, где  определены в (12), а

 .

, где  определены в (13), а

.

§ 4. Организация вычислительного эксперимента.

Шаг по времени  выбран так, чтобы максимальное значение числа Куранта  () было ограничено сверху величиной  (). Для всех приведенных расчетов 0.25, 10. Пространственные шаги ,  выбраны единичными. Область вычисления содержит  узлов. Начальный профиль при движении по заданной прямой на конечный момент времени пробегает 17 (9 по одному направлению) своих длин.

Точность вычислений определена через нормы разности точного и приближенного решений вида:

,, .

Эти нормы дают ошибку численного решения интегрально по всему промежутку . Кроме того, применены их аналоги, соответствующие моменту времени :

          Для сравнительного анализа новой схемы с известными расчеты проведены также по схеме с «лимитерами» [8], по схеме с направленными разностями и по схеме П. Лакса [9] (вычислительные алгоритмы содержатся в приложении).

§ 5. Описание результатов численных расчетов.

Результаты численных расчетов представлены в виде рисунков и таблиц. На каждом из рисунков представлены численное решение  (справа) и его линии уровня (слева) на конечный момент времени  для определенной схемы, определенного начального условия , поля скоростей ,  и заданного . Также указано максимальное и минимальное значение численного решения на конечный момент времени. Из всей пространственной области для рисования вырезана подобласть, содержащая ненулевые значения численного решения. Её начальные и конечные координаты можно видеть на рисунке.

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

Таблица норм точного решения задачи.

1.     Монотонизованная схема К. И. Бабенко.

Численное решение, полученное по предложенной схеме на конечный момент времени, «размазывается» на небольшое число узлов сетки. Наилучшим образом сохраняется форма начального профиля для третьего начального условия. Амплитуда численного решения практически неизменна. Для двух других начальных условий характерно некоторое понижение амплитуды численного решения.

Рис. 1. Монотонизованная схема К. И. Бабенко. Решение задачи с начальным профилем «конус» на момент времени = 120. Начальный профиль переместился на 9 своих длин по диагонали квадрата без вращения. .

                                    Ошибка численного решения.

Рис. 2. Монотонизованная схема К. И. Бабенко. Решение задачи с начальным профилем «конус без сектора» на момент времени = 120. Начальный профиль переместился на 9 своих длин по диагонали квадрата без вращения. .

                                    Ошибка численного решения.

Рис. 3. Монотонизованная схема К. И. Бабенко. Решение задачи с начальным профилем «М» на момент времени = 120. Начальный профиль переместился на 9 своих длин по диагонали квадрата без вращения. .

                                    Ошибка численного решения.

Рис. 1.1. Монотонизованная схема К. И. Бабенко. Решение задачи с начальным профилем «конус» на момент времени  = 480. Начальный профиль переместился на 17 своих длин (траектория движения в приложении 2) без вращения. .

                                    Ошибка численного решения.

Рис. 2.1. Монотонизованная схема К. И. Бабенко. Решение задачи с начальным профилем «конус без сектора» на момент времени = 480. Начальный профиль переместился на 17 своих длин (траектория движения в приложении 2) без вращения. .

                                    Ошибка численного решения.

Рис. 3.1. Монотонизованная схема К. И. Бабенко. Решение задачи с начальным профилем «М» на момент времени = 480. Начальный профиль переместился на 17 своих длин (траектория движения в приложении 2) без вращения. .

                                    Ошибка численного решения.

2.     Схема с «лимитерами» [8].

Для численного решения, полученного по данной схеме, характерно некоторое понижение амплитуды, более сильно выраженное, чем в предыдущем случае. Линии уровня начальных профилей «конус» и «конус без сектора» приобретают более ярко выраженную «квадратную» форму, по сравнению с формой в начале расчета. Наилучшим образом передано численное решение задачи с начальным профилем «М».

Рис. 4. Схема с «лимитерами». Решение задачи с начальным профилем «конус» на момент времени = 120. Начальный профиль переместился на 9 своих длин по диагонали квадрата без вращения. .

                                    Ошибка численного решения.

Рис. 5. Схема с «лимитерами». Решение задачи с начальным профилем «конус без сектора» на момент времени = 120. Начальный профиль переместился на 9 своих длин по диагонали квадрата без вращения. .

                                    Ошибка численного решения.

Рис. 6. Схема с «лимитерами». Решение задачи с начальным профилем «М» на момент времени = 120. Начальный профиль переместился на 9 своих длин по диагонали квадрата без вращения. .

                                    Ошибка численного решения.

Рис. 6.1. Схема с «лимитерами». Решение задачи с начальным профилем «М» на момент времени = 480. Начальный профиль переместился на 17 своих длин (траектория движения в приложении 2) без вращения. .

                                    Ошибка численного решения.

3.  Схема с направленными разностями.

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

Рис. 7. Схема с направленными разностями. Решение задачи с начальным профилем «М» на момент времени = 120. Начальный профиль переместился на 9 своих длин по диагонали квадрата без вращения. .

                                    Ошибка численного решения.

4. Схема П. Лакса [9].

          Качество численного решения в этом случае близко по своим характеристикам к качеству решения, полученного по схеме с направленными разностями (начальный профиль «размазывается», амплитуда решения уменьшается). Вдобавок к негативным характеристикам появляются характерные коротковолновые осцилляции. Сравнение норм численных решений, представленных в таблицах к рис.7 и рис.8, показывает, что метод Лакса дает худший результат по всем показателям.

Рис. 8. Схема П. Лакса. Решение задачи с начальным профилем «М» на момент времени = 120. Начальный профиль переместился на 9 своих длин по диагонали квадрата без вращения. .

                                    Ошибка численного решения.

Заключение.

1.     Предложенная монотонизованная схема К. И. Бабенко («квадрат») показала наилучшие результаты среди всех рассмотренных в работе схем (при движении начального профиля без вращения),  в особенности, в случае негладкого начального профиля «М».

2.     Численные решения, полученные по схеме с направленными разностями и схеме П. Лакса, существенно уступают по точности решениям, полученным с помощью новой монотонизованной схемы К. И. Бабенко и схемы с «лимитерами».

Приложение 1.

1)     Схема с направленными разностями.

где

2)     Схема Лакса. [9]

где .

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

 , .

3)     Схема с «лимитерами». [8]

где  , ,

,

функция  определяется следующим образом:

Приложение 2.

Траектория движения начального профиля, по которой он проходит 17 своих длин:

Литература.

1.     Галанин М. П., Еленина Т. Г. Сравнительный анализ разностных схем для линейного уравнения переноса. // Препринт ИПМ им. М. В. Келдыша РАН, № 52, 1998 г., 33 с.

2.     Галанин М. П., Еленина Т. Г. Тестирование разностных схем для линейного уравнения переноса. // Препринт ИПМ  им. М. В. Келдыша РАН, № 40, 1999 г., 42 с.

3.     Галанин М. П., Еленина Т. Г. Нелинейная монотонизация разностных схем для линейного уравнения переноса. // Препринт ИПМ им. М. В. Келдыша РАН, № 44, 1999 г., 30 с.

4.     А. А. Самарский. Теория разностных схем. М.: Наука, 1977, 656 с.

5.     К. И. Бабенко, Г. П. Воскресенский. Численный метод расчета пространственного обтекания тел сверхзвуковым потоком газа. // ЖВМиМФ, 1961, т. 1, № 6, с. 1051 - 1060.

6.     С. К. Годунов. Уравнения математической физики. М. : Наука, Физматлит, 1979, 392 с.

7.     Еленина Т. Г. Решение нелинейной монотонизированной схемы К.И.Бабенко («квадрат»). // Препринт ИПМ им. М. В. Келдыша РАН, 2001 г., в печати.

8.     К. В. Вязников, В. Ф. Тишкин, А. П. Фаворский, М. Ю. Шашков. Квазимонотонные разностные схемы повышенного порядка точности. // М.: Препринт ИПМ АН СССР, 1987, №36, 27 с.

9.     Ю. Н. Днестровский, Д. П. Костомаров. Математическое моделирование плазмы. М.: Наука, 1982, 319 с.