Геометрическая интерпретация устойчивости




Дадим геометрическую интерпретацию устойчивости по начальным данным. Для простоты рассмотрим однородное уравнение (3), т.е. при f (t, x) = 0. Общее решение однородного уравнения переноса имеет вид (4): , т.е. значения решения переносятся вдоль характеристики x - ct = const без изменения.

Рассмотрим разностную схему (7а) с шаблоном на рис.6,а (аналогичный шаблон приведен на рис.2,а). На рис.6,а стрелкой обозначена характеристика, приходящая в точку (tm +1, xn). Характеристика пересекает слой tm в точке . Точка пересечения находится согласно определению характеристики: = const.

 

Рис.6,а. Перенос значения функции по характеристике при выполнении условия Куранта (для схемы (7а)) Рис.6,б. Перенос значения функции по характеристике при нарушении условия Куранта (для схемы (7а)) Рис.6,в. Перенос значения функции по характеристике при любых значениях t и h (для схемы (7в))

 

Найдем решение в точке с помощью линейной интерполяции разностного решения между двумя узлами шаблона рис.6,а на слое tm, тогда

.

Далее полученное значение из точки переносим без изменения по характеристике в узел (tm +1, xn), т.е. положим .

Выполнение условия устойчивости для схемы (7а) ct £ h обеспечивает попадание точки в полуинтервал [ xn -1, xn), т.е. . Если условие устойчивости нарушается, точка становится меньше xn -1, т.е. при ct > h. В этой ситуации решение в точке находится с помощью экстраполяции. Другими словами, схема (7а) устойчива, если вычисляется по значениям с предыдущего слоя при помощи интерполяции (рис.6,а). Схема (7а) неустойчива, если при вычислении величины используется экстраполяция (рис.6,б).

Разностные схемы (7б), (7в) можно также интерпретировать как линейные интерполяции по двум уже вычисленным значениям с переносом полученного значения по характеристике. Так, безусловная устойчивость схемы (7в) обязана тому, что приходящая в узел (tm +1, xn) характеристика при любых t и h пересекает отрезок (пунктирная линия на рис.6,в), соединяющий узлы интерполяции значения решения, которое по характеристике переносится в узел (tm +1, xn).

Разностная схема (7г) также может быть истолкована как интерполяция, но не двухточечная, а трехточечная. Учитывая рис.2,г видно, что при любом выборе шагов t и h приходящая в узел (tm +1, xn) характеристика переносит значение, которое вычисляется с помощью интерполяции по ранее вычисленным значениям.

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

Явно-неявная схема. Будем считать, что шаг по времени и по пространству не постоянны, а скорость переноса переменная величина, т.е. c = c (t, x). При вычислении проверим условие Куранта (8) в текущей ячейке. Если условие Куранта выполняется, то используем схему (7а):

, (15)

если условие Куранта в форме (8) нарушается, используем схему (7б):

. (15¢)

Явно-неявная схема (15), (15¢) безусловно устойчива, причем невязка этой схемы меньше, чем у безусловно устойчивой схемы (7в). Схему (15), (15¢) используют в том случае, когда искомое решение является недостаточно гладким или быстропеременным.

Сравним разностные схемы (15), (15¢) и (7в) на примере решения задачи:

(16)

Решение задачи (16) можно легко найти

. (17)

На листинге_№4 приведен код программы численного решения задачи (16) и сравнения полученного решения с аналитическим решением (17). Сравнение производится в норме , т.е. error = .

 

Листинг_№4

%Программа сравнения явно-неявной схемы для решения

%уравнения переноса с безусловно устойчивой схемой (7в)

function obvious

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

%по времени и по пространству, [0,T]x[0,a]

T=1; a=1; kmax=300; l=1;

for k=50:50:kmax

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

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

Nt=k; Nx=k;

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

t(1)=0;

for m=1:(Nt-1)

t(m+1)=t(m)+(T/(Nt-1))*(1+0.1*randn);

end

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

x(1)=0;

for n=1:(Nx-1)

x(n+1)=x(n)+(a/(Nx-1))*(1+0.1*randn);

