MKR_05

ГЛАВА: Метод конечных разностей.

Лекция №5: Устойчивость разностных схем (10 слайдов, 6 рисунков)

Слайд №1: Классификация РС по типам устойчивости.

По типам устойчивости выделяют следующие РС:

•абсолютно (безусловно) устойчивые;

•абсолютно (безусловно) неустойчивые;

•условно устойчивые.

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

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

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

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

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

Для анализа устойчивости РС разработано несколько методов. Наиболее известными из них являются:

•Метод дискретного возмущения;

•Метод Неймана.

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

1

Слайд №2: Динамическая неустойчивость.

Рассмотрим РС ВВЦП для уравнения теплопроводности

возмущенное решение Uin, можем записать

Члены, относящиеся к истинному решению с точностью до величины погрешности аппроксимации сокращаются (погрешность аппроксимации обуславливает лишь появление погрешности εni , но не ее динамику). Оставшиеся члены показывают, как поведет себя при переходе к следующему временному слою возмущение εni

εi = εni +1 − εni = aτh2 ³εni+1 − 2εni + εni−1´ .

Поскольку εni+1 > 0 и εni−1 > 0, а в i-ом узле εni < 0, то поправка положительна. Аналогично можно показать, что для узлов с положительным отклонением от истинного решения поправки будут отрицательны. Поправки корректируют предыдущие значения, уменьшая погрешность. Однако с увеличением шага по времени τ поправки εi растут. При некотором шаге τ поправка оказывается чрезмерной, и вместо коррекции наблюдается смена знака у начального возмущения (на рисунке показано стрелками), а при еще несколько большем шаге возмущения εi начинают не просто менять знак, но еще и усиливаться при переходе с одного временного слоя на следующий.

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

Явная схема ВВЦП для уравнения теплопроводности условно устойчива.

Явные РС, как правило, условно устойчивы. Исключение — РС Дюфорта–Франкела.

2

Слайд №3: Статическая неустойчивость.

возмущенное решение Uin, можем записать

Члены, относящиеся к истинному решению с точностью до величины погрешности аппроксимации сокращаются (погрешность аппроксимации обуславливает лишь появление погрешности εni , но не ее динамику). Оставшиеся члены показывают, как поведет себя при переходе к следующему временному слою возмущение εni

εi = εni +1 − εni = −2vhτ ³εni+1 − εni−1´ .

Если v > 0 (поток идет слева направо) и если возмущения растут в том же направлении: εi+1 > εi−1, то поправка отрицательна для всех узлов, т.е. начальное возмущение будет монотонно сдвигаться вниз. Это — статическая неустойчивость. Причем эту неустойчивость невозможно уничтожить, измельчая шаги сетку по переменным τ и h.

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

3

Слайд №4: Метод дискретного возмущения.

При сформулированной маршевой задаче РС можно интерпретировать как некоторый

оператор ˆ который используя значения искомой функции на некоторых слоях по мар

M, -

шевой координате (Uik, k ≤ n + 1) определяет значения функции на очередном слое (Uin+1):

где εki — погрешность. Подстановка последнего выражения позволяет определить оператор

вообще говоря отличный от ˆ описывающий динамику роста погрешности

( , M),

{ n+1} ˆ ³{ n+1} { n+1} { n} { n} { n−1} { n−1} ´

εi = N Ui , εi , Ui , εi , Ui , εi , . . . .

Далее, в некоторой произвольной точке (i, n) вводится дискретное возмущение величины ε, и анализируется влияние этого возмущения на последующие слои. РС устойчива, если

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

4

Слайд №5: Метод дискретного возмущения: схема ВВЦП.

Динамика погрешности определяется оператором

εni +1 = λεni+1 + (1 − 2λ)εni + λεni−1 , λ = aτ/h2 .

Вводим дисктетное возмущение величины ε в точке (i, n). Для (n + 1)-го слоя имеем

εni−+11 = λε , εni +1 = (1 − 2λ)ε , εni+1+1 = λε.

