Решение уравнения Бюргерса с использованием вейвлет-базисов

We present, numerical simulation of Burgers’ equation using the combination of standard Galerkin techniques and multiscale compactly supported orthonormal wavelet bases of Daubechies. The numerical investigations for different values of viscosity and for various 'initial data have been performed. Daubechies" basis is found to stably represent the solution for small, and even vanishing viscosity.

Уравнение Бюргерса [2|

UiD = vu_n:, Цх,0) — uo(x), (1)

представляет собой простейшее волновое уравнение, сочетающее диссипативные и нелинейные эффекты, и появляется поэтому во многих физических приложениях. При нулевой вязкости и за время t* = mni(l/uo,$) образуются разрывы (скачки), и решение приобретает характерный “пилообразный⁷ профиль ударных волн. Для малой вязкости решение является сглаженной версией невязкого случая, т.с. возникают сильные градиенты, которые затем медленно диссипируют, и при t -п ос решение вырождается до нуля. Несмотря на. свою простоту, уравнение обладает достаточно богатой динамикой, связанной с нелинейным взаимодействием мод различных масштабов, что выражается в образовании упорядоченных крупномасштабных структур, “схлопывании” мелкомасштабных возмущений и т.п. Решения уравнения Бюргерса со случайными начальными данными показывают два. механизма, которые свойственны реальной турбулентности: нелинейный перенос энергии по спектру и вязкая диссипация в области мелких масштабов. Все это позволяет предполагать, что подходы, успешно применяемые для разрешения уравнения Бюргерса, могут быть эффективными и для моделирования реальной гидродинамической турбулентности.

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

Уравнение Бюргерса является хорошим тестовым примером для численных методов вследствие его простоты и предсказуемости динамики. Более того, используя преобразование Хопфа-Коула, можно получить точное решение, что позволяет оценить погрешность различных численных методов, сравнивая их с точным решением.

Вейвлеты

Напомним, что вейвлет-функции локализованы как в реальном пространстве, так и в пространстве частот, и получаются друг из друга, путем масштабного преобразования (растяжения/сжатия) и сдвига (см., напрмер. [4. 9]). Кроме того, вейвлет-функция ф обязана удовлетворять так называемому условию допустимости (admissibility condition)

J ф(х) dx — 0.

Некоторые классы вейвлет-функций образуют полные ортонормированные базисы пространства. I? на прямой или окружности. При этом все функции, образующие' ортонормированный вейвлет-базис, должны быть самоподобны и получаются переносом и растяжением/сжатием только одной “материнской функции, т.е.

Нт,» - 2^т/2П(2’М - пф т, п е Z.

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

/ V           ) Ф^ГП12 ,«2 (Д ) Фа Аг; ‘

Метод Галеркина

Представим решение гЦх, /) в виде ряда, по ортонормированной системе вейвлет-функций

АхЛ) = An^Wx) = ANAty2N^ - nV где A.yn(t) есть амплитуда n-й базисной функции JV-го уровня (масштаба). Подставляем теперь решение в исходное уравнение Бюргерса (1) и, используя стандартную технику метода. Галеркина, проектируем уравнение на пространство вейвлет-функций и получаем систему обыкновенных дифференциальных уравнений относительно коэффициентов Ап;,ФА)

— У/ RNnMm,LtAM_mAu + v У/                      (2)

где массив

Кк.иМтДЛ " — / Л\пфаАФмтАА--,--da                  рЗ I dx описывает нелинейные взаимодействия, обусловленные конвективной производной. и матрица

AN.i,Mm = / WnW --ГМ-- их                   (4)

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

Конечно, остается определенный произвол в выборе диапазонов изменения индексов А⁷ и п. Очевидно, что максимальное значение индекса. Л⁷ (минимальный масштаб) определяет пространственно!' разрешение метода; максимальный же масштаб можно выбирать, исходя из конкретной постановки задачи (например, начальных условий). Диапазон изменения индекса //. тесно связан с граничными условиями задачи. В случае периодического решения число функций данного уровня определяется из условий полноты базиса: функции путем сдвига должны покрывать всю окружность (очевидно, на каждом следующем более мелком масштабе число базисных функций удваивается). Для непериодического локализованного в пространстве решения в силу того, что система вейвлет-функций не образует полный базис для конечного интервала, неминуемо возникает проблема, краевых эффектов и. следовательно, точно удовлетворить граничные условия невозможно (хотя, в принципе, можно построить полный ортонормированный вейвлет-базис для отрезка). Тем не менее, для функций с компактным носителем краевые эффекты будут локализованы в окрестностях концов рассматриваемого интервала и в большинстве случаев можно пренебречь их влиянием .

Вейвлеты Добеши

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

Вычисление элементов матриц R и К

Остановимся подробно на важном для приложений вопросе, связанном с вычислением элементов матриц R и К, т.е. с вычислением интегралов от произведений самих вейвлет-функций и их производных. Вейвлеты с компактным носителем имеют конечное число производных, которые, особенно высших порядков, являются сильно осциллирующими и нерегулярными. Это делает любые¹ численные оценки интегралов весьма трудоемкими и неустойчивыми. Тем не менее, для вейвлетов Добеши существует точный алгоритм получения значений интегралов вида (3). (4) и им подобных