end

%Задаем левое граничное условие u(t,0)=-t^2

for m=1:Nt

y(m,1)=-t(m)^2;

end

%Задаем начальные данные u(0,x)=x^2

for n=1:Nx

y(1,n)=x(n)^2;

end

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

%схеме (15), (15')

for m=1:(Nt-1)

tau=t(m+1)-t(m);

for n=2:Nx

h=x(n)-x(n-1); c=t(m+1)/x(n); r=(tau*c)/h;

if r<=1

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

tau*f(t(m)+0.5*tau,x(n)-0.5*h);

else

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

(tau/r)*f(t(m)+0.5*tau,x(n)-0.5*h);

end

end

end

%Находим различие между численным и

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

for m=1:Nt

for n=1:Nx

yerror(m,n)=abs(y(m,n)-ya(t(m),x(n)));

end

end

%Выводим ошибку численного решения в норме C

error_obvious(l)=max(max(yerror));

%Организуем расчет по безусловно устойчивой схеме (7в)

for m=1:(Nt-1)

tau=t(m+1)-t(m);

for n=2:Nx

h=x(n)-x(n-1); c=t(m+1)/x(n); r=(tau*c)/h;

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

(tau/(1+r))*f(t(m)+0.5*tau,x(n)-0.5*h);

end

end

%Находим различие между численным и

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

for m=1:Nt

for n=1:Nx

yerror(m,n)=abs(y(m,n)-ya(t(m),x(n)));

end

end

%Выводим ошибку численного решения в норме C

error_non_obvious(l)=max(max(yerror));

l=l+1;

end

plot(50:50:kmax,error_obvious,'-*',...

50:50:kmax,error_non_obvious,'-p','MarkerSize',12);

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

function y=f(t,x)

y=x^2+2*t^2;

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

function y=ya(t,x)

y=x^2-t^2+x^2*t;

 

На рис.7 приведен итоговый график сравнения явно-неявной схемы (15), (15¢) и безусловно устойчивой неявной схемы (7в) на примере решения задачи (16), (17). На графике по оси абсцисс отложено число точек сетки по времени и по пространству, по оси ординат — ошибка как разность численного и аналитического решений в норме . Из графика отчетливо видно, что явно-неявная схема более точна, чем безусловно устойчивая неявная схема.

 

Рис.7. Сравнение ошибок при расчетах по явно-неявной (15), (15¢) и
безусловно устойчивой неявной (7в) разностных схем

 

Схема без шаблона. Рассмотрим характеристику, которая приходит в узел (tm +1, xn) и которая пересекает слой tm в точке . Найдем пару узлов на слое tm, между которыми попадает точка . Определим с помощью линейной интерполяции по значениям :

, (18)

при этом . Переносим значение по характеристике в узел (tm +1, xn), т.е. полагаем . В соответствии с геометрическим смыслом устойчивости, разностная схема (18) безусловно устойчива.

В разностной схеме (18) положение пары узлов и относительно узла (tm +1, xn) не определено, когда скорость c = c (t, x) не фиксирована или сетка по пространству не равномерна. Это означает, что разностная схема (18) является схемой без шаблона.

Изучим поведение точности схемы (18) на примере решения задачи:

(19)

Решение задачи (19) можно легко найти

. (20)

На листинге_№5 приведен код программы, которая с помощью схемы без шаблона численно решает задачу (19). На выходе работы программы строится график ошибки численного решения error = в зависимости от числа узлов сетки по времени и по пространству (Nt = Nx). Данный график приведен на рис.8. График на рис.8 демонстрирует “правильное” поведение ошибки, она уменьшается по мере роста числа узлов одновременно по времени и по пространству.

 

Листинг_№5

%Программа тестирования схемы без шаблона для

%численного решения уравнения переноса (19)

function witht_template

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

%переноса по времени и по пространству, [0,T]x[0,a]

T=1; a=1; kmax=400; l=1;

for k=50:50:kmax

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

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

Nt=k; Nx=k;

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

t(1)=0;

for m=1:(Nt-1)

t(m+1)=t(m)+(T/(Nt-1))*(1+0.1*randn);

end

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

x(1)=0;

for n=1:(Nx-1)

x(n+1)=x(n)+(a/(Nx-1))*(1+0.1*randn);

end

%Задаем левое граничное условие u(t,0)=-t^2

for m=1:Nt

y(m,1)=-t(m)^2;

end

%Задаем начальные данные u(0,x)=x^2

for n=1:Nx

y(1,n)=x(n)^2;

end

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

for m=1:(Nt-1)

tau=t(m+1)-t(m);

for n=2:Nx

c=t(m+1)/x(n); xa=x(n)-c*tau;

p=1;

while (~((xa>=x(p))&(xa<x(p+1))))&((p+1)~=Nx)

p=p+1;

end

if xa>0

%интерполяция

y(m+1,n)=((x(p+1)-xa)/(x(p+1)-x(p)))*...

y(m,p)+((xa-x(p))/(x(p+1)-x(p)))*y(m,p+1);

else

y(m+1,n)=-(t(m)-xa/c)^2;

end

end

end

%Находим различие между численным и

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

for m=1:Nt

for n=1:Nx

yerror(m,n)=abs(y(m,n)-ya(t(m),x(n)));

end

end

%Запоминаем ошибку численного решения в норме C

error_without_template(l)=max(max(yerror));

l=l+1;

end

%Рисуем зависимость ошибки численного решения от

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

plot(50:50:kmax,error_without_template,...

'-*','MarkerSize',12);

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

function y=ya(t,x)

y=x^2-t^2;

 

Знакопеременная скорость переноса c (t, x). В этом случае задача поставлена корректно, когда условия поставлены на тех границах, характеристики из которых идут внутрь области G (t, x).

Пусть скорость c (t, x) непрерывна в области G (t, x) = [0, T ]´[0, a ] и меняет знак на линии , т.е. , при этом будем считать , и , .

 

Рис.8. Зависимость ошибки численного решения задачи (19), (20)
от числа узлов сетки по времени и пространству
для схемы без шаблона

 

Для примера рассмотрим уравнение переноса следующего вида

(21)

где . Найдем вид характеристик для уравнения (21). Для этого необходимо решить обыкновенное дифференциальное уравнение

. (22)

Решение уравнения (22) легко найти. Решение содержит две ветви, из которых выбираем положительную (t > 0), т.е.

, A = const. (23)

Нарисуем характеристики (23) средствами MATLAB. На листинге_№6 приведен короткий код программы рисования характеристик (23). Итог работы кода программы листинга_№6 приведен на рис.9,а. В целях упрощения кода на рис.9,а нарисованы лишь те характеристики, которые выпущены из левой и правой границ области G (t, x) = [0, T ]´[0, a ]. Стрелкой на рис.9,а обозначена линия, на которой скорость переноса обращается в ноль.

 

Листинг_№6

%Программа рисования характеристик

%t=sqrt(A-ln(xv-x)^2) уравнения

%du/dt+t(xv-x)du/dx=0

%очищаем рабочее пространство

clear all

%задаем размеры области G=[0,T]x[0,a]

T=1; a=1;

%определяем координату xv, на которой

%скорость переноса меняет знак

xv=0.5;

t=0:0.1:T;

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

for m=1:length(t)

A(m)=t(m)^2+log(xv^2);

end

%Рисуем характеристики, выходящие из

%левой границы области G=[0,T]x[0,a]

x=0:0.001:(xv-1e-4);

for m=1:length(A)

for n=1:length(x)

y(n)=sqrt(A(m)-log((xv-x(n))^2));

end

plot(x,y);

hold on

end

%Рисуем характеристики, выходящие из

%правой границы области G=[0,T]x[0,a]

x=(xv-1e-4):0.001:a;

for m=1:length(A)

for n=1:length(x)

y(n)=sqrt(A(m)-log((xv-x(n))^2));

end

plot(x,y);

hold on

end

 

Рис.9,а. Характеристики (23) уравнения (21) Рис.9,б. Шаблон разностной схемы для уравнения (24)

 

Возвращаясь к вопросу о корректности уравнения переноса со скоростью переноса меняющей знак и учитывая пример (21) — (23), можно констатировать, что необходимо задавать два граничных условия на левой и правой границах области G (t, x) = [0, T ]´[0, a ], т.е.

(24)

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

Есть и другой способ. Построим для шаблона рис.9,б неявную разностную схему:

(25)

С учетом направления характеристик (стрелки на рис.9,а) и принимая во внимание шаблон на рис.9,б видно, что при любом знаке скорости c и при любом соотношении шагов по времени t и пространству h значение вычисляется интерполяцией. Методом разделения переменных можно убедиться, что схема устойчива при любом знаке c.

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

. (26)

Условие Куранта (26) для применимости метода прогонки является достаточным, поэтому при его нарушении схема (25) может все еще давать разумные численные решения. Исследуем этот вопрос на примере решения дачи (24) в виде:

(27)

Задача (27) имеет аналитическое решение

.

На листинге_№7 приведен код программы численного решения задачи (27) с помощью разностной схемы (25) для сеток разной длины по времени и пространству.

 

Листинг_№7

%Программа численного решения задачи (27)

%Очищаем рабочее пространство

clear all

%Определяем область интегрирования G=[0,T]x[0,a] и

%прямую x=xv, на которой скорость в уравнении

%переноса меняет знак

T=1; a=1; xv=0.5;

k=1;

%Организуем цикл решения задачи (27) для различных

%сеток: Nt - число узлов в сетке по времени,

%Nx - число узлов в сетке по пространству

for Nt=10:10:100

for Nx=10:10:100

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

tau=T/(Nt-1); h=a/(Nx-1);

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

for n=1:Nx

y(1,n)=-log((xv-h*(n-1))^2);

end

%программируем процедуру прогонки решения

%задачи (27) по разностной схеме (25)

for m=1:(Nt-1)

alpha(2)=0;

beta(2)=-(tau*m)^2-log(xv^2);

for n=2:(Nx-1)

c=tau*m*(xv-h*(n-1));

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

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

beta(n+1)=(y(m,n)+r*beta(n))/...

(1-r*alpha(n));

end

y(m+1,Nx)=-(tau*m)^2-log((xv-a)^2);

for n=Nx:-1:2

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

end

end

%оцениваем ошибку численного решения, как

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

%решения в норме l2

s=0;

for m=2:(Nt-1)

for n=2:(Nx-1)

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

log((xv-h*(n-1))^2))^2*tau*h;

