Двухслойная акустическая схема

Лекция №13

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

Схема “крест”

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

Типичной одномерной задачей является задача описания малых колебаний натянутой струны с распределенной по длине нагрузкой f(t,x):

; (1)

; (2)

. (3)

Уравнение колебания струны (1) дополняется начальными условиями (2), которые, в отличие от параболического уравнения требуют двух условий: начального смещения относительно положения равновесия u₀(x) и начальную скорость движения u₁(x). Кроме того, задаются краевые условия (3), которые называются краевыми условиями первого рода, они описывают смещение концов струны относительно положений равновесия. Краевые условия могут быть и иного рода.

Составим несложную, но эффективную разностную схему для численного решения задачи (1) — (3), выбирая для простоты равномерные по t и x сетки. В качестве шаблона разностной схемы возьмем, представленный на рис.1 шаблон в форме “креста”. Аппроксимируя производные в (1) конечными разностями, получим трехслойную схему следующего вида

, (4)

с граничными условиями

. (4¢)

Рис.1. Разностная схема (4) в форме креста

По форме шаблона схему (4) называют схемой “крест”. Исследуем схему крест.

Процедура вычисления решения выглядит следующим образом. На нулевом слое решение известно из начального условия

. (5)

На первом слое решение можно также вычислить по начальным данным. Простейший способ получения решения на первом слое состоит в аппроксимации второго начального условия в (2) согласно представлению:

. (6)

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

. (6¢)

Подставляя в (6¢) u_tt из (1), найдем

. (6¢¢)

Схема крест (4) является явной, поскольку позволяет выразить решение на следующем слое через значения y_n и на двух предыдущих слоях. Поэтому, зная значения решения на нулевом (5) и первом слое (6¢¢), можно вычислить решения на всех последующих слоях. Итак, разностное решение существует и единственно.

Для изучения аппроксимации схемы крест, разложим точное решение по формуле Тейлора с центром в узле (t_m,x_n), считая все появляющиеся в разложении производные непрерывными:

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

(7)

а также невязку начального условия (6):

,

или более точного начального условия (6¢¢):

(7¢)

Начальные (2) и краевые условия (3) аппроксимируются точно. В итоге разностная схема (4) с начальными условиями (5), (6¢¢) имеет порядок аппроксимации , а та же разностная схема с начальным условием (6¢) имеет несколько худший порядок — .

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

. (8)

Подставляя представление (8) в (4), для множителя роста гармоник находим квадратное уравнение

. (9)

По теореме Виета произведение корней квадратного уравнения (9) равно единице, т.е. . Условие устойчивости может таким образом быть выполнено при . Это означает, что корни уравнения (9) должны образовывать комплексно сопряженную пару. Для этого необходимо, чтобы дискриминант уравнения (9) был неположительным для всех гармоник, т.е.

. (10)

Для выполнения неравенства (10) необходимо и достаточно соблюдение условия Куранта:

ct < h. (11)

В (11) стоит строгое неравенство, т.к., если верно равенство Куранта ct = h для некоторых гармоник, то появляется слабая неустойчивость из-за того, что корень квадратного уравнения (9) становится кратным.

В итоге разностная схема крест (4) с начальными условиями (5), (6¢) при выполнении условия Куранта (11) сходится с порядком аппроксимации . Из наших рассуждений следует, что сходимость имеет место в норме . Схема (4) обеспечивает хорошую точность расчета решения u(t,x), имеющих четвертую непрерывную производную. Схема позволяет также рассчитывать и менее гладкие решения.

Изучим аппроксимацию и сходимость схемы крест на численном примере решения уравнения (1) с правой частью, начальными и граничными условиями вида:

(12)

Нетрудно проверить, что решение задачи (1), (12) имеет следующее аналитическое решение:

. (13)

На листинге_№1 приведен код программы численного решения задачи (1), (12) с помощью разностной схемы крест (4) и начальными условиями (5), (6¢¢), т.е. с порядком аппроксимации .

Листинг_№1

%Программа решения волнового уравнения (1), (12) с

%помощью схемы крест (4)

function cross

global c

%Определяем габариты области интегрирования по времени T

%и пространству a, а также параметр c

T=1; a=2*pi; c=1;

%Определяем максимальное число удвоений числа узлов по

%времени и пространству smax

smax=6; N=2;

%Организуем основной цикл расчетов с различными сетками

for s=1:smax

%Определяем число узлов по пространству и удвоенное

%число узлов по времени

N=2*N; M=2*N;

%Определяем шаги по времени и пространству

tau=T/(M-1); h=a/(N-1);

r=(c^2*tau^2)/h^2;

%Определяем сетки по времени и пространству

t=0:tau:T; x=0:h:a;

%Используем начальные данные из (12) для

%определения численного решения на первом и

%втором слое по формуле (6'')

for n=1:N

y(1,n)=-x(n)*sin(x(n));

y(2,n)=y(1,n)+tau*c*x(n)*cos(x(n))+...

0.5*tau^2*(c^2*(-2*cos(x(n))+...

x(n)*sin(x(n)))+f(0,x(n)));

end

%Определяем левое и правое граничные условия,

%взятые из (12)

for m=1:M

y(m,1)=0; y(m,N)=a*sin(c*t(m)-a);

end

%Применяем схему крест (4) к нашей задаче

for m=2:(M-1)

for n=2:(N-1)

y(m+1,n)=2*y(m,n)-y(m-1,n)+r*(y(m,n+1)-...

2*y(m,n)+y(m,n-1))+tau^2*f(t(m),x(n));

end

end

%Оцениваем зависимость предстепенной константы

%const=||y-u||/h^2 ошибки численного решения на