Этот алгоритм существенным образом опирается на дополнительные соотношения, даваемые так называемым анализом с переменным разрешением (multiresolntion analysis) [8], которому подчиняется подавляющее большинство вейвлет-функции. образующих ортонормированные базисы. Из условий, налагаемых анализом с переменным разрешением, вытекает существование помимо “материнской" вейвлет-функции -ф и так называемой масштабной функции (scaling function) ф. При этом о удовлетворяет следующим условиям:

I ф(х) dx = 1,                                 (5)

I Ффх — пфффх — m) dx = Ф„„, п, т € Z.

Кроме того, масштабная функция удовлетворяет масштабному соотношению ффхф " 52 Ь.кф(2х - кф                            (6)

ку.

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

= 52 9кФ(2х - кф где д_к =            .                (7)

fcez

Очевидно, что в случае функций с компактным носителем отлично от нуля лишь конечное число коэффициентов h_k и д_к. А именно, если отличны от нуля А значений коэффициентов фильтра. h_k (скажем для к = 0,. . . , N — 1). то соответствующие масштабная и вейвлет-функции отличны от нуля на отрезке [О, А — 1]. На рис. 1 приведены масштабная и вейвлет-функции для шеститочечного фильтра к.

Нетрудно видеть, что соотношения (6) и (7) дают возможность вычислить коэффициенты Ry_nMm,u и KNn,Mm Для любых значений индексов, зная лишь интегралы от произведения соответствующих масштабных функций и их производных:

причем для функций Добеши лишь конечное их число отлично от нуля. Используя масштабное соотношение (6) и .замену переменных интегрирования, можно получить систему линейных однородных алгебраических уравнений относительно искомых значений интегралов. Так. например, для интегралов вида (8) и шеститочечной функции Добеши (DG) система уравнений будет иметь вид

К - 2 J' hJi₃Ivj_TlHi_₃ = 0,     -4 <  т < 4.                  (9)

где. очевидно. К_т = 0. если т < —4 или т > 4. Получающиеся таким образом системы будут сингулярными.

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

52 ф(х — k) = 1,     52 кф(х — к) = х + .4,

52 FoU x2 + Bx + С, кх7, где А, В, С - некоторые константы. Используя последнее тождество и равенство (5 ). можно замкнуть сис,тему уравнений (9)

IF - 2 £           = 0,    -4 < т < 3,

*    4                                                                  ДО)

к —— 4

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

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

Численные результаты

Для численного моделирования были использованы шеститочечные вейвлеты Добеши (D6)- функции с минимально возможной гладкостью (и тем самым с минимальной шириной носителя) для данной задачи. Для интегрирования системы (2) применялся стандартный метод Рунге-Кутта с постоянным шагом по времени. В качестве тестовой задачи была рассмотрена эволюция начального периодического поля, представляющего собой синусоиду, для различных значений вязкости (рис. 2 5). Вейвлет-базис строился периодическим образом. Для численного счета мы использовали 8 масштабов, что с учетом периодичности дает 255 базисных функций на. период (128 функций самого мелкого масштаба). Для рис. 2 4 ширина полосы взаимодействующих масштабов оставалась максимально возможной (8 масштабов), т.е. все масштабы взаимодействуют друг с другом.

На. рис. 2 приведены результаты численного счета для экстремально малой вязкости и = 0.0001 (соответствующее число Рейнольдса Re to 6 • 10°). Возмущение решения вследствие феномена Гиббса имеет место, но метод остается устойчивым и на больших временах решение стремится к точному.

Для несколько большей вязкости (рис. 3) и = 0.001, Re to 6 ■ 10;’ отклонение от точного решения еще заметно, но оно локализованно в окрестности скачка, и не вызывает осцилляций, характерных для спектральных методов.

При увеличении вязкости (рис. 4) и = 0.01, Re to 6 • 10⁴ метод показывает полностью корректное поведение решения.

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

На рис. 6 9 продемонстрировано сохранение на больших временах крупномасштабных структур, присутствующих в начальных условиях [1].

Рис. 6 показывает начальное поле вида sin 4т+ 0.1 sinx (и = 0.001, Re to 6-10°), а рис. 7 демонстрирует, что на достаточно больших временах происходит исчезновение мелкомасштабных возмущений, и доминируют структуры, обусловленные’ крупномасштабной составляющей (0.1 sin т) начального поля.

Наконец, на рис. 8,9 показано аналогичное явление для начального случайного поля. На рис. 8 представлено начальное поле, возмущениям подвергались два. масштаба (1-й и 5-й). Рис. 9 демонстрирует сохранение с течением времени только достаточно крупномасштабных мод.

Рис. 4                                                  Рис.

Рис. 8                                                     Рис.

Выводы

Метод Галеркина с использованием вейвлет-базисов Дебоши остается устойчивым для любой вязкости, и решение сходится к точному при исчезающе малой вязко сти (тля Rc 10^s - 1(F).

Осцилляции в окрестности сильных градиентов, связанные с явлением Гиббса, незначительны и локализованны в непосредственной близости от скачков. Даже' для достаточно малой вязкости (Re ^ 104 — 1(Г) возмущения, вызываемые явлением Гиббса, пренебрежимо малы, и вейвлет-аппроксимации корректно описывают как гладкие компоненты решешта, так и участки сильных градиентов.

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

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

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

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант 94-0Т00951-а). Численные расчеты стали возможны благодаря грантам Шведской Королевской Академии наук, фонда Ерпкссона (Швеция) и материальной поддержке Департамента численного анализа, и информатики Короле вского технологического института. (Стокгольм).