Нетрудно убедиться, что условие устойчивости (|εki | ≤ ε) на первом шаге сводится к:

λ ≤ 1 .

Но схема устойчива, если устойчивость выполняется на любом шаге по времени. Для (n + 2)-го слоя имеем

Это ужесточает условие устойчивости

λ ≤ 23 .

После этого надо проанализировать, что дают все последующие шаги по времени, т.е.

Анализ решений, получающихся неравенств дает в пределе n → ∞ условие: λ ≤ 1/2. Т.е. для того, чтобы РС ВВЦП для уравнения теплопроводности была устойчива, требуется выполнение определенного соотношения между параметром a в уравнении и шагами сетки

h2 τ ≤ τc = 2a .

τc — ”критический” размер шага по времени, который превращает РС ВВЦП из устойчивой (τ ≤ τc) в неустойчивую (τ > τc).

5

Слайд №6: Метод Неймана.

При сформулированной маршевой задаче РС можно интерпретировать как некоторый

оператор ˆ который используя значения искомой функции на некоторых слоях по мар

M, -

шевой координате (Uik, k ≤ n + 1) определяет значения функции на очередном слое (Uin+1):

{ n+1} ˆ ³{ n+1} { n} { n−1} { n−2} ´

Ui = M Ui , Ui , Ui , Ui , . . . .

где εki — погрешность. Подстановка последнего выражения позволяет определить оператор

вообще говоря отличный от ˆ описывающий динамику роста погрешности

( , M),

{ n+1} ˆ ³{ n+1} { n} { n−1} ´

εi = L εi , εi , εi , . . . .

Для линейных ДУ (а именно для них строго обоснован метод Неймана) получающийся

оператор ˆ линейный и не содержит величин k

L — Ui .

Какой бы ни была зависимость погрешности ε(x) на интервале интегрирования (0, L), ее всегда можно протранслировать за этот интервал, как показано на рисунке и считать

функцию ε(x) периодической. Периодическую функцию можно разложить в ряд Фурье

Cm — комплексная амплитуда Фурье-компоненты (или гармоники).

Отдельные слагаемые в сумме называют: фурье-компонента, гармоника, мода.

Фурье гармоники экспоненты являются собственными функциями оператора ˆ для ли

