и единственную на данном шаблоне схему третьего порядка аппроксимации с
порядком :
( 3.16) |
Остальные коэффициенты находятся с использованием (3.15) и (3.16).
Введем пространство коэффициентов . Тогда
любая точка в этом пространстве есть разностная схема с порядком аппроксимации . Прямая (3.15) отделяет в пространстве множество схем с порядком
(рис.
3.2), на ней лежит единственная точка — аппроксимация . Должна быть также и точка с порядком аппроксимации
.
Зафиксируем какое — либо число Куранта, например, .
Применим к разностной схеме (3.14) с неопределенными коэффициентами
спектральный признак устойчивости (фон Неймана). Получится кривая, которая определяет границу устойчивости разностных схем в пространстве неопределенных коэффициентов.
Рис.
3.2.
Для схем первого порядка выпишем первое дифференциальное приближение:
Можно также выделить множество таких схем, что .
Это — монотонные схемы (заштрихованный многоугольник на рис. 3.2). Среди монотонных схем можно найти схему с наименьшей ошибкой аппроксимации. Это точка многоугольника, которая при данном лежит ближе всего к прямой со схемами второго порядка аппроксимации.
Закончим рассмотрение примера с модельным линейным уравнением переноса:
На выбранном шаблоне любая разностная схема, как указывалось ранее,
представляется в виде
В случае монотонной схемы можно оценить норму погрешности. Заметим, что
погрешность v определяется тем же разностным уравнением (3.14),
тогда с использованием первой нормы (максимум абсолютной величины)
в силу аппроксимации
Отсюда следует, что монотонные разностные схемы всегда устойчивы. В общем
случае можно рассматривать многослойные шаблоны для уравнения переноса (рис. 3.3)
и записывать условия порядка для аппроксимации соответствующего порядка:
Рис.
3.3.
Исключая два коэффициента из условий порядка, можно от пространства неопределенных коэффициентов перейти к пространству
, размерность
которого на 2 меньше, например:
где, конечно, точки (0; — 1) и (0; 1) не включаются в
суммирование. Условие устойчивости в пространстве неопределенных коэффициентов имеет вид
где q есть спектр оператора послойного перехода. Эта величина
определяется из условия
и дополнительного требования
Основная гипотеза: Разностным схемам, которым в пространстве соответствуют близкие друг к другу точки (в смысле
) по своим свойствам также близки.
Расширяя шаблон, как и в случае пятиточечного шаблона, можно строить области схем высокого порядка аппроксимации, монотонные схемы ( ) и т.д.
Монотонные разностные схемы в пространстве неопределенных коэффициентов
занимают некий выпуклый многоугольник, вершины которого определяются довольно просто: это все возможные при данном числе Куранта трехточечные разностные схемы, причем для характеристики, проходящей через
, одна точка схемы лежит выше (левее) характеристики, а другая ниже (правее), см. рис. 3.3.
Метод построения разностных схем в пространстве неопределенных коэффициентов для квазилинейных систем уравнений гиперболического типа (к ним относятся системы уравнений механики сплошных сред, в частности, газовой динамики, механики деформируемого твердого тела (МДТТ) и т.п.) допускает обобщение и на многомерные случаи. Подробное описание можно найти в монографии [13.9]. Здесь же многомерные обобщения рассматриваться не будут. Они приводят к эффективным численным методам для нестационарных многомерных задач.
Исследовав схемы (3.14) на устойчивость по спектральному признаку, получаем множество устойчивых схем, а потребовав выполнения условия для всех точек шаблона, получаем множество схем с положительной аппроксимацией (монотонных по Фридрихсу схем). На рассматриваемом
шаблоне устойчивые схемы существуют при (для a > 0 ).
Множество схем с положительной аппроксимацией не пересекается с множеством
схем с порядком аппроксимации выше первого, как это следует из теоремы С.К. Годунова.
Первое дифференциальное приближение. Дисперсионная и диссипативная
ошибки
Поскольку решения дифференциальной задачи и разностного уравнения
принадлежат разным функциональным пространствам, что порождает определенные трудности при теоретическом анализе свойств разностных схем, для такого исследования возможно рассматривать разностные операторы в том же пространстве. Будем считать, что разностные схемы удовлетворяются функциями непрерывного аргумента в каждой точке рассматриваемой области.
Обычно ограничиваются рассмотрением уравнений, в которых оставлены члены в
разложении в ряд Тейлора проекции точного решения на сетку по и h, порядок которых совпадает с порядком погрешности аппроксимации схемы. Получающиеся при этом уравнения называют первым дифференциальным
приближением (ПДП).
Для схемы первого порядка (5) при выполнении условия (6) первым дифференциальным приближением будет
( 3.17) |
Из уравнения (3.17) исключены члены со второй производной u’tt, с использованием так называемой продолженной системы:
Иногда уравнение (3.17) называют — формой (параболической формой) первого дифференциального приближения. Если производные по времени не исключаются из ПДП, то имеем Г — форму (гиперболическую
форму) ПДП, которая, как правило, не применяется в исследованиях как малоинформативная.
При уравнение (3.17) можно трактовать как
присутствие в схеме некоторой диссипации ( схемной вязкости ). Ее наличие проявляется в расчетах в виде размазывания точного решения, причем его интенсивность увеличивается при ухудшении аппроксимации (увеличении шага h ). В этом случае говорят, что ошибка схемы носит диссипативный характер. Если схемная вязкость получается отрицательной, то приходим к обратной задаче теплопроводности. Как известно из курсов математической физики, такая задача поставлена некорректно. А соответствующая разностная схема при исследовании по спектральному признаку оказывается неустойчивой — по ПДП можно сделать вывод об устойчивости схемы.
Для более высокого (второго) порядка ПДП имеет вид
( 3.18) |
Уравнение (3.18) обладает дисперсией, т.е. разные пространственные гармоники разложения начального возмущения в ряд Фурье распространяются по сетке с разными скоростями. Говорят, что ошибка носит дисперсионный характер. Сеточная дисперсия легко получается, если искать частное решение последнего уравнения в виде комплексной экспоненты: . Подробнее о ПДП в [13.4].
Интересна связь ПДП и исследования свойств схем в пространствах неопределенных коэффициентов. Так, для схем первого порядка расстояние от точки в пространстве неопределенных коэффициентов до прямой схем высокого порядка (3.15) по абсолютной
величине равен коэффициенту в уравнении (3.17), а знак определяется положением внутри области устойчивости.
Понятие о гибридных схемах
Численные расчеты по разностным схемам высокого порядка показывают, что
осцилляции нефизического характера появляются в окрестности разрывов
решения или его первой производной. В связи с этим возникает идея
построения численного метода, имеющего высокий порядок аппроксимации на
участках гладкого решения, в то время как в окрестности разрывов функций
или их производных применяется монотонная схема (с положительной
аппроксимацией) первого порядка.
Можно формализовать этот подход:
где — коэффициенты первой схемы (высокого
порядка аппроксимации), применяемой в области гладкого решения, — коэффициенты второй (монотонной) схемы,
— весовой коэффициент, вспомогательный параметр b характеризует гладкость решения (очевидно, что b = 0, при
), k — коэффициент гибридности — целое число из диапазона
.
При этом реализовано достаточно гладкое переключение со схем высокого порядка аппроксимации на монотонные схемы.
|
|
#1 |
New Member Join Date: Apr 2013 Posts: 15 Rep Power: 12 |
Hi all, We know that for the spatial derivative using a finite difference method, if the leading error term is odd, then the results contain dispersive error, while if the leading error term is even, then it has dissipative error. I want to compare the degree of such errors. Say, the leading error term is 7th order (FD7), compared with the leading error term is 5th order (FD5), which one has more severe dispersive error? For sure, FD7 is more accurate, what’s its performance with regard to dispersive error? Thanks. Shu |
|
|
|
|
|
#2 |
Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,595 Rep Power: 70 |
Quote:
Originally Posted by shubiaohewan Hi all, We know that for the spatial derivative using a finite difference method, if the leading error term is odd, then the results contain dispersive error, while if the leading error term is even, then it has dissipative error. I want to compare the degree of such errors. Say, the leading error term is 7th order (FD7), compared with the leading error term is 5th order (FD5), which one has more severe dispersive error? For sure, FD7 is more accurate, what’s its performance with regard to dispersive error? Thanks. Shu This is a classical numerical analysis task, you can analyze the local truncation error of the discretization formula, it gives much information about the character of the error. |
|
|
|
|
|
#3 |
Senior Member cfdnewbie Join Date: Mar 2010 Posts: 557 Rep Power: 19 |
There is a standard procedure of analysing dispersive and dissipative behavior of FD schemes. Discretize a linear advection equation, plug in a wave with frequency k and amplitude a and check the resulting frequency and amplitude response — that gives you the dissipation and dispersion error of your scheme. |
|
|
|
|
|
#4 |
Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,595 Rep Power: 70 |
Quote:
Originally Posted by cfdnewbie There is a standard procedure of analysing dispersive and dissipative behavior of FD schemes. Discretize a linear advection equation, plug in a wave with frequency k and amplitude a and check the resulting frequency and amplitude response — that gives you the dissipation and dispersion error of your scheme. Such analysis can be extended, in a suitable way, also for the non linear equation (e.g., Burgers) |
|
|
|
|
|
#5 |
Senior Member cfdnewbie Join Date: Mar 2010 Posts: 557 Rep Power: 19 |
Hello Prof. Denaro, |
|
|
|
|
|
#6 |
Senior Member
Paolo Lampitella Join Date: Mar 2009 Location: Italy Posts: 2,120 Blog Entries: 29 Rep Power: 38 |
I don’t know if Prof. Denaro is referencing what i have in mind; however, a possible practical (i.e., non analytical) approach is this: http://www.sciencedirect.com/science…21999111001148 |
|
|
|
|
|
#7 |
Senior Member cfdnewbie Join Date: Mar 2010 Posts: 557 Rep Power: 19 |
Thank you Paolo, |
|
|
|
← Назад Далее →
Все построенные на равномерной сетке разностные операторы являются однородными — их коэффициенты в вычислительных шаблонах инвариантны относительно сдвигов шаблона по сетке. Для однородных дискретных операторов можно построить частные решения в виде бегущих волн:
$$u_{j}^{n} =D\cdot \exp \left\{i\left[{\rm\omega }\cdot n{\rm\tau }-kjh\right]\right\}=A\cdot q^{n} \cdot \exp \left\{-i\left(kjh\right)\right\};{\rm\; \; \; } (2.38)$$
здесь \({\rm\; \; \; }q=\exp \left\{i\left({\rm\omega \tau }\right)\right\}\) — множитель переходана следующий временной слой.
Подставим (2.38) в разностный оператор схемы «уголок» (2.26). Получаем:
$$\begin{array}{l} {q=r\cdot e^{ikh} +\left(1-r\right)\cdot ,\; \Rightarrow \; } \\{\Rightarrow \left|q\left(r,kh\right)\right|{\rm =}\sqrt{\left[\left(1- r\right)+r\cdot \cos \left(kh\right)\right]^{2} {\rm +}\left[{\rm sin}\left(kh\right)\right]^{2} } {\rm\; \; }} \end{array} (2.39)$$
Соотношение, выражающее множитель перехода через параметры дискретного оператора, будем называть характеристическим уравнениемразностной схемы.
Модуль множителя перехода \(\left|q\right|\) является функцией от числа Куранта — Фридрихса-Леви \(r\) и приведенного волнового числа \(kh\), изменяющегося в пределах от \(-{\rm\pi }\): \(\left(-{\rm\pi }\le kh\le {\rm\pi }\right)\). Он показывает, как изменяется амплитуда бегущей волны при переходе на следующий слой. Если модуль перехода меньше или равен единице, то амплитуда волны не возрастает. Если для заданного \(r\) амплитуды всех волн при \(\left(-{\rm\pi }\le kh\le {\rm\pi }\right)\) не возрастают, то говорят, что данное число \(r\) принадлежит области устойчивости дискретного оператора. Опираясь на (2.39) можно показать, что область устойчивости разностной схемы «уголок» представляет собой единичный отрезок \(r\in \left[0,1\right]\). Если в области устойчивости модуль перехода любой волны равен единице, то дискретный оператор называется бездиссипативным, в противном случае имеет место численная диссипация.
Диссипативные свойства дискретного оператора для всех приведенных волновых чисел наиболее наглядно можно представить в виде двумерной поверхности \(\left|q\right|=f\left(r,kh\right)\). Для схемы «уголок» диссипативная поверхность изображена на рис. 6а, где по оси \(z\) отложены значения модуля перехода, по оси \(x\) — значения числа \(r\in \left[0,1\right]\), по оси \(y\) — значения приведенного волнового числа \(kh\in \left[-{\rm\pi },{\rm\pi }\right]\). Видно, что модуль перехода стремится к единице для длинноволновых гармоник, для которых \(kh\approx 0\), и равен единице на границе области устойчивости, при \(r=1,{\rm\; }\left(-{\rm\pi }\le kh\le {\rm\pi }\right)\).
Из (2.39) также следует дисперсионное соотношение:
$$q=e^{i{\rm\omega \tau }} =r\cdot e^{ikh} +\left(1-r\right)\cdot ,\; \Rightarrow \; {\rm\omega }=-\frac{i}{{\rm\tau }} \arg \left[r\cdot e^{ikh} +\left(1-r\right)\right]$$
откуда получается выражение для приведенной фазовой скорости бегущих волн:
$$\begin{array}{l} {\gamma \left(r,kh\right){\rm =}\frac{{\rm\omega }}{kc} {\rm\; }=-\frac{i}{{\rm\tau }kc} \arg \left[r\cdot e^{ikh} +\left(1-r\right)\right]{\rm =}} \\{{\rm =}-\frac{i}{r\cdot kh} \arg \left[r\cdot e^{ikh} +\left(1-r\right)\right]{\rm\; \; }} \end{array}$$
Приведенная фазовая скорость бегущих волн, описываемых исходным дифференциальным оператором, не зависит от номера гармоники и равна единице. Как ранее уже отмечалось, зависимость фазовой скорости от длины волны называется дисперсией. В разностном случае приведенная фазовая скорость может принимать значения как меньшие, так и большие единицы. В первом случае дисперсия называется «нормальной», во втором — «аномальной». В случае нормальной дисперсии гармоника отстает от «эталонной» гармоники, двигающейся с фазовой скоростью \(c\), в случае аномальной дисперсии — опережает ее.
Дисперсионные свойства дискретного оператора представляются его «дисперсионной поверхностью». На рис. 6в приведена дисперсионная поверхность разностного оператора «уголок». Здесь по оси \(z\) отложены значения приведенной фазовой скорости \({\rm\gamma }\), по оси \(x\)- значения числа \(r\in \left[0,1\right]\), по оси \(y\) — значения приведенного волнового числа \(kh\in \left[-{\rm\pi },{\rm\pi }\right]\). Дисперсия стремится к нулю, что соответствует стремлению \({\rm\gamma }\) к единице, для длинноволновых гармоник, для которых \(kh\approx 0\) и на границе области устойчивости, при \(r=1,{\rm\; }\left(-{\rm\pi }\le kh\le {\rm\pi }\right)\).
Стремление к нулю дисперсии и диссипации длинноволновых гармоник является общим свойством всех дискретных операторов, аппроксимирующих исходный дифференциальный оператор. Характер касания дисперсионных и диссипативных поверхностей единичной плоскости определяется порядком аппроксимации разностного оператора. Тот факт, что при \(r = 1\)все гармоники не диссипируют и не диспергируют, говорит о том, что разностный оператор при данном числе Куранта переносит начальный профиль без каких либо искажений.
Из непрерывности дисперсионной и диссипативной поверхностей следует, что в некоторой окрестности этого числа точность численного решения также будет аномально высокой и эту область можно назвать «каналом высокой точности» дискретного оператора.
Из рис. 6в видно, что дисперсия схемы «уголок» при \(r < 0.5\) является нормальной, а при \(r > 0.5\) — аномальной. При \(r = 0.5\)дисперсия отсутствует.
Диссипативная и дисперсионная поверхности дают полное представление о свойствах разностной схемы, позволяют предсказывать поведение дискретного решения на разных начальных условиях. Рассмотрим модельную задачу (2.1) с периодическими начальными условиями на отрезке \(x\in \left[0,2{\rm\pi }\right]\):
$$\begin{array}{l} {\frac{\partial u}{\partial t} +c\cdot \frac{\partial u}{\partial x} =0,{\rm\; \; }c=1;{\rm\; \; }t\in \left[0,\infty \right);} \\{u\left(x,0\right)=f\left(x\right);{\rm\; \; \; }u\left(0,t\right)=u\left(1,t\right);x\in \left[0,2{\rm\pi }\right];} \end{array} (2.40)$$
и двумя видами начальных условий — достаточно гладкими и разрывными.
В качестве гладкого начального условия возьмем «двойной гауссиан» с параметрами:
$$\begin{array}{l} {f\left(x\right)=\exp \left[-\left(x-x_{1} \right)^{2} /\Delta ^{2} \right]+\exp \left[-\left(x-x_{2} \right)^{2} /\Delta ^{2} \right];} \\{x_{1} =0.5\cdot {\rm\pi },{\rm\; \; \; }x_{2} =0.7\cdot {\rm\pi },{\rm\; \; \; }\Delta =0.1\cdot {\rm\pi }} \end{array} (2.41)$$
в качестве разрывного — «ступеньку»:
$$f\left(x\right)=\left\{\begin{array}{l} {1,{\rm\; \; \; \; if\; \; \; \; \; \; }0.2\cdot {\rm\pi }\le x\le 0.4\cdot {\rm\pi }} \\{0,{\rm\; \; \; \; if\; \; \; \; \; \; }\left(0\le x \le 0.2\cdot {\rm\pi }\right)\& \left(0.4\cdot {\rm\pi } \le x \le 2{\rm\pi }\right)} \end{array}\right. (2.42)$$
Аппроксимируем дифференциальную задачу (2.40) — (2.42) разностной схемой «уголок». Разобьём отрезок \(x\in \left[0,2{\rm\pi }\right]\) на N=100 равных интервалов с шагом \(h=0.2{\rm\pi }\) и спроектируем на нее начальные данные (2.41) — (2.42):
$$\begin{array}{l} {f_{i}^{\left(1\right)} =\exp \left\{-\left[\left(i-1\right)h-x_{1} \right]^{2} /\left(0.1\cdot {\rm\pi }\right)^{2} \right\}+} \\{+\exp \left\{-\left[\left(i-1\right)h-x_{2} \right]^{2} /\left(0.1\cdot {\rm\pi }\right)^{2} \right\};} \\{f_{i}^{\left(2\right)} =\left\{\begin{array}{l} {1,{\rm\; \; \; \; if\; \; \; \; \; \; }11\le i\le 21} \\{0,{\rm\; \; \; \; if\; \; \; \; \; \; }\left(1\le i \le 11\right)\& \left(21 \le i\le 101\right)} \end{array}\right. ;{\rm\; }\quad i=1,…,101} \end{array} (2.43)$$
Разложим эти сеточные функции в ряд Фурье:
$$f_{n}^{\left(m\right)} =\sum _{k=0}^{N}a_{k}^{\left(m\right)} \cdot \exp \left[-i\left(kh\right)\cdot n\right];{\rm\; \; }m=1,2;{\rm\; \; }n=1,…,N+1 (2.44)$$
где коэффициенты разложения определяются по формуле:
$$a_{k}^{\left(m\right)} =\frac{1}{N+1} \cdot \sum _{n=1}^{N+1}f_{n}^{\left(m\right)} \cdot \exp \left[i\left(kh\right)\cdot n\right] (2.45)$$
На рис. 7 представлены графики начальных данных (2.43), на рис. 8 — модули их коэффициентов разложения \(a_{k}^{\left(m\right)}\).
В начальный момент времени все гармоники разложения (2.44) подогнаны друг к другу так, что в сумме составляют начальные данные. В последующие моменты времени каждая гармоника будет уменьшаться в амплитуде в соответствии с ее местом на диссипативной поверхности, и распространяться со своей фазовой скоростью, определяемой дисперсионной поверхностью, так что «пакет» гармоник начнет разъезжаться. Полная картина будет еще зависеть от числа Куранта-Фридрихса-Леви \(r\).
Из рис. 6 видно, что наибольшая диссипация высоких гармоник имеет место при \(r = 0.5\). На рис. 9 приведен результаты расчета переноса начальных профилей (2.43), рассчитанные по схеме «уголок» при разных числах Куранта, до одного и того же момента времени. Сплошной линией показано аналитическое решение дифференциальной задачи, линией с маркерами — решение, рассчитанное по разностной схеме.
← Назад Далее →