Конечно разностная схема для волнового уравнения

Волновое уравнение

Одним из наиболее распространенных в инженерной практике уравнений с частными производными второго порядка является волновое уравнение, описывающее различные виды колебаний. Поскольку колебания — процесс нестационарный, то одной из независимых переменных является время t. Кроме того, независимыми переменными в уравнении являются также пространственные координаты х, у, z. В зависимости от их количества различают одномерное, двумерное и трехмерное волновые уравнения.

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

Рассмотрим одномерное волновое уравнение, которое можно записать в виде

(2.63)

Для поперечных колебаний струны искомая функция U(x,t) описывает положение струны в момент t. В этом случае а2 = Т/ρ, где Т — натяжение струны, ρ — ее линейная (погонная) плотность. Колебания предполагаются малыми, т.е. амплитуда мала по сравнению с длиной струны. Кроме того, уравнение (2.63) записано для случая свободных колебаний. В случае вынужденных колебаний в правой части уравнения добавляют некоторую функцию f(x,t), характеризующую внешние воздействия, при этом сопротивление среды колебательному процессу не учитывается.

Простейшей задачей для уравнения (2.63) является задача Коши: в начальный момент времени задаются два условия (количество условий равно порядку входящей в уравнение производной по t):

(2.64)

Эти условия описывают начальную форму струны и скорость ее точек .

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

(2.65)

Рассмотрим некоторые разностные схемы для решения задачи (2.63)-(2.65). Простейшей является явная трехслойная схема типа крест (шаблон показан на рис. 2.21). Заменим в уравнении (2.63) вторые производные искомой функции Uпо tи х их конечно-разностными соотношениями с помощью значений сеточной функции в узлах сетки :

Рис. 2.21. Шаблон явной схемы

Отсюда можно найти явное выражение для значения сеточной функции на ( j + 1)-ом слое:

(2.66)

Здесь, как обычно в трехслойных схемах, для определения неизвестных значений на (j + 1)-ом слое нужно знать решения на j-ом и (j — 1)-ом слоях. Поэтому начать счет по формулам (2.66) можно лишь для второго слоя, а решения на нулевом и первом слоях должны быть известны. Их находят с помощью начальных условий (2.64). На нулевом слое имеем

(2.67)

Для получения решения на первом слое воспользуемся вторым начальным условием (2.64). Производную заменим конечно-разностной аппроксимацией. В простейшем случае полагают

(2.68)

Из этого соотношения можно найти значения сеточной функции на первом временном слое:

(2.69)

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

(2.70)

Вместо нужно взять . А выражение для второй производной можно найти с использованием исходного уравнения (2.63) и первого начального условия (2.64). Получим

Тогда (2.70) примет вид:

(2.71)

Разностная схема (2.66) с учетом (2.71) обладает погрешностью аппроксимации порядка

При решении смешанной задачи с граничными условиями вида (2.65), т.е. когда на концах рассматриваемого отрезка заданы значения самой функции, второй порядок аппроксимации сохраняется. В этом случае для удобства крайние узлы сетки располагают в граничных точках (х0=0, xI = l). Однако граничные условия могут задаваться и для производной.

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

(2.72)

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

Рассмотренная разностная схема (2.66) решения задачи (2.63) — (2.65) условно устойчива. Необходимое и достаточное условие устойчивости:

(2.73)

Следовательно, при выполнении этого условия и с учетом аппроксимации схема (2.66) сходится к исходной задаче со скоростью O(h2+τ2). Данная схема часто используется в практи-ческих расчетах. Она обеспечивает приемлемую точность получения решения U(x,t), которое имеет непрерывные производные четвертого порядка.

Рис. 2.22. Алгоритм решения волнового уравнения

Алгоритм решения задачи (2.63)-(2.65) с помощью данной явной разностной схемы приведен на рис. 2.22. Здесь представлен простейший вариант, когда все значения сеточной функции, образующие двумерный массив, по мере вычисления хранятся в памяти компьютера, а после решения задачи выводятся результаты. Можно было бы предусмотреть хранение решения лишь на трех слоях, что сэкономило бы память. Результаты в таком случае можно выводить в процессе счета (см. рис. 2.13).

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

Построим простейшую неявную схему. Вторую производную по tв уравнении (2.63) аппроксимируем, как и ранее, по трехточечному шаблону с помощью значений сеточной функции на слоях j— 1, j, j + 1. Производную до х заменяем полусуммой ее аппроксимации на (j + 1)-ом и (j— 1)-ом слоях (рис. 2.23):

Рис. 2.23. Шаблон неявной схемы

Из этого соотношения можно получить систему уравнений относительно неизвестных значений сеточной функции на (j+ 1)-ом слое:

(2.74)

Полученная неявная схема устойчива и сходится со скоростью . Систему линейных алгебраических уравнений (2.74) можно, в частности, решать методом прогонки. К этой системе следует добавить разностные начальные и граничные условия. Так, выражения (2.67), (2.69) или (2.71) могут быть использованы для вычисления значений сеточной функции на нулевом и первом слоях по времени.

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

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

Конечно разностная схема для волнового уравнения

В последние 10 лет наблюдается повышенный интерес к исследованию диффузионных процессов, не подчиняющихся закону Фика и не описываемых классическим уравнением диффузии [1,2]. Такие процессы наблюдаются на практике, в частности, в высокоэнергетической плазме, при переносе во фрактальных средах и в аморфных полупроводниках. Оказалось, что для математического описания таких процессов удобно использовать аппарат дробного интегродифференцирования. Однако, наличие в уравнении дробных производных приводит к существенным трудностям при построении как аналитических, так и численных решений. В связи с этим является актуальной задача разработки эффективных алгоритмов численного решения краевых задач для дифференциальных уравнений в частных производных дробного порядка.

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

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

производная Римана — Лиувилля, и . Уравнение (1.1) дополняется граничными условиями первого рода (условия Дирихле)

и начальными значениями

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

Для аппроксимации дробной производной Римана-Лиувилля воспользуемся формулой Грюнвалда [9] и L1-аппроксимацией [7]. В результате в первом случае получаем метод конечных разностей (МКР) для дробных производных, а во-втором случае приходим к L1-схеме.

Если непрерывна, а интегрируема на отрезке , тогда для любого порядка производные Римана — Лиувилля и Грюнвалда — Летникова существуют и совпадают для любого момента времени из [2].

Определение Грюнвалда — Летникова позволяет численно находить производную Римана — Лиувилля:

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

Для справедливы рекуррентные соотношения:

Важно отметить, что аппроксимация (1.6) справедлива, либо при
, либо когда достаточна гладкая в окрестности [8].

Для простоты рассмотрим случай const . Уравнение (1.1), согласно (1.4) и (1.6), заменяется явным конечно-разностным уравнением :

Следуя [9], отметим, что при указанных коэффициентах схема будет иметь первый порядок точности.

Утверждение 1. Для устойчивости схемы (1.8) достаточно выполнения неравенства

Д о к а з а т е л ь с т в о. Для анализа устойчивости разностной схемы воспользуемся методом фон Неймана. Предположим, что решение имеет вид , где — вещественное пространственное волновое число. Подставим это представление в выражение (1.8) и полученное выражение разделим на . C учетом того факта, что