end

end

%запоминаем норму ошибки и число узлов NtNx

%разностной схемы в области интегрирования G

error_l2(k)=sqrt(s);

gr(k)=Nt*Nx;

k=k+1;

end

end

t=0:tau:T;

x=0:h:a;

%Рисуем численное решение задачи (27) в области G

%с максимальным количеством узлов NtNx

subplot(1,2,1); surf(x,t,y);

%Рисуем график зависимости ошибки численного

%решения задачи (27) в норме l2 от произведения NtNx

subplot(1,2,2); plot(gr,error_l2);

 

На рис.10 приведен итог работы кода программы листинга_№7. На рис.10 слева изображено численное решение задачи (27) согласно разностной схеме (25) при выборе области интегрирования G (t, x) = [0,1]´[0,1]. По времени и пространству выбирались равномерные сетки: 0 = t 1 < … < t 100 = 1 и 0 = x 1 < … < x 100 = 1. На рис.10 справа приведен график зависимости разности численного и аналитического решений в норме от числа узлов сетки в области G (t, x) = [0,1]´[0,1]. Если число узлов по времени Nt, а по пространству — Nx, то общее число узлов сетки Nt ´ Nx. Из графика видна тенденция уменьшения ошибки в норме по мере роста числа узлов сетки Nt ´ Nx.

