Лабораторная работа: Симуляция водной поверхности.

Задания

Рассмотрим волновое уравнение:

\[\frac{\partial^2\eta}{\partial t^2}=v^2\Delta\eta, \text{ где } \Delta\eta=\frac{\partial^2\eta}{\partial x^2}+\frac{\partial^2\eta}{\partial y^2},\]

$x,y$ обозначают декартовы координаты на плоскости дна, $\eta(x,y)$ равна высоте столба жидкости над точкой дна $(x,y)$, $v$ задает скорость распространения волн. Из каких соображений выводится это уравнение? При каких предположениях волновое уравнение моделирует водную поверхность? Будет ли буй на поверхности воды двигаться со скоростью $v$? Почему уравнение называется волновым? Какие еще системы описываются волновым уравнением?

Если рассматриваемый объем воды ограничен не только по вертикали, то на границе поверхности нужно задать граничные условия, выражающие закон взаимодействия с соприкасающейся средой. Если рассматриваемая область имеет вид $-a\leq x\leq a$, $-a\leq y\leq a$, то в качестве граничных условий можно взять условия Неймана:

\[\frac{\partial\eta}{\partial x}\bigg|_{x=\pm a}=0,\quad \frac{\partial\eta}{\partial y}\bigg|_{y=\pm a}=0.\]

Что означают условия Неймана в физических терминах? Насколько адекватно отражают условия Неймана поведение воды на границе? Можно ли обойтись вообще без граничных условий? Какая физическая система может быть описана волновым уравнением с условиями Дирихле: $\eta=0$ во всех точках границы?

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

\[x_k:=-a+2ak/N,\; y_j:=-a+2aj/N,\text{ где }k,j=0\ldots N.\]

В дискретной модели функция $\eta$ задается значениями в узлах решетки: $\eta_{kj}=\eta(x_k,y_j)$. Тогда значение лапласиана может быть приближенно выражено через вторые конечные разности:

\[\Delta\eta_{kj}:=\Delta\eta(x_k,y_j)\approx \frac1{\delta x^2}(\eta_{k+1,j}+\eta_{k-1,j}+\eta_{k,j+1}+\eta_{k,j-1} -4\eta_{kj}),\]

где $\delta x=2a/N$ обозначает шаг решетки.

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

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

\[\frac{\partial\eta(x,y)}{\partial x}\approx \frac{\eta(x+\delta x/2,y)-\eta(x-\delta x/2,y)}{\delta x}.\]

Насколько точна эта оценка? Как изменится формула для лапласиана и ее точности, если вместо симметричной конечной разности использовать прямую конечную разность:

\[\frac{\partial\eta(x,y)}{\partial x}\approx \frac{\eta(x+\delta x,y)-\eta(x,y)}{\delta x}?\]

Для оценки точности используйте формулу Тейлора:

\[f(x+\delta)=f(x)+f'(x)\delta+f''(x)\frac{\delta^2}{2} +f'''(x)\frac{\delta^3}{3!}+\ldots\]

Получите более точную формулу для $\Delta\eta$, используя значения в большем числе соседних узлов.

Рассмотрим теперь частный случай волны, высота которой не зависит от координаты $y$, тогда функция $\eta=\eta(x)$ зависит только от $t$ и $x$, следовательно волновое уравнение принимает вид:

\[\frac{\partial^2\eta}{\partial t^2} =v^2\frac{\partial^2\eta}{\partial x^2},\]

и сохраняется только половина граничных условий:

\[\frac{\partial\eta}{\partial x}\bigg|_{x=\pm a}=0.\]

Покажите, что в этом случае волновое уравнение имеет решения в виде стоячих волн:

\[\eta(x;t)=\cos\frac{\pi n x}{a}\cos\frac{v\pi n t}{a},\]

и в виде бегущих волн:

\[\eta(x;t)=\phi(x\pm vt),\]