%последнем шаге по времени в норме C от h

for n=1:N

z1(n)=abs(y(M,n)-u(t(M),x(n)));

end

const(s)=max(z1)/h^2;

step(s)=h;

end

mi=0;

for m=1:5:M

mi=mi+1; ti(mi)=t(m); ni=0;

for n=1:5:N

ni=ni+1; xi(ni)=x(n); z(mi,ni)=y(m,n);

end

end

%Рисуем трехмерную поверхность решения в координатах

%время-пространство

subplot(1,2,1); surf(xi,ti,z);

%Рисуем график зависимости предстепенной константы от

%шага сетки

subplot(1,2,2); loglog(step,const);

%Определяем функцию правой части

function y=f(t,x)

global c

y=2*c^2*cos(c*t-x);

%Определяем аналитическое решение

function y=u(t,x)

global c

y=x*sin(c*t-x);

Рис.2. Численное решение волнового уравнения (1), (12) с помощью
схемы крест (4), (5), (6¢¢)

На рис.2 приведен итог работы кода программы листинга_№1. На левом рисунке приведен внешний вид численного решения , m = 1,…,M. На правом рисунке приведена зависимость предстепенной константы const в представлении зависимости ошибки численного решения от шага сетки, при этом шаг по времени выбирался порядка шага по пространству (при выполнении условия Куранта), т.е. t ~ h. Видно, что по мере уменьшения шага h величина const действительно выходит на некоторое постоянное значение, что подтверждает второй порядок аппроксимации схемы крест (4) с начальными условиями (5), (6¢¢).

Неявная схема

В схеме крест условие Куранта (11), обеспечивающее устойчивость, может быть достаточно обременительным. По этой причине построим неявную схему для задачи (1) — (3), которая будет безусловно устойчивой.

На рис.3 приведен шаблон искомой схемы. Составим по шаблону рис.3 следующую разностную схему с весами:

(14)

Рис.3. Шаблон разностной схемы (14)

Веса схемы (14) веса неотрицательны при 0 £ s £ ½. При s = 0 неявная схема (14) переходит в схему крест. Граничные и начальные условия можно определять по типу (4¢) и (5), (6¢¢) соответственно. На третьем и последующих слоях разностная схема представляет собой линейную систему уравнений с трехдиагональной матрицей, в которой имеет место диагональное преобладание, т.е. решение системы уравнений существует, единственно и вычисляется методом прогонки.

Аналогично схеме крест (4) можно показать, что неявная схема (14) на решениях с четвертыми непрерывными производными аппроксимирует задачу (1) — (3) с порядком при любых значениях параметра s.

Устойчивость схемы (14) изучим методом разделения переменных. Подставляя представление (8) в (14), получим уравнение для оценки множителя роста гармоники следующего вида:

(15)

Проводя рассуждения аналогичные тем, которые проводились при изучении уравнения (9), можно сделать вывод о том, что условие устойчивости выполняется, когда корни уравнения (15) комплексно сопряженные, что возможно только при . Из последнего неравенства следует условие устойчивости разностной схемы (14):

. (16)

Из неравенства (16) следует, что при s ³ ¼ схема (14) безусловно устойчива, когда s < ¼ схема (14) условно устойчива при .

В итоге, разностная схема (14) при выборе параметра веса ¼ £ s £ ½, безусловно сходится с точностью .

Решим численно задачу (1), (12), (13) с помощью неявной разностной схемы (14). На листинге_№2 приведен код соответствующей программы.

Листинг_№2

%Программа решения волнового уравнения (1), (12), (13)

%с помощью неявной схемы (14)

function implicit

global c

%Определяем габариты области интегрирования по времени T

%и пространству a, а также параметр c

T=10; a=2*pi; c=1;

%Определяем максимальное число удвоений числа узлов по

%времени и пространству smax

smax=6; N=4;

%Определяем весовую константу

sigma=0.25;

%Организуем основной цикл расчетов с различными сетками

for s=1:smax

%Определяем число узлов по пространству и равное

%число узлов по времени

N=2*N; M=N;

%Определяем шаги по времени и пространству

tau=T/(M-1); h=a/(N-1);

r=(sigma*c^2*tau^2)/h^2;

r2=((1-2*sigma)*c^2*tau^2)/h^2;

%Определяем сетки по времени и пространству

t=0:tau:T; x=0:h:a;

%Используем начальные данные из (12) для

%определения численного решения на первом и

%втором слое по формуле (6'')

for n=1:N

y(1,n)=-x(n)*sin(x(n));

y(2,n)=y(1,n)+tau*c*x(n)*cos(x(n))+...

0.5*tau^2*(c^2*(-2*cos(x(n))+...

x(n)*sin(x(n)))+f(0,x(n)));

end

%Определяем левое и правое граничные условия,

%взятые из (12)

for m=1:M

y(m,1)=0; y(m,N)=a*sin(c*t(m)-a);

end

%Применяем неявную схему (14) к нашей задаче

for m=2:(M-1)

alpha(2)=0; beta(2)=y(m+1,1);

for n=2:(N-1)

alpha(n+1)=r/(1+r*(2-alpha(n)));

beta(n+1)=(r2*y(m,n-1)+(2-2*r2)*y(m,n)+...

r2*y(m,n+1)+r*y(m-1,n-1)-(1+2*r)*...

y(m-1,n)+r*y(m-1,n+1)+tau^2*...

f(t(m),x(n))+r*beta(n))/...

(1+r*(2-alpha(n)));

end

for n=N:-1:2

y(m+1,n-1)=alpha(n)*y(m+1,n)+beta(n);

end

end

%Оцениваем зависимость предстепенной константы