Вернемся теперь к условию Куранта в форме (26). Для параметров задачи (27), выбранных в листинге_№7, соотношение Куранта варьировалось в довольно широком диапазоне (), при этом проблем с устойчивостью численного решения методом прогонки не наблюдалось.

 

Рис.10.Численное решение задачи (27) (рисунок слева) и зависимость ошибки
численного решения в норме от числа узлов сетки в области
интегрирования G = [0, T ]´[0, a ] (рисунок справа)

 

Квазилинейное уравнение

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

, (28)

которое возникает при одномерном описании движения жидкости и газа. В уравнении переноса (28) скорость переноса определяется самим решением, т.е. c = u (t, x).

В дальнейшем будем рассматривать простейший случай уравнения (28), когда решение u знакопостоянно u (t, x) > 0 и граничное и начальное значения на положительных полуосях системы координат (t, x) полностью определяют решение в первом квадранте. Поскольку граничное и начальное значения передаются по характеристикам x - ut = const, а их наклон зависит от решения, постольку возможны такие начальные данные, при которых характеристики могут пересекаться, что приводит к появлению разрывных решений. С этой точки зрения существует четыре типа начальных данных.

Первый тип начальных данных. Начальные и граничные являются непрерывными функциями такими, что u (0, x) монотонно не убывает, а u (t,0) монотонно не возрастает и они непрерывно согласованы в начале координат.