где $\phi$ - любая гладкая функция, такая что $\eta$ принимает непостоянные значения лишь на некотором удалении от границы.

Пусть дана некоторая дискретизация волнового уравнения по переменной $x$, тогда

\[\ddot\eta(t)=v^2L\eta(t),\]

где $\eta(t)$ представляет собой вектор конечного размера (например, значения в узлах сетки), через $L\eta$ записана дискретизация лапласиана $\Delta\eta$, а точки обозначают производные по времени, $\ddot\eta=\partial^2\eta/\partial t^2$. Эволюция системы становится однозначной, если задать начальные условия вида:

\[\eta\bigg|_{t=0}=\psi^0,\quad \dot\eta\bigg|_{t=0}=\psi^1,\]

где $\psi^0$ - высота водной поверхности над нулевым уровнем в начальный момент времени, а $\psi^1$ - вертикальная компоненты скорости точки поверхности в начальный момент времени (не путать с $v$). Какие начальные условия соответствуют решениям в виде стоячих и бегущих волн?

Для численного решения волнового уравнения его нужно дискретизовать, что можно сделать, например, задавшись постоянным шагом по времени $\delta t$, $\eta^q=\eta|_{t=q\delta t}$, и заменив частную производную по времени центральной конечной разностью:

\[\frac{\eta^{q+1}-2\eta^q+\eta^{q-1}}{\delta t^2} =v^2 L\eta^t.\]

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

Выразив $\eta^{q+1}$ через значения $\eta$ в предыдущие моменты времени, и приближенно вычислив

\[\eta^1=\eta^0+\delta t\cdot\dot\eta^0 =\psi^0+\delta t\cdot \psi^1,\]

найдите численно решения с начальными условиями для стоячих волн и для бегущей волны. Вычислите, насколько отличается $\eta^2$ от аналитического решения? Постройте график зависимости ошибки вычисления $\eta^2$ в зависимости от $\delta t$. Как по графику найти порядок точности численного интегрирования.

Найдите $\eta^{100}$ и, сравнив найденное численно значение с аналитическим, постройте график ошибки интегрирования как функции от $\delta t$. Как видно на графике, при значениях шага интегрирования $\delta t$ больше некоторого критического, погрешность драматически растет. Объясните причину неустойчивости метода интегрирования при больших шага $\delta t$. Как можно найти критическое значения шага аналитически? Как зависит это критическое значение от шага сетки в пространстве $\delta x$ и от скорости волн $v$?
Энергия поверхности воды имеет вид:

\[E[\eta,\dot\eta]=\int_{\Omega}\bigg[ \bigg(\frac{\partial\eta}{\partial t}\bigg)^2+ v^2\bigg(\frac{\partial\eta}{\partial x}\bigg)^2+ v^2\bigg(\frac{\partial\eta}{\partial y}\bigg)^2 \bigg]dx\,dy.\]

где интеграл берется по всему объему воды $\Omega$. Покажите, что волновое уравнение выражает закон сохранения энергии:

\[\frac{d E}{dt}=0.\]

Вычислите некоторую дискретизацию энергии конечными разностями, и проверьте, насколько сохраняется энергия при численном интегрировании. Почему сохранение энергии важно? Как зависит точность сохранения энергии от шага интегрирования по времени $\delta t$? Можно ли предложить альтернативную дискретизацию волнового уравнения, чтобы энергия сохранялась точно?

Рассмотрим волновое уравнение с периодическими граничными условиями

\[\eta\bigg|_{x=-a}=\eta\bigg|_{x=a},\quad \frac{\partial\eta}{\partial x}\bigg|_{x=-a}=\frac{\partial\eta}{\partial x}\bigg|_{x=a}.\]

