Волновое уравнение
Одним из наиболее распространенных в инженерной практике уравнений с частными производными второго порядка является волновое уравнение, описывающее различные виды колебаний. Поскольку колебания — процесс нестационарный, то одной из независимых переменных является время 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) могут быть использованы для вычисления значений сеточной функции на нулевом и первом слоях по времени.
При двух или трех независимых пространственных переменных волновые уравнения принимают вид
Для них также могут быть построены разностные схемы по аналогии с одномерным волновым уравнением. Разница состоит в том, что нужно аппроксимировать производные по двум или трем пространственным переменным, что, естественно, усложняет алгоритм и требует значительно больших объемов памяти и времени счета. Подробнее двумерные задачи будут рассмотрены ниже для уравнения теплопроводности.
Двухслойная акустическая схема
Лекция №13
Волновое уравнение
Гиперболическое волновое уравнение описывает колебания струны, движение сжимаемого газа, распространение электромагнитных волн и ряд других явлений.
Типичной одномерной задачей является задача описания малых колебаний натянутой струны с распределенной по длине нагрузкой f(t,x):
; (1)
; (2)
. (3)
Уравнение колебания струны (1) дополняется начальными условиями (2), которые, в отличие от параболического уравнения требуют двух условий: начального смещения относительно положения равновесия u0(x) и начальную скорость движения u1(x). Кроме того, задаются краевые условия (3), которые называются краевыми условиями первого рода, они описывают смещение концов струны относительно положений равновесия. Краевые условия могут быть и иного рода.
Составим несложную, но эффективную разностную схему для численного решения задачи (1) — (3), выбирая для простоты равномерные по t и x сетки. В качестве шаблона разностной схемы возьмем, представленный на рис.1 шаблон в форме “креста”. Аппроксимируя производные в (1) конечными разностями, получим трехслойную схему следующего вида
, (4)
с граничными условиями
. (4¢)
Рис.1. Разностная схема (4) в форме креста
По форме шаблона схему (4) называют схемой “крест”. Исследуем схему крест.
Процедура вычисления решения выглядит следующим образом. На нулевом слое решение известно из начального условия
. (5)
На первом слое решение можно также вычислить по начальным данным. Простейший способ получения решения на первом слое состоит в аппроксимации второго начального условия в (2) согласно представлению:
. (6)
Аппроксимация (6) имеет первый порядок точности, тогда как в схеме крест разностный оператор второй производной по времени имеет второй порядок точности. Поэтому для более точной аппроксимации необходимо учесть еще один член разложения, т.е.
. (6¢)
Подставляя в (6¢) utt из (1), найдем
. (6¢¢)
Схема крест (4) является явной, поскольку позволяет выразить решение на следующем слое 
через значения yn и 
на двух предыдущих слоях. Поэтому, зная значения решения на нулевом (5) и первом слое (6¢¢), можно вычислить решения на всех последующих слоях. Итак, разностное решение существует и единственно.
Для изучения аппроксимации схемы крест, разложим точное решение по формуле Тейлора с центром в узле (tm,xn), считая все появляющиеся в разложении производные непрерывными:
Используя данные разложения, легко можно найти невязку схемы крест:
(7)
а также невязку начального условия (6):
,
или более точного начального условия (6¢¢):
(7¢)
Начальные (2) и краевые условия (3) аппроксимируются точно. В итоге разностная схема (4) с начальными условиями (5), (6¢¢) имеет порядок аппроксимации 
, а та же разностная схема с начальным условием (6¢) имеет несколько худший порядок — 
.
Устойчивость схемы крест исследуем методом разделения переменных, полагая
. (8)
Подставляя представление (8) в (4), для множителя роста гармоник находим квадратное уравнение
. (9)
По теореме Виета произведение корней квадратного уравнения (9) равно единице, т.е. 
. Условие устойчивости 
может таким образом быть выполнено при 
. Это означает, что корни уравнения (9) должны образовывать комплексно сопряженную пару. Для этого необходимо, чтобы дискриминант уравнения (9) был неположительным для всех гармоник, т.е.
. (10)
Для выполнения неравенства (10) необходимо и достаточно соблюдение условия Куранта:
ct 2 (t,x) > 0. Будем считать, что функции k(t,x), f(t,x) в (17) кусочно-непрерывны, причем разрывы неподвижны, т.е. лежат на линиях x = const. Предполагается, что на разрывах выполняется условие сопряжения [u] = [kux] = 0, т.е. решение u и поток kux являются непрерывными функциями.
Выберем по времени равномерную сетку, а по пространству специальную неравномерную сетку, у которой все точки разрыва коэффициентов являются узлами. Построим аналог наилучшей консервативной схемы (11.19) — (11.21¢), разобранной в лекции №11:
, (18)
(18¢)
Известно[1], что при сделанных предположениях и при достаточно гладких начальных и граничных условиях схема (18), (18¢) равномерно сходится со скоростью 
при выполнении условия устойчивости (16).
Проведем численный расчет по разностной схеме (18), (18¢) задачи (17), (2), (3) в прямоугольной области G(t,x) = [0,T]´[0,a] при условии, что
(19)
(20)
где k1, k2, k3, f1, f2, f3 — некоторые постоянные величины. С учетом положений разрывов в (19), (20) определим равномерную сетку по пространству xn = nh, n = 1,…,N, где N = 4l + 1, l = 1,2,… Определим также равномерную сетку по времени, т.е. tm = t m, 1,…,M. С учетом сделанных предположений рассмотрим следующую схему аппроксимации для интегралов в (18¢):
. (21)
В качестве начальных и граничных условий положим
. (22)
На листинге_№3 приведен код программы численного решения задачи (17), (19), (20), (22) с помощью наилучшей разностной схемы (18), (18¢), (21).
%Программа решения волнового уравнения (17), (19),
%(20) с помощью разностной схемы (18), (18′), (21)
global a k1 k2 k3 f1 f2 f3
%Определяем габариты области интегрирования по
%времени T и пространству a
%Определяем константы k1,k2,k3,f1,f2,f3
k1=1; k2=10; k3=1; f1=1; f2=10; f3=1;
%Определяем весовую константу
%Определяем число узлов по пространству и времени
%Определяем шаги по времени и пространству
%Определяем сетки по времени и пространству
%Используем начальные данные из (22) для
%определения численного решения на первом и
%втором слое по формуле (6»)
%Определяем левое и правое граничные условия,
%Применяем неявную наилучшую схему (18), (18′)
%Определяем коэффициенты прогонки
%Рисуем трехмерную поверхность решения в координатах
%Определяем функцию квадрата скорости звука k(t,x)
global a k1 k2 k3
if (x>=0)&(x 0.25*a)&(x =0.75*a)&(x =0)&(x 0.25*a)&(x =0.75*a)&(x 0), инвариант s — уравнению переноса влево. Для однородной задачи, когда F = 0, m1 = m2 =0, величины r, s переносятся по своим характеристикам без изменения, с чем и связано их наименование.
Для инвариантов можно составить разностные схемы, аналогичные разностным схемам для уравнения переноса. При этом шаблон каждой из схем должен учитывать направление характеристики данного уравнения. Рассмотрим простейшую явную разностную схему для уравнений (39):
. (42)
Схема (42) является схемой бегущего счета. Нетрудно показать, что при выполнении условия Куранта ct £ h схема устойчива и сходится с порядком точности 
на дважды непрерывно дифференцируемых функциях.
Изучим разностную схему для инвариантов на примере численного решения задачи (39) с помощью разностной схемы (42), выбирая правую часть, начальные и граничные условия согласно (35), (36). В этом случае при учете (40), (41), находим следующие выражения для правой части, начальных и граничных условий:
(43)
Нетрудно проверить, что аналитическим решением задачи (39), (43) являются выражения вида:
. (44)
На листинге_№5 приведен код программы расчета задачи акустики в инвариантах (39) со специальной правой частью, начальными и граничными условиями (43).
%Программа решения уравнений акустики (39) в
%инвариантах с правой частью, начальными и
%граничными условиями вида (43)
%Определяем габариты области интегрирования по
%времени T и пространству a, а также определяем
%величину скорости c
%Определяем число удвоений числа узлов сетки по
%времени и пространству dmax
%Организуем цикл расчетов на различных сетках
%Определяем число узлов сетки по пространству
%определяем шаги по времени и пространству
%Определяем сетки по времени и пространству
%Определяем начальные условия (43) для
%инвариантов r и s
%Организуем цикл расчетов по времени
%Осуществляем бегущий расчет инварианта r,
%осуществляемый слева направо
%Учитываем правое граничное условие (43)
%Осуществляем бегущий расчет инварианта s,
%осуществляемый справа налево
%Учитываем левое граничное условие (43)
%Ошибку численного расчета инвариантов r и s
%оцениваем как разность численного и
%аналитического решений в норме C. Делим
%полученную ошибку на шаг по пространству h и
%находим предстепенную константу скорости
%сходимости схемы с инвариантами
%Рисуем численное решение инварианта r
%Рисуем численное решение инварианта s
%Рисуем график зависимости предстепенной
%константы от шага сетки h
%Определяем функцию правой части
%Определяем аналитический вид инварианта r
%Определяем аналитический вид инварианта s
Рис.8. Решение задачи акустики (39), (43) в инвариантах
На рис.8 приведен итог работы кода программы листинга_№5. На левом рисунке приведено изображение численного решения инварианта r, на среднем рисунке — инварианта s. Правый график демонстрирует зависимость const от h в представлении для ошибки численного решения следующего вида:
,
где r(t,x), s(t,x) — аналитические решения (44). Из графика видно, что по мере уменьшения шага сетки h величина const действительно выходит на некоторое постоянное значение, при этом шаг по времени выбирался порядка шага по пространству, т.е. t
h. Это подтверждает также то, что скорость сходимости разностной схемы (42) — .
Многомерные схемы
Рассмотрим волновое уравнение в p-мерном пространстве следующего вида:
. (45)
Решение задачи (45) предполагает определение начальных и краевых условий:
(46)
Схема “крест” строится аналогично одномерной схеме (4) на шаблоне, вид которого для случая двух измерений показан на рис.9. При произвольном числе измерений разностная схема крест может быть записана в виде:
, (47)
при этом в случае переменных коэффициентов ka операторы La определяются также как и для наилучшей схемы (18), (18¢).
Рис.9. Шаблон схемы крест в двумерном случае
Трехслойная схема (47) является явной. Вычисления с помощью схемы (47) одинаково просты для любого числа измерений. Легко проверить, что схема имеет порядок аппроксимации 
.
Исследуем устойчивость схемы (47) методом разделения переменных, считая коэффициенты ka постоянными. Для этого подставим в (47) многомерную гармонику, имеющую следующий вид на трех временных слоях:
.
В итоге для множителя роста r получим квадратное уравнение
. (48)
Анализ корней уравнения (48) показывает, что разностная схема (47) устойчива при условии, что
. (49)
Условие устойчивости (49) является аналогом условия Куранта (11).
Изучим схему крест на примере численного решения двумерной задачи (45), (46), когда k1 = k2 = c 2 = const > 0, т.е. решим уравнение
. (50)
Положим вид правой части, начальные и краевые условия в прямоугольной области G(x1,x2) = [0,a]´[0,b] следующего вида:
(51)
Нетрудно проверить, что задача (50), (51) в области G имеет аналитическое решение
. (52)
(53)
На листинге_№6 приведен код программы решения задачи (50), (51) с помощью разностной схемы (53).
%Программа решения двумерного волнового уравнения
%(50), (51) с помощью явной схемы крест (53)
%Определяем габариты области интегрирования
%(t,x1,x2)=[0,T]x[0,a]x[0,b] и скорость переноса c
%Определяем число сгущений сеток по времени и
%Определяем цикл расчетов с различными сетками
%Определяем число узлом по времени Nt и по
%пространству N и M
%Определяем шаги по времени и пространству
%Определяем сетки по времени и пространству
t=0:tau:T; x1=0:h1:a; x2=0:h2:b;
%Используем начальное условие в (51) для
%определения численного решения на первом слое
%Используем начальное условие в (51) для
%определения численного решения на втором слое
%со вторым порядком точности по времени
%Удовлетворяем граничным условиям (51) при
%Удовлетворяем граничным условиям (51) при
%Организуем основной цикл расчета по времени с
%помощью явной схемы крест (53)
%Определяем разницу между численным решением и
%аналитическим решениями в норме C, т.е. определяем
%ошибку численного решения
%Делим ошибку численного решения на квадрат шага
%сетки по пространству
%Рисуем численное решение на момент времени T
%Рисуем график зависимости предстепенной константы
%Определяем функцию правой части уравнения (50)
%Определяем аналитическое решение (52)
На рис.10 приведен итог работы кода программы листинга_№6. На левом рисунке изображено численное решение yn,m в момент времени t = T, полученное с помощью численного расчета по разностной схеме (53). На правом рисунке изображена динамика зависимости предстепенной константы cost от шага сетки h1 в оценке ошибки численного решения вида:
,
где 
— аналитическое решение (52), при этом считается, что t
h1. Анализ правого графика показывает, что предстепенная константа немного растет. К сожалению, имеющиеся вычислительные ресурсы не позволяют за разумное время провести расчет для более подробных сеток и добиться более очевидного выхода величины const на некоторое постоянное значение при стремлении шага сетки к нулю. Это явилось бы подтверждением скорости сходимости схемы (53) — 
.
Рис.10. Численное решение двумерного волнового уравнения (50), (51) с
помощью явной схемы крест (53)
Факторизованные схемы. Ограничение на шаг по времени, возникающее в связи с условием Куранта (49), в явной схеме крест предыдущего численного примера было достаточно обременительным. Это ограничение может быть преодолено путем использования безусловно устойчивых многомерных экономичных схем, которые называются также факторизованными схемами или схемами с расщеплением. Построим такие схемы.
Для численного решения многомерного волнового уравнения (45) рассмотрим следующую разностную схему, называемую в дальнейшем исходной:
. (54)
Операторы La в (54) — трехточечные и вычисляются по формуле, аналогичной (18¢). Поскольку схема (54) симметрична по пространству и времени, постольку она имеет порядок аппроксимации при любых значениях веса s. Методом разделения переменных можно показать, что схема (54) является безусловно устойчивой при s ³ ¼.
Перепишем схему (54) в виде
, (55)
. (56)
Оператор B в (56) в той же форме (11.47) уже встречался у нас при изучении локально-одномерного метода решения многомерного параболического уравнения. Обращение оператора B сводится к обращению ленточной матрицы, внешний вид которой приведен на рис.6 лекции №12. Для двумерного случая данная матрица является пятидиагональной, к которой метод прогонки не применим. Таким образом, оператор B является труднообратимым, а разностная схема (54), (55) неэкономичной.
Оператор (56) можно приближенно заменить факторизованным оператором
(57)
т.е. приближенно расщепить на произведение одномерных сомножителей. Заменяя в исходной схеме (55) оператор B на оператор C, получим факторизованную схему:
(58)
которая отличается от исходной (54). Исследуем схему (58).
Аппроксимация. Учитывая (57), раскроем произведение в левой части (58), тогда, после несложных преобразований, получим
(59)
Сравнивая схемы (54) и (59) убеждаемся, что они различаются на члены порядка 
. Поскольку исходная схема (54) имеет второй порядок аппроксимации по времени и пространственным координатам, постольку и схема (58) также будет иметь аппроксимацию 
.
Устойчивость, как обычно, исследуем методом разделения переменных, подставляя в (58) многомерную гармонику, которая ранее уже была использована в схеме крест. Для множителя роста r можно получить квадратное уравнении
, (60)
(60¢)
Уравнение (60) аналогично уравнению (30), поэтому оба его корня по модулю не превышают единицы, если выполняются неравенства (32):
.
Согласно (60¢) первое из неравенств выполняется всегда. Второе неравенство заменим на более жесткое |m| £ n, которое, как нетрудно проверить, выполняется при условии, что
. (61)
Неравенство (61) является достаточным условием устойчивости разностной схемы (58). В частности, когда s ³ ¼ неравенство (61) выполняется при любых шагах по времени и пространству, т.е. схема (58) является безусловно устойчивой.
Безусловная сходимость факторизованной схемы (58) со скоростью следует из доказанных выше свойств аппроксимации и безусловной устойчивости при s ³ ¼.
Вычисление разностного решения сводится к последовательности одномерных прогонок по всем пространственным направления xa, a = 1,…,p. В итоге, можно констатировать, что для многомерного волнового уравнения существуют экономичные разностные схемы, сходящиеся со скоростью .
Изучим факторизованную разностную схему (58) на примере численного решения задачи (50), (51). Перепишем схему (58) для данного случая в виде системы двух уравнений, раскрывая тем самым преимущества факторизации:
(62)
В систему (62) добавлено вспомогательное решение w. Для численного решения систему (62) необходимо преобразовать к более подробному виду:
(63)
. (63¢)
Учитывая (51), (63¢), граничные условия без потери точности для вспомогательного решения w можно представить в следующем виде:
(64)
Расчет по схеме (63), (63¢) состоит из двух этапов. На первом этапе прогонками по направлению x1 решается уравнение (63) и находится w. На втором этапе прогонками по направлению x2 решается уравнение (63¢) и находится 
. На листинге_№7 приведен код программы, которая решает задачу (50), (51) с помощью факторизованной схемы (63), (63¢), (64) и оценивает скорость сходимости схемы расщепления.
%Программа решения двумерного волнового уравнения
%(50), (51) с помощью факторизованной схемы
%Определяем габариты области интегрирования
%(t,x1,x2)=[0,T]x[0,a]x[0,b] и скорость переноса c
%Определяем число сгущений сеток по времени и
%Определяем цикл расчетов с различными сетками
%Определяем число узлом по времени Nt и по
%пространству N и M
%Определяем шаги по времени и пространству
%Определяем сетки по времени и пространству
t=0:tau:T; x1=0:h1:a; x2=0:h2:b;
%Используем начальное условие из (51) для
%определения численного решения на первом слое
%Используем начальное условие из (51) для
%определения численного решения на втором слое
%со вторым порядком точности по времени
%Удовлетворяем граничным условиям из (51) при
%Удовлетворяем граничным условиям из (51) при
%Организуем основной цикл расчета по времени с
%помощью факторизованной схемы (63), (63′)
%Решаем уравнение (63) путем осуществления
%прогонок по направлению x1
%Используем левое граничное условие из (64)
http://helpiks.org/7-57615.html