Наклон (тангенс угла наклона к оси x) характеристик в каждой точке плоскости (t, x) равен 1/ u (t, x). При данном типе начальных данных наклон монотонно и непрерывно убывает слева направо (рис.11,а). Первый квадрант всюду плотно покрыт характеристиками и через каждую точку проходит одна и только одна характеристика. По этим характеристикам в данную точку переносится либо граничное, либо начальное значения. Решение однозначно определено и непрерывно во всем квадранте.

Для примера возьмем следующие граничные и начальные условия

,

тогда характеристики, выходящие из точек (t,0), t Î [0, T ] и (0, x), x Î [0, a ], в окрестности этих точек удовлетворяют уравнениям:

. (29)

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

, (30)

где t 0, x 0 — константы, пробегающие значения по отрезку [0, T ] и отрезку [0, a ] соответственно.

На листинге_№8 приведен код небольшой программы, которая изображает приближенные характеристики (30) в области G = [0, T ]´[0, a ]. Итог работы кода программы листинга_№8 приведен на рис.11,а. Стрелка на рис.11,а указывает общее направление переноса граничных и начальных условий.

 

Листинг_№8

%Программа рисования приближенных

%характеристик (30)

%очищаем рабочее пространство

clear all

%задаем размеры области G=[0,T]x[0,a] и

%константу alpha

T=1; a=1; alpha=0.3;

%определяем набор значений констант t0 и x0

t0=0:0.1:T; x0=0:0.1:a;

t=0:0.1:T; x=0:0.1:a;

%рисуем характеристики, выходящие из

%левого граничного условия

for m=1:length(t0)

for n=1:length(x)

y(n)=(exp(alpha*x(n))-1)/alpha+...

t0(m)*exp(alpha*x(n));

end

plot(x,y);

hold on

end

%рисуем характеристики, выходящие из

%начального условия

for n=1:(length(x0)-1)

x=x0(n):((a-x0(n))/10):a;

for nn=1:length(x)

y(nn)=log((1+alpha*x(nn))/...

(1+alpha*x0(n)))/alpha;

end

plot(x,y);

hold on

end

 

Рис.11,а. Характеристики первого типа начальных данных Рис.11,б. Характеристики второго типа начальных данных с пустой областью Рис.11,в. Характеристики второго типа начальных данных с доопределенной пустой областью

 

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

(31)

На рис.11,б приведены соответствующие характеристики. В силу (31) характеристики приближенно можно представить прямыми линиями, которые имеют одинаковый наклон левее и правее точки разрыва x 0. Из точки разрыва начальных данных x 0 выходит две характеристики, выделенные на рис.11,б жирными стрелками красного цвета. Между этими линиями нет ни одной характеристики и решение в этой области не определено.

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

(32)

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

Разрыв производных называют слабым разрывом решения. Слабые разрывы квазилинейного уравнения переносятся по характеристикам.

Третий тип начальных данных. Пусть нарушено условие монотонности, предполагаемое в двух предыдущих случаях. Положим, как и в (31),