%const=||y-u||/h^2 ошибки численного решения на

%последнем шаге по времени в норме C от h

for n=1:N

z1(n)=abs(y(M,n)-u(t(M),x(n)));

end

const(s)=max(z1)/h^2;

step(s)=h;

end

mi=0;

for m=1:10:M

mi=mi+1; ti(mi)=t(m); ni=0;

for n=1:10:N

ni=ni+1; xi(ni)=x(n); z(mi,ni)=y(m,n);

end

end

%Рисуем трехмерную поверхность решения в координатах

%время-пространство

subplot(1,2,1); surf(xi,ti,z);

%Рисуем график зависимости предстепенной константы от

%шага сетки

subplot(1,2,2); loglog(step,const);

%Определяем функцию правой части

function y=f(t,x)

global c

y=2*c^2*cos(c*t-x);

%Определяем аналитическое решение

function y=u(t,x)

global c

y=x*sin(c*t-x);

На рис.4 приведен итог работы кода программы листинга_№2. На левом рисунке приведен внешний вид численного решения , m = 1,…,M. Обратим внимание на то, что использование неявной безусловно устойчивой схемы (14) позволило при том же приблизительно количестве узлов увеличить отрезок интегрирования по времени в десять раз по сравнению со схемой крест, которая является условно устойчивой. На правом рисунке приведена зависимость предстепенной константы const в представлении зависимости ошибки численного решения от шага сетки, при этом шаг по времени выбирался порядка шага по пространству, т.е. t ~ h. Видно, что по мере уменьшения шага h величина const действительно выходит на некоторое постоянное значение, что подтверждает второй порядок аппроксимации неявной схемы (14) с начальными условиями (5), (6¢¢), (12).

Рис.4. Решение волнового уравнения (1), (12), (13) с помощью неявной
разностной схемы (14)

Обобщим неявную разностную схему (14) на случай, когда скорость звука с является переменной. В этом случае уравнение (1) переписывается в виде:

, (17)

где k(t,x) º c²(t,x) > 0. Будем считать, что функции k(t,x), f(t,x) в (17) кусочно-непрерывны, причем разрывы неподвижны, т.е. лежат на линиях x = const. Предполагается, что на разрывах выполняется условие сопряжения [u] = [ku_x] = 0, т.е. решение u и поток ku_x являются непрерывными функциями.

Выберем по времени равномерную сетку, а по пространству специальную неравномерную сетку, у которой все точки разрыва коэффициентов являются узлами. Построим аналог наилучшей консервативной схемы (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)

где k₁, k₂, k₃, f₁, f₂, f₃ — некоторые постоянные величины. С учетом положений разрывов в (19), (20) определим равномерную сетку по пространству x_n = nh, n = 1,…,N, где N = 4l + 1, l = 1,2,… Определим также равномерную сетку по времени, т.е. t_m = t m, 1,…,M. С учетом сделанных предположений рассмотрим следующую схему аппроксимации для интегралов в (18¢):

. (21)

В качестве начальных и граничных условий положим

. (22)

На листинге_№3 приведен код программы численного решения задачи (17), (19), (20), (22) с помощью наилучшей разностной схемы (18), (18¢), (21).

Листинг_№3

%Программа решения волнового уравнения (17), (19),