Решения в виде стоячей и бегущей волн сохраняет свою форму на любых промежутках времени. Проверьте, изменяется ли форма волны этих волн при численном интегрировании. Выберете форму бегущей волны в виде синусоиды, и сделав несколько шагов численного интегрирования, найдите скорость движения волны. Постройте график зависимости скорости волны от ее частоты, сравните со скоростью волны $v$ из волнового уравнения. Объясняет ли зависимость скорости волны от ее длины (дисперсия волн) изменение формы волнового пакета? Почему при численном интегрировании возникает зависимость скорости волны от частоты, хотя согласно волновому уравнению скорости всех волн одинаковы? Как аналитически вычислить зависимость скорости от частоты?

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

\[E[\eta,\dot\eta]=T[\dot\eta]+U[\eta],\]

это позволяет применить простейший симлектический интегратор: метод Верле. Реализуйте метод Верле и выполните интегрирование стоячих и бегущих волн. Оцените устойчивость метода, сохранение энергии и выраженность дисперсии волн, аналогично тому, как это было проделано выше.

Точность интегрирования можно заметно повысить, перейдя к методам интегрирования более высокого порядка. Наиболее популярным семейством интеграторов, допускающих сколь угодно высокий порядок интегрирования, является семейство методов Рунге-Кутты. Реализуйте классический метод Рунге-Кутты четвертого порядка. Проанализируйте устойчивость, сохранение энергии, дисперсию волн. Метод Верле является методом второго порядка, означает ли это, что метод Рунге-Кутты четвертого порядка предпочтительнее метода Верле? Что означает порядок метода? В каких случаях дополнительные затраты на вычисления в методе более высокого порядка оправданы?
Среди методов Рунге-Кутты выделяют подсемейство симплектических интеграторов, сохраняющих квадратичные формы и энергию в частности. Эти симплектические интеграторы относятся к неявным методам Рунге-Кутты. Что означает “неявный метод”? Как сделать шаг интегрирования неявным методом, если решение в следующий момент времени трудно явно выразить через предыдущие моменты времени? Если решение неявным методом предполагает значительно больший объем вычислений, чем явные методы, почему неявные методы вообще используются? Реализуйте метод Radau, проведите анализ устойчивости, сохранения энергии и наличия дисперсии волн.
Реализуйте один из вышеизученных методов для интегрирования двухмерного волнового уравнения. Объедините интегратор с системой визуализации, разработанной в предыдущей лабораторной работе. Как переход к двухмерной задаче повлиял на устойчивость? На изменение энергии? На скорость волн?
Изучите уравнение мелкой воды. Как оно выводится? При каких условиях оно верно описывает движение воды? Как оно соотносится с волновым уравнением?
Проведите дискретизацию уравнения мелкой воды методом конечных разностей. Уравнение мелкой воды обычно записывается с производными от произведений. Почему производные не посчитаны до конца? Как лучше: заменить производные от произведений на конечные разности, или сначала раскрыть производную произведения, а затем заменить производные на конечные разности? На каких решетках будут заданы входящие в уравнение мелкой жидкости величины, если при дискретизации использовать центральные конечные разности? Можно ли совместить узлы решеток для разных компонент скорости? Как отобразить поле скоростей, если проекции скорости на разные оси известны в разных точках?
Реализуйте какой-либо метод численного интегрирования для уравнения мелкой воды. Отобразите динамику поверхности воды, используя визуализатор из предыдущей лабораторной работы.
Часто используемой альтернативой дискретизации конечными разностями, является дискретизация конечными элементами. Почему метод конечных элементов весьма популярен, хотя он более сложен, чем метод конечных разностей? Почему для решения уравнений в частных производных в областях со сложной границей почти всегда используется метод конечных элементов? В каких случаях метод конечных разностей более предпочтителен?
Проведите дискретизацию волнового уравнения конечными элементами. Сравните точность дискретизации конечными элементами и конечными разностями.
Любую волну можно разложить в сумму плоских волн, использую преобразование Фурье. Найдите, как эволюционируют плоские волны под действием волнового уравнения. Используйте преобразование Фурье для нахождения аналитического решения волнового уравнения. Как можно использовать аналитического решение такого рода для моделирования водной поверхности?