но теперь a > b > 0. В этом случае характеристики будут иметь вид, представленный на рис.12,а.

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

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

 

Рис.12,а. Характеристики третьего типа начальных данных Рис.12,б. Построение обобщенного, разрывного решения Рис.12,в. Обобщенное, разрывное решение и поле характеристик

 

Квазилинейное уравнение переноса (28) легко переписать в одной из дивергентных форм:

. (33)

Нас интересует решение, имеющее единственный разрыв. Пусть наклон линии разрыва соответствует скорости V, и разрыв двигается как волна. Из вида характеристик на рис.12,а следует, что искомое разрывное решение должно иметь следующий вид:

Проинтегрируем дивергентную форму (33) по площади прямоугольника со сторонами t и h = Vt, и левый нижний угол которого находится в точке разрыва x 0 (рис.12,б), тогда

откуда находим скорость распространения разрыва

. (34)

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

Четвертый тип начальных данных. В этом случае функция начального условия u (0, x) непрерывна и убывает на некотором интервале, что переадресует нас к третьему типу. Пересечение характеристик также приводит к образованию сильного разрыва (рис.12,в), скорость которого описывается уравнением, подобным (34). Однако в этом случае скорость сильного разрыва может меняться.

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

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

Рассмотрим пример. Добавим в исходное квазилинейное уравнение (28) соответствующий член, тогда

, (35)

где третье слагаемое выступает в качестве псевдовязкости, а коэффициент e 2 считается малым.

Видно, что для гладких, дважды непрерывно-дифференцируемых функций решения уравнения (28) близки к решениям уравнения (35).

Поищем среди гладких решений уравнения (35) решение, напоминающее ударную волну:

где b < a и скорость ударной волны . Рассмотрим автомодельное решение в виде бегущей волны

. (36)

После подстановки (36) в (35), находим

. (37)

Приравнивая каждый из сомножителей в (37) к нулю, получаем два решения:

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

(38)

Решение (38) непрерывно-дифференцируемо и при e ® 0 сходится к разрывному решению — ударной волне. На листинге_№9 приведен код программы рисования последовательности профилей (38), приводящих к образованию ударной волны при e ® 0. Итоговый рисунок приведен на рис.13.

 

Листинг_№9

%Программа рисования по формуле (38)

%последовательности формирования

%ударной волны при eps -> 0

%очищаем рабочее пространство

clear all

%задаем параметры функции (38)

a=2; b=0.2; eps=2;

%определяем переменную xi=x-Vt

xi=-2:0.01:2;

%Формируем цикл построения графиков

%ueps при различных значениях eps

for i=1:10

%уменьшаем текущее значение

%значение eps вдвое

eps=eps/2;

for j=1:length(xi)

if xi(j)<-0.5*pi*eps

ueps(j)=a;

end

if (xi(j)>=(-0.5*pi*eps))&...

(xi(j)<=(0.5*pi*eps))

ueps(j)=0.5*(a+b)-...

0.5*(a-b)*sin(xi(j)/eps);

end

if xi(j)>0.5*pi*eps

ueps(j)=b;

end

end

axis([-2 2 0 2]);

plot(xi,ueps,'LineWidth',2);

hold on

end

 

Обычно коэффициент псевдовязкости связывается с шагом сетки по пространству, т.е. считают, что

e = n h, n = const, (39)

тогда любой сильный разрыв согласно (38) размазывается на одно и то же число шагов сетки pe / h = pn. В этом случае при h ® 0 уравнение с псевдовязкостью переходит в исходное квазилинейное уравнение (28), а размазанная ударная волна (38) переходит в разрывную ударную волну.

Составим для уравнения с псевдовязкостью (35) явную разностную схему на равномерной сетке:

(40)

 

Рис.13. Образование разрывного решения — ударной
волны при e ® 0 согласно представлению (38)

 

Исследуем вопрос об устойчивости разностной схемы (40). Поскольку разностная схема (40) является нелинейной, линеаризуем ее, определяя уравнение для погрешности dy:

(41)