- ( ) L ( -

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

Таким образом, для установления устойчивости достаточно определить собственное зна-

чение оператора ˆ соответствующее произвольной Фурье гармонике и убедиться что

G L, - , ,

|G| ≤ 1.

6

Слайд №7: Метод Неймана: схема ВВЦП.

Динамика погрешности определяется оператором

Представим возмущение на n-ом временном слое в виде гармоники с произвольным волновым числом km:

εn(x) = exp(jkmx) .

Тогда решение на временном (n + 1)-ом слое будет иметь вид:

εn+1(x) = G exp (jkmx) ,

где G — в общем случае комплексное число, называемое множителем перехода (коэффициентом перехода, коэффициентом усиления). Условием устойчивости РС является соотношение

|G| ≤ 1 .

Имеем:

εni = exp(jkmxi) , εni+1 = exp(jkmxi+1) , εni−1 = exp(jkmxi−1) , εni +1 = G exp(jkmxi) .

Подставляя в первое выражение, получаем соотношение, определяющее коэффициент перехода,

G exp(jkmxi) = λ exp(jkmxi+1) + (1 − 2λ) exp(jkmxi) + λ exp(jkmxi−1) .

Учитывая, что xi+1 = xi + h, xi−1 = xi − h, после сокращения на exp (jkmxi) приходим к

G = 1 − 4λ sin2 km2h .

Теперь условие устойчивости дает систему неравенств

1 − 4λ sin2 (kmh/2) ≤ 1 ,1 − 4λ sin2 (kmh/2) ≥ −1 .

Первое неравенство выполняется автоматически, т.к. λ — положительно. Второе неравен-

Т.к. это должно выполняться для любых волновых чисел, то окончательно получаем усло-

вие устойчивости РС ВВЦП в виде

λ ≤ 12

или

h2 τ ≤ τc = 2a .

При любом τ > τc проявляется динамическая неустойчивость схемы ВВЦП.

7

Слайд №8: Метод Неймана: неявная схема для уравнения теплопроводности.

Динамика погрешности определяется оператором

Представим возмущение на n-ом временном слое в виде гармоники с произвольным волновым числом km:

εn(x) = exp(jkmx) .

Тогда решение на временном (n + 1)-ом слое будет иметь вид:

εn+1(x) = G exp (jkmx) ,

где G — в общем случае комплексное число, называемое множителем перехода (коэффициентом перехода, коэффициентом усиления). Условием устойчивости РС является соотношение

|G| ≤ 1 .

Имеем:

εni = exp(jkmxi) , εni +1 = G exp(jkmxi) , εni−+11 = G exp(jkmxi−1) , εni+1+1 = G exp(jkmxi+1) .

Подставляя в первое выражение, получаем соотношение, определяющее коэффициент перехода,

G exp(jkmxi) = λG [exp(jkmxi+1) + −2 exp(jkmxi) + exp(jkmxi−1)] + exp (jkmxi) .

Учитывая, что xi+1 = xi + h, xi−1 = xi − h, после сокращения на exp (jkmxi) приходим к

1

G = 1 + 2λ (1 − cos(kmh)) .

Видим, что условие устойчивости выполняется при любых значениях km и λ. Неявная схема для уравнения теплопроводности абсолютно устойчива.

•Неявные РС, как правило, абсолютно устойчивы. В этом их основное преимущество перед явными РС.

8

Слайд №9: Устойчивость РС для волнового уранения первого порядка.

∂∂tU + v · ∂∂xU = 0

описывает перенос (распространение вдоль оси x) со скоростью v некоторой величины U, характеризующей сплошную среду. Проанализируем для этого уравнения по методу

Неймана устойчивость некоторых РС.

|G| = q12 + C2 sin2(kmh) , где C = vhτ

— число Куранта. Поскольку при любых C и km коэффициент перехода |G| > 1, схема ВВЦП для одномерного волнового уравнения абсолютно неустойчива.

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

Для коэффициента перехода имеем

q

G = (1 − C + C cos(kmxi))2 + C2 sin2(kmxi) .

Анализ показывает, что условие |G| ≤ 1 выполняется, если 0 ≤ C ≤ 1, т.е. РС устойчива

при выполнении

0 ≤ vhτ ≤ 1 .

Это — условие Куранта – Фридрихса – Леви (КФЛ).

Таким образом, данную РС можно применять только в случае, если v > 0, т.е. когда волна движется в положительном направлении оси Ox.

Схема вперед по времени и с правой разностью по координате

условно устойчива для v < 0 при выполнении условия:

|vh|τ ≤ 1 .

9

Слайд №10: Критерии устойчивости РС для волновых уранений.

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

Во-первых:

Как можно интерпретировать величину h/v? h/v = tp, где tp — время распространения физического свойства U на расстояние, равное одному пространственному шагу сетки. Условие КФЛ требует, чтобы vτ/h ≤ 1, т.е. τ ≤ tp. Это означает, что выбираемый шаг по времени должен быть меньше, чем время распространения физического процесса на один пространственный шаг сетки.

Во-вторых:

Изобразим схематически направление потоков физической величины U в случае v > 0 и в случае v < 0. На этой же схеме изобразим узлы пространственной сетки. Если v > 0, то

устойчивой является схема, в которой пространственная производная аппроксимируется левой разностью, т.е. для вычисления используются значения величины U в узлах i и i−1. Если v < 0, то устойчивой является схема, в которой пространственная производная аппроксимируется правой разностью, т.е. для вычисления используются значения величины U в узлах i и i + 1. Таким образом, устойчивой всегда оказывается РС, в которой для вычислений используется информация из области ”вверх по потоку”.

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

10