%(20) с помощью разностной схемы (18), (18'), (21)

function best_scheme

global a k1 k2 k3 f1 f2 f3

%Определяем габариты области интегрирования по

%времени T и пространству a

T=5; a=1;

%Определяем константы k1,k2,k3,f1,f2,f3

k1=1; k2=10; k3=1; f1=1; f2=10; f3=1;

%Определяем весовую константу

sigma=0.25;

%Определяем число узлов по пространству и времени

M=201; N=161;

%Определяем шаги по времени и пространству

tau=T/(M-1); h=a/(N-1);

r=(sigma*tau^2)/h^2;

r2=((1-2*sigma)*tau^2)/h^2;

%Определяем сетки по времени и пространству

t=0:tau:T; x=0:h:a;

%Используем начальные данные из (22) для

%определения численного решения на первом и

%втором слое по формуле (6'')

for n=1:N

y(1,n)=0;

y(2,n)=y(1,n)+0.5*tau^2*f(x(n));

end

%Определяем левое и правое граничные условия,

%взятые из (22)

for m=1:M

y(m,1)=0; y(m,N)=0;

end

%Применяем неявную наилучшую схему (18), (18')

%к нашей задаче

for m=2:(M-1)

%Определяем коэффициенты прогонки

alpha(2)=0; beta(2)=y(m+1,1);

for n=1:(N-1)

kapa(n)=0.5*(k(x(n))+k(x(n+1)));

end

for n=2:(N-1)

alpha(n+1)=(r*kapa(n))/(1+r*(kapa(n)+...

kapa(n-1)*(1-alpha(n))));

beta(n+1)=(r2*kapa(n-1)*y(m,n-1)+...

(2-r2*(kapa(n-1)+kapa(n)))*y(m,n)+...

r2*kapa(n)*y(m,n+1)+r*kapa(n-1)*...

y(m-1,n-1)-(1+r*(kapa(n-1)+kapa(n)))*...

y(m-1,n)+r*kapa(n)*y(m-1,n+1)+0.25*...

tau^2*(f(x(n-1))+2*f(x(n))+f(x(n+1)))+...

r*kapa(n-1)*beta(n))/(1+r*(kapa(n)+...

kapa(n-1)*(1-alpha(n))));

end

for n=N:-1:2

y(m+1,n-1)=alpha(n)*y(m+1,n)+beta(n);

end

end

%Рисуем трехмерную поверхность решения в координатах

%время-пространство

surf(x,t,y);

%Определяем функцию квадрата скорости звука k(t,x)

function y=k(x)

global a k1 k2 k3

if (x>=0)&(x<=0.25*a)

y=k1;

end

if (x>0.25*a)&(x<0.75*a)

y=k2;

end

if (x>=0.75*a)&(x<=a)

y=k3;

end

%Определяем функцию правой части

function y=f(x)

global a f1 f2 f3

if (x>=0)&(x<=0.25*a)

y=f1;

end

if (x>0.25*a)&(x<0.75*a)

y=f2;

end

if (x>=0.75*a)&(x<=a)

y=f3;

end

Рис.5. Численное решение волнового уравнения (17), (19), (20), (22) с
помощью наилучшей разностной схемы (18), (18¢), (21)

Итог работы кода программы листинга_№3 приведен на рис.5. На графике представлено решение в виде волны, что, конечно же, характерно для волнового уравнения (17). Плоскости разрывов x = a/4 и x = 3a/4 просматриваются на поверхности решения.

Двухслойная акустическая схема

Волновое уравнение второго порядка может быть представлено в виде пары уравнений первого порядка. Рассмотрим систему уравнений

(23)

Покажем, что система (23), состоящая из пары уравнений первого порядка, сводится к волновому уравнению (1) второго порядка. Продифференцируем по x второе уравнение в (23) и подставим v_x из первого уравнения, тогда получим

.

Полное совпадение с волновым уравнением (1) наступит при условии, что v_x = u_t и F_x = f. Если проинтегрировать последние два уравнения по x, то

. (24)

Величины (24) называются потенциалами скорости и правой части.

Начальные данные (2) с учетом (24) можно переписать в виде:

, (25)

где V — некоторая константа.

Граничные условия (3) остаются без изменения, т.е.

. (26)

В ряде случаев задача акустики (23) — (26) оказывается более удобной для численного решения, чем волновое уравнение (1). Построим и исследуем неявную разностную схему для решения уравнений акустики.

Рис.6. Шаблон разностной схемы (27), (27¢)

Возьмем в общем случае неравномерную сетку по пространству 0 = x₀ < x₁ < … < x_N = a, в узлах которой определим величины , n = 0,1,…,N, а в серединах интервалов — величины , n = 0,1,…,N-1. Выберем шаблон разностной схемы, представленный на рис.6 и составим по нему разностную схему с весами:

, (27)

, (27¢)

где , и считается, что s₁, s₂ Î [0,1]. Шаблон схемы (27) выделен на рис.6 красным цветом (сплошные линии), шаблон схемы (27¢) — зеленым цветом (точечные линии).

Если положить , то схема (27) будет симметричной по переменной x и иметь порядок аппроксимации . Если положить s₁ = s₂ = ½ и , то порядок аппроксимации по времени схемы (27), (27¢) станет вторым, т.е. .

Как обычно устойчивость схемы (27), (27¢) изучим с помощью метода разделения переменных, представляя возмущения функций y и z в виде гармоник следующего вида:

. (28)

Гармоники для y и z в (28) взяты с одной и той же частотой и множителем роста, но с различными амплитудами. Подставляя (28) в (27), (27¢) и полагая , для амплитуд a и b получим однородную систему уравнений

(29)

Система линейных уравнений (29) имеет нетривиальное решение, когда ее определитель равен нулю. Это условие дает квадратное уравнение для определения множителя роста r:

, (30)

где

(31)

Каждый из двух корней уравнения (30) меньше по модулю единицы, если и только если

. (32)

Первое из неравенств в (32) следует из теоремы Виета , второе можно доказать с помощью несложных, но несколько громоздких выкладок. Неравенства (32) выполняются для всех гармоник, когда

, (33)

неравенства (33) являются условиями устойчивости для разностной схемы (27), (27¢).

Из неравенства (33) следует, что если s₁ ³ ½ и s₂ ³ ½, то схема (27), (27¢) является безусловно устойчивой. Если s₁ + s₂ ³ 1, но один из весов меньше ½, то схема условно устойчива . При s₁ + s₂ < 1 схема безусловно неустойчива.

Особо выделим симметричную схему, которая реализуется при s₁ = s₂ = ½. В этом случае разностная схема (27), (27¢) является безусловно устойчивой и сходится со скоростью . Данная схема является одной из лучших схем для задач акустики.

Разностная схема (27), (27¢) при произвольных значениях весов s₁ и s₂ решается путем определения из уравнения (27), т.е.

(34)

и подстановки и в (27¢), что приводит к уравнению

(34¢)

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

Численно изучим разностную схему (34), (34¢) на примере решения задачи акустики (24) при

, (35)

а также при начальных и граничных значениях вида:

(36)

Как нетрудно проверить уравнения акустики (23) при условии (35), (36) имеют следующее аналитическое решение:

. (37)

Аналитические решения (37) получены с помощь преобразования (24) при v(t,0) = c×cos(ct), исходя из решений задачи (1), (12), (13).

На листинге_№4 приведен код программы численного решения задачи акустики (23), (35), (36) с помощью разностной схемы (34), (34¢).

Листинг_№4

%Программа численного решения задачи акустики

%(23) - (26) с помощью разностной схемы (34), (34')

function acoustics

global c

%Определяем габариты области интегрирования

%G=[0,T]x[0,a] и скорость звука c

T=8; a=2*pi; c=1;

%Выбираем для анализа симметричную схему, у которой

%весовые коэффициенты

sigma1=0.5; sigma2=0.5;

%Количество удвоений smax числа узлов сетки

smax=6; N=4;

%Организуем цикл расчетов по разностной схеме

%(34), (34') с различным числом узлов по

%пространству и времени

for s=1:smax

%Удваиваем число узлов по пространству

N=2*N;

%Считаем, что число узлов по времени равно

%числу узлов по пространству

M=N;

%Определяем шаги по времени и пространству

tau=T/(M-1); h=a/(N-1);

%Определяем сетки по времени и пространству

t=0:tau:T; x=0:h:a;

r=(sigma1*sigma2*c^2*tau^2)/h^2;

r2=((1-sigma1)*sigma2*c^2*tau^2)/h^2;

%Определение начального распределения для u-y

%согласно (36)

for n=1:N

y(1,n)=u(t(1),x(n));

end

%Определение начального распределения для v-z

%согласно (36)

for n=1:(N-1)

z(1,n)=v(t(1),x(n)+0.5*h);

end

%Определение граничных условий для u-y

%согласно (36)

for m=1:M

y(m,1)=u(t(m),x(1));

y(m,N)=u(t(m),x(N));

end

%Организуем цикл расчетов по времени

for m=1:(M-1)

%Находим коэффициенты прогонки для вычисления u-y

alpha(2)=0; beta(2)=y(m+1,1);

for n=2:(N-1)

alpha(n+1)=r/(1+r*(2-alpha(n)));

beta(n+1)=(r2*y(m,n-1)+(1-2*r2)*y(m,n)+...

r2*y(m,n+1)+(tau/h)*(z(m,n)-z(m,n-1))+...

((sigma2*tau^2)/h)*(F(t(m)+0.5*tau,...

x(n)+0.5*h)-F(t(m)+0.5*tau,x(n-1)+...

0.5*h))+r*beta(n))/(1+r*(2-alpha(n)));

end

for n=N:-1:2

y(m+1,n-1)=alpha(n)*y(m+1,n)+beta(n);

end

%Вычисляем функцию v-z

for n=1:(N-1)

z(m+1,n)=z(m,n)+((c^2*tau)/h)*(sigma1*...

(y(m+1,n+1)-y(m+1,n))+(1-sigma1)*...

(y(m,n+1)-y(m,n)))+tau*F(t(m)+...

0.5*tau,x(n)+0.5*h);

end

end

%Находим ошибку численного решения в норме C

for m=1:M

for n=1:N

ery(m,n)=abs(y(m,n)-u(t(m),x(n)));

end

end

for m=1:M

for n=1:(N-1)

erz(m,n)=abs(z(m,n)-v(t(m),x(n)+0.5*h));

end

end

%Делим ошибку численного решения на h^2 и находим

%предстепенную константу в зависимости ошибки

%численного решения от шага сетки

const(s)=max(max(max(ery)),max(max(erz)))/h^2;

step(s)=h;

end

mi=0;

for m=1:10:M

mi=mi+1; ni=0;

for n=1:10:N

ni=ni+1; yi(mi,ni)=y(m,n);

end

end

%Рисуем численное решение функции u-y

subplot(1,3,1); surf(x([1:10:N]),t([1:10:M]),yi);

mi=0;

for m=1:10:M

mi=mi+1; ni=0;

for n=1:10:(N-1)

ni=ni+1; zi(mi,ni)=z(m,n);

end

end

%Рисуем численное решение функции v-z

subplot(1,3,2); surf(x([1:10:(N-1)]),t([1:10:M]),zi);

%Рисуем зависимость предстепенной константы

%const = max{||y-u||,||z-v||}/h^2 от шага сетки h

subplot(1,3,3); semilogx(step,const);

%Определяем функцию правой части

function y=F(t,x)

global c

y=-2*c^2*sin(c*t-x);

%Определяем аналитическое решение для u

function y=u(t,x)

global c

y=x*sin(c*t-x);

%Определяем аналитическое решение для v

function y=v(t,x)

global c

y=-c*x*sin(c*t-x)+c*cos(c*t-x);

На рис.7 приведен итог работы кода программы листинга_№4. На левом рисунке изображено численное решение u = u(t_m,x_n) » y(m,n), на среднем рисунке — v = v(t_m,x_n) » z(m,n). На правом рисунке изображена динамика предстепенной константы const в оценке ошибки численного решения симметричной схемы (s₁ = s₂ = ½) от шага по пространству h, при этом полагалось, что t ~ h. Видно, что при уменьшении шага сетки h предстепенная константа действительно выходит на некоторое постоянное значение.

Рис.7. Численное решение задачи акустики (24), (35), (36) с помощью
разностной схемы (34), (34¢)

Инварианты. Перепишем уравнения акустики с помощью инвариантов:

. (38)

Умножим первое уравнение в системе уравнений акустики (23) на c и в начале сложим оба уравнения друг с другом, а затем вычтем из второго уравнения первое, тогда получим

(39)

С учетом (25) нетрудно найти начальные данные для инвариантов:

. (40)

С учетом (26) определим краевые условия для инвариантов:

. (41)

Инвариант r удовлетворяет уравнению переноса вправо (c > 0), инвариант s — уравнению переноса влево. Для однородной задачи, когда F = 0, m₁ = m₂ =0, величины r, s переносятся по своим характеристикам без изменения, с чем и связано их наименование.

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

. (42)

Схема (42) является схемой бегущего счета. Нетрудно показать, что при выполнении условия Куранта ct £ h схема устойчива и сходится с порядком точности на дважды непрерывно дифференцируемых функциях.

Изучим разностную схему для инвариантов на примере численного решения задачи (39) с помощью разностной схемы (42), выбирая правую часть, начальные и граничные условия согласно (35), (36). В этом случае при учете (40), (41), находим следующие выражения для правой части, начальных и граничных условий:

(43)

Нетрудно проверить, что аналитическим решением задачи (39), (43) являются выражения вида:

. (44)

На листинге_№5 приведен код программы расчета задачи акустики в инвариантах (39) со специальной правой частью, начальными и граничными условиями (43).

Листинг_№5

%Программа решения уравнений акустики (39) в

%инвариантах с правой частью, начальными и

%граничными условиями вида (43)

function invariants

global c

%Определяем габариты области интегрирования по

%времени T и пространству a, а также определяем

%величину скорости c

T=2*pi; a=2*pi; c=1;

%Определяем число удвоений числа узлов сетки по

%времени и пространству dmax

dmax=5; N=4;

%Организуем цикл расчетов на различных сетках

for d=1:dmax

%Определяем число узлов сетки по пространству

%и времени

N=2*N; M=N;

%определяем шаги по времени и пространству

tau=T/(M-1); h=a/(N-1); kur=(c*tau)/h;

%Определяем сетки по времени и пространству

t=0:tau:T; x=0:h:a;

%Определяем начальные условия (43) для

%инвариантов r и s

for n=1:N

r(1,n)=2*c*x(n)*sin(x(n))+c*cos(x(n));

s(1,n)=c*cos(x(n));

end

%Организуем цикл расчетов по времени

for m=1:(M-1)

%Осуществляем бегущий расчет инварианта r,

%осуществляемый слева направо

for n=2:N

r(m+1,n)=(1-kur)*r(m,n)+kur*r(m,n-1)+...

tau*F(t(m),x(n));

end

%Учитываем правое граничное условие (43)

s(m+1,N)=r(m+1,N)+2*c*a*sin(c*t(m+1)-a);

%Осуществляем бегущий расчет инварианта s,

%осуществляемый справа налево

for n=(N-1):-1:1

s(m+1,n)=(1-kur)*s(m,n)+kur*s(m,n+1)+...

tau*F(t(m),x(n));

end

%Учитываем левое граничное условие (43)

r(m+1,1)=s(m+1,1);

end

%Ошибку численного расчета инвариантов r и s

%оцениваем как разность численного и

%аналитического решений в норме C. Делим

%полученную ошибку на шаг по пространству h и

%находим предстепенную константу скорости

%сходимости схемы с инвариантами

for m=1:M

for n=1:N

er_r(m,n)=abs(r(m,n)-ra(t(m),x(n)));

er_s(m,n)=abs(s(m,n)-sa(t(m),x(n)));

end

end

const(d)=max(max(max(er_r)),max(max(er_s)))/h;

step(d)=h;

end

mi=0;

for m=1:5:M

mi=mi+1; ni=0;

for n=1:5:N

ni=ni+1;

ri(mi,ni)=r(m,n);

si(mi,ni)=s(m,n);

end

end

%Рисуем численное решение инварианта r

subplot(1,3,1); surf(x([1:5:N]),t([1:5:M]),ri);

%Рисуем численное решение инварианта s

subplot(1,3,2); surf(x([1:5:N]),t([1:5:M]),si);

%Рисуем график зависимости предстепенной

%константы от шага сетки h

subplot(1,3,3); semilogx(step,const);

%Определяем функцию правой части

function y=F(t,x)

global c

y=-2*c^2*sin(c*t-x);

%Определяем аналитический вид инварианта r

function y=ra(t,x)

global c

y=-2*c*x*sin(c*t-x)+c*cos(c*t-x);

%Определяем аналитический вид инварианта s

function y=sa(t,x)

global c

y=c*cos(c*t-x);

Рис.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)

при этом в случае переменных коэффициентов k_a операторы L_a определяются также как и для наилучшей схемы (18), (18¢).

Рис.9. Шаблон схемы крест в двумерном случае

Трехслойная схема (47) является явной. Вычисления с помощью схемы (47) одинаково просты для любого числа измерений. Легко проверить, что схема имеет порядок аппроксимации .

Исследуем устойчивость схемы (47) методом разделения переменных, считая коэффициенты k_a постоянными. Для этого подставим в (47) многомерную гармонику, имеющую следующий вид на трех временных слоях:

.

В итоге для множителя роста r получим квадратное уравнение

. (48)

Анализ корней уравнения (48) показывает, что разностная схема (47) устойчива при условии, что

. (49)

Условие устойчивости (49) является аналогом условия Куранта (11).

Изучим схему крест на примере численного решения двумерной задачи (45), (46), когда k₁ = k₂ = c² = const > 0, т.е. решим уравнение

. (50)

Положим вид правой части, начальные и краевые условия в прямоугольной области G(x₁,x₂) = [0,a]´[0,b] следующего вида:

(51)

Нетрудно проверить, что задача (50), (51) в области G имеет аналитическое решение

. (52)

Определим сетку по пространству: (x_1,n,x_2,m) = (h₁n,h₂m), n = 0,1,…,N, m = 0,1,…,M, где шаги по направлениям x₁ и x₂ определяются согласно формулам: h₁ = a/N, h₂ = b/M. Запишем разностную схему для уравнения (50), тогда

(53)

На листинге_№6 приведен код программы решения задачи (50), (51) с помощью разностной схемы (53).

Листинг_№6

%Программа решения двумерного волнового уравнения

%(50), (51) с помощью явной схемы крест (53)

function multi_cross

global c

%Определяем габариты области интегрирования

%(t,x1,x2)=[0,T]x[0,a]x[0,b] и скорость переноса c

T=0.5; a=1; b=1; c=1;

%Определяем число сгущений сеток по времени и

%пространству smax

smax=5; N=4;

%Определяем цикл расчетов с различными сетками

for s=1:smax

%Определяем число узлом по времени Nt и по

%пространству N и M

N=N+10; M=N; Nt=N;

%Определяем шаги по времени и пространству

tau=T/(Nt-1); h1=a/(N-1); h2=b/(M-1);

r1=(c^2*tau^2)/h1^2; r2=(c^2*tau^2)/h2^2;

%Определяем сетки по времени и пространству

t=0:tau:T; x1=0:h1:a; x2=0:h2:b;

%Используем начальное условие в (51) для

%определения численного решения на первом слое

for n=1:N

for m=1:M

y(1,n,m)=x1(n)^2*x2(m)^2;

end

end

%Используем начальное условие в (51) для

%определения численного решения на втором слое

%со вторым порядком точности по времени

for n=1:N

for m=1:M

y(2,n,m)=x1(n)^2*x2(m)^2-...

2*c*tau*x1(n)*x2(m)*(x1(n)+x2(m))+...

0.5*tau^2*(2*c^2*(x1(n)^2+x2(m)^2)+...

f(t(1),x1(n),x2(m)));

end

end

%Удовлетворяем граничным условиям (51) при

%x2=0 и x2=b

for tm=1:Nt

for n=1:N

y(tm,n,1)=(c*t(tm)-x1(n))^2*c^2*t(tm)^2;

y(tm,n,M)=(c*t(tm)-x1(n))^2*(c*t(tm)-b)^2;

end

end

%Удовлетворяем граничным условиям (51) при

%x1=0 и x1=a

for tm=1:Nt

for m=1:M

y(tm,1,m)=c^2*t(tm)^2*(c*t(tm)-x2(m))^2;

y(tm,N,m)=(c*t(tm)-a)^2*(c*t(tm)-x2(m))^2;

end

end

%Организуем основной цикл расчета по времени с

%помощью явной схемы крест (53)

for tm=2:(Nt-1)

for n=2:(N-1)

for m=2:(M-1)

y(tm+1,n,m)=2*y(tm,n,m)-y(tm-1,n,m)+...

r1*(y(tm,n-1,m)-2*y(tm,n,m)+...

y(tm,n+1,m))+r2*(y(tm,n,m-1)-...

2*y(tm,n,m)+y(tm,n,m+1))+...

tau^2*f(t(tm),x1(n),x2(m));

end

end

end

%Определяем разницу между численным решением и

%аналитическим решениями в норме C, т.е. определяем

%ошибку численного решения

for tm=1:Nt

for n=1:N

for m=1:M

er(tm,n,m)=abs(y(tm,n,m)-...

u(t(tm),x1(n),x2(m)));

end

end

end

%Делим ошибку численного решения на квадрат шага

%сетки по пространству

const(s)=max(max(max(er)))/h1^2;

step(s)=h1;

end

for n=1:N

for m=1:M

z(n,m)=y(Nt,n,m);

end

end

%Рисуем численное решение на момент времени T

subplot(1,2,1); surf(x2,x1,z);

%Рисуем график зависимости предстепенной константы

%от шага сетки

subplot(1,2,2); plot(step,const);

%Определяем функцию правой части уравнения (50)

function y=f(t,x1,x2)

global c

y=8*c^2*(c*t-x1)*(c*t-x2);

%Определяем аналитическое решение (52)

function y=u(t,x1,x2)

global c

y=(c*t-x1)^2*(c*t-x2)^2;

На рис.10 приведен итог работы кода программы листинга_№6. На левом рисунке изображено численное решение y_n,m в момент времени t = T, полученное с помощью численного расчета по разностной схеме (53). На правом рисунке изображена динамика зависимости предстепенной константы cost от шага сетки h₁ в оценке ошибки численного решения вида:

,

где — аналитическое решение (52), при этом считается, что t ~ h₁, h₂ ~ h₁. Анализ правого графика показывает, что предстепенная константа немного растет. К сожалению, имеющиеся вычислительные ресурсы не позволяют за разумное время провести расчет для более подробных сеток и добиться более очевидного выхода величины const на некоторое постоянное значение при стремлении шага сетки к нулю. Это явилось бы подтверждением скорости сходимости схемы (53) — .

Рис.10. Численное решение двумерного волнового уравнения (50), (51) с
помощью явной схемы крест (53)

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

Для численного решения многомерного волнового уравнения (45) рассмотрим следующую разностную схему, называемую в дальнейшем исходной:

. (54)

Операторы L_a в (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 ³ ¼.

Вычисление разностного решения сводится к последовательности одномерных прогонок по всем пространственным направления x_a, a = 1,…,p. В итоге, можно констатировать, что для многомерного волнового уравнения существуют экономичные разностные схемы, сходящиеся со скоростью .

Изучим факторизованную разностную схему (58) на примере численного решения задачи (50), (51). Перепишем схему (58) для данного случая в виде системы двух уравнений, раскрывая тем самым преимущества факторизации:

(62)

В систему (62) добавлено вспомогательное решение w. Для численного решения систему (62) необходимо преобразовать к более подробному виду:

(63)

. (63¢)

Учитывая (51), (63¢), граничные условия без потери точности для вспомогательного решения w можно представить в следующем виде:

(64)

Расчет по схеме (63), (63¢) состоит из двух этапов. На первом этапе прогонками по направлению x₁ решается уравнение (63) и находится w. На втором этапе прогонками по направлению x₂ решается уравнение (63¢) и находится . На листинге_№7 приведен код программы, которая решает задачу (50), (51) с помощью факторизованной схемы (63), (63¢), (64) и оценивает скорость сходимости схемы расщепления.

Листинг_№7

%Программа решения двумерного волнового уравнения

%(50), (51) с помощью факторизованной схемы

%(63), (63')

function factorization

global c

%Определяем габариты области интегрирования

%(t,x1,x2)=[0,T]x[0,a]x[0,b] и скорость переноса c

T=0.75; a=1; b=2; c=1;

%Определяем вес

sigma=0.25;

%Определяем число сгущений сеток по времени и

%пространству smax

smax=5; N=4;

%Определяем цикл расчетов с различными сетками

for s=1:smax

%Определяем число узлом по времени Nt и по

%пространству N и M

N=N+10; M=N; Nt=N;

%Определяем шаги по времени и пространству

tau=T/(Nt-1); h1=a/(N-1); h2=b/(M-1);

r1=(sigma*c^2*tau^2)/h1^2;

r2=(sigma*c^2*tau^2)/h2^2;

r3=((1-2*sigma)*c^2*tau^2)/h1^2;

r4=((1-2*sigma)*c^2*tau^2)/h2^2;

%Определяем сетки по времени и пространству

t=0:tau:T; x1=0:h1:a; x2=0:h2:b;

%Используем начальное условие из (51) для

%определения численного решения на первом слое

for n=1:N

for m=1:M

y(1,n,m)=x1(n)^2*x2(m)^2;

end

end

%Используем начальное условие из (51) для

%определения численного решения на втором слое

%со вторым порядком точности по времени

for n=1:N

for m=1:M

y(2,n,m)=x1(n)^2*x2(m)^2-...

2*c*tau*x1(n)*x2(m)*(x1(n)+x2(m))+...

0.5*tau^2*(2*c^2*(x1(n)^2+x2(m)^2)+...

f(t(1),x1(n),x2(m)));

end

end

%Удовлетворяем граничным условиям из (51) при

%x2=0 и x2=b

for tm=1:Nt

for n=1:N

y(tm,n,1)=(c*t(tm)-x1(n))^2*c^2*t(tm)^2;

y(tm,n,M)=(c*t(tm)-x1(n))^2*(c*t(tm)-b)^2;

end

end

%Удовлетворяем граничным условиям из (51) при

%x1=0 и x1=a

for tm=1:Nt

for m=1:M

y(tm,1,m)=c^2*t(tm)^2*(c*t(tm)-x2(m))^2;

y(tm,N,m)=(c*t(tm)-a)^2*(c*t(tm)-x2(m))^2;

end

end

%Организуем основной цикл расчета по времени с

%помощью факторизованной схемы (63), (63')

for tm=2:(Nt-1)

%Решаем уравнение (63) путем осуществления

%прогонок по направлению x1

for m=2:(M-1)

alpha(2)=0;

%Используем левое граничное условие из (64)

beta(2)=y(tm+1,1,m)-...

2*sigma*tau^2*c^4*t(tm+1)^2;

for n=2:(N-1)

alpha(n+1)=r1/(1+r1*(2-alpha(n)));

beta(n+1)=(2*y(tm,n,m)+r3*(y(tm,n-1,m)-...

2*y(tm,n,m)+y(tm,n+1,m))+...

r4*(y(tm,n,m-1)-2*y(tm,n,m)+...

y(tm,n,m+1))-y(tm-1,n,m)+...

r1*(y(tm-1,n-1,m)-2*y(tm-1,n,m)+...

y(tm-1,n+1,m))+r2*(y(tm-1,n,m-1)-...

2*y(tm-1,n,m)+y(tm-1,n,m+1))+...

tau^2*f(t(tm),x1(n),x2(m))+...

r1*beta(n))/(1+r1*(2-alpha(n)));

end

%Используем правое граничное условие из (64)

w(N,m)=y(tm+1,N,m)-...

2*sigma*tau^2*c^2*(c*t(tm+1)-a)^2;

%Вычисляем промежуточное решение w

for n=N:-1:2

w(n-1,m)=alpha(n)*w(n,m)+beta(n);

end

end

%Решаем уравнение (63') путем осуществления

%прогонок по направлению x2

for n=2:(N-1)

alpha(2)=0; beta(2)=y(tm+1,n,1);

for m=2:(M-1)

alpha(m+1)=r2/(1+r2*(2-alpha(m)));

beta(m+1)=(w(n,m)+r2*beta(m))/...

(1+r2*(2-alpha(m)));

end

for m=M:-1:2

y(tm+1,n,m-1)=...

alpha(m)*y(tm+1,n,m)+beta(m);

end

end

end

%Определяем разницу между численным решением и

%аналитическим решениями в норме C, т.е. определяем

%ошибку численного решения

for tm=1:Nt

for n=1:N

for m=1:M

er(tm,n,m)=abs(y(tm,n,m)-...

u(t(tm),x1(n),x2(m)));

end

end

end

%Делим ошибку численного решения на квадрат шага

%сетки по пространству

const(s)=max(max(max(er)))/h1^2;

step(s)=h1;

end

for n=1:N

for m=1:M

z(n,m)=y(Nt,n,m);

end

end

%Рисуем численное решение на момент времени T

subplot(1,2,1); surf(x2,x1,z);

%Рисуем график зависимости предстепенной константы

%от шага сетки

subplot(1,2,2); plot(step,const);

%Определяем функцию правой части уравнения (50)

function y=f(t,x1,x2)

global c

y=8*c^2*(c*t-x1)*(c*t-x2);

%Определяем аналитическое решение (52)

function y=u(t,x1,x2)

global c

y=(c*t-x1)^2*(c*t-x2)^2;

На рис.11 приведен итог работы кода программы листинга_№7. На левом рисунке представлено трехмерное изображение численного решения y_n,m в момент времени t = T, полученное с помощью расчета по факторизованной разностной схеме (63), (63¢), (64). На правом рисунке изображена динамика зав