Поскольку коэффициенты при dy являются переменными, для дальнейшего исследования вопроса устойчивости схемы (40) воспользуемся принципом “замороженных” коэффициентов. Согласно этому принципу коэффициенты при dy считаются постоянными величинами. Проведем в (41) замены и т.д., тогда

(42)

Изучим рост ошибки, которая имеет вид . Подставим в (42) выражения

,

тогда множитель роста гармоники предстанет в виде:

. (43)

Если коэффициент псевдовязкости выбран согласно (39), то величина в квадратных скобках в (43) равномерно ограничена по шагу h. В этом случае последний член в (43) имеет порядок O (t) и не влияет на устойчивость разностной схемы (40). Первые два члена в правой части уравнения (43) для обеспечения условия устойчивости разностной схемы (40) приводят к условию, аналогичному условию Куранта:

.

Схема (40) является примером так называемой однородной разностной схемы для решения задач с произвольным числом движущихся разрывов, число которых может меняться со временем. Отметим, что помимо псевдовязкости уравнения (35), применяют также псевдовязкость, именуемую линейной. Она имеет следующий вид:

, (44)

где , n 1 = const. Уравнение (44) напоминает уравнение теплопроводности все решения которого достаточно гладкие.

Для уравнения (44) определим явную разностную схему вида:

. (45)

Проведем сравнительный численный анализ схем (40) и (45). На листинге_№10 приведен код соответствующей программы расчета образования ударной волны по схемам (40) и (45).

 

Листинг_№10

%Программа сравнения разностных схем (40), (45)

%на примере динамики образования ударной волны

%Очищаем рабочее пространство

clear all

%Задаем пределы области интегрирования G=[0,T]x[0,a]

T=4; a=2;

%Задаем параметры сетки по времени и пространству

tau=0.001; h=0.01; r=tau/h;

%Определяем параметры псевдовязкости схем (40), (45)

nu=1; nu1=0.1;

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

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

Nt=length(t); Nx=length(x);

%Определяем левое граничное условие u(t,0)=0

for m=1:Nt

y(m,1)=0;

end

%Определяем правое граничное условие u(t,0)=0

for m=1:Nt

y(m,Nx)=0;

end

%Определяем начальное условие u(0,x)=sin(pi x)^2,

%0<=x<=1; u(0,x)=0, 1<x<=2

for n=1:Nx

if x(n)<=1

y(1,n)=sin(pi*x(n))^2;

else

y(1,n)=0;

end

end

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

%моменты времени

for m=1:(Nt-1)

for n=2:(Nx-1)

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

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

(y(m,n+1)-2*y(m,n)+y(m,n-1));

end

end

%Рисуем профили волны в различные моменты времени на

%одном и том же графике в координатах решение -

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

for m=1:35:Nt

subplot(1,2,1);

plot(x,y(m,:));

hold on

end

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

%моменты времени

for m=1:(Nt-1)

for n=2:(Nx-1)

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

r*nu1*(y(m,n+1)-2*y(m,n)+y(m,n-1));

end

end

%Рисуем профили волны в различные моменты времени

for m=1:35:Nt

subplot(1,2,2);

plot(x,y(m,:));

hold on

end

 

Итоговый результат работы кода листинга_№10 представлен на рис.14. Слева на рис.14 приведены профили волны на различные моменты времени, рассчитанные по разностной схеме (40). Жирными линиями на графике выделено начальное распределение и уже сформировавшаяся ударная волна в последний момент времени. Заметно, что ударная волна слегка размазана в соответствии с выбором величины параметра n. Аналогичная картина представлена на рис.14 справа, где приведены результаты расчета по “линейной” разностной схеме (45). Визуальное сравнение левого и правого рисунков показывает их неплохое совпадение. Это удалось добиться за сче



Поделиться:




Поиск по сайту

©2015-2024 poisk-ru.ru
Все права принадлежать их авторам. Данный сайт не претендует на авторства, а предоставляет бесплатное использование.
Дата создания страницы: 2016-04-11 Нарушение авторских прав и Нарушение персональных данных


Поиск по сайту: