Устойчивость

Конечно из того, что (21) при стремится к (18) еще не следует, что их решения будут сближаться. Как известно, для сходимости, по крайней мере, в линейном случае (), необходима еще устойчивость разностной схемы. А поскольку в линейном случае система (18) переходит в совокупность простейших уравнений типа (1) (толькоможет быть и комплексным, т. е.), при анализе устойчивости можно ограничиться изучением (1). Чтобы избежать громоздких выкладок и иметь возможность графического представления некоторых результатов, ограничимся также случаемK= 2 в (21).

Итак, в случае модельного уравнения для обеспечения второго порядка аппроксимации из (23) – (25) имеем:

(28)

Отсюда

, (29)

т. е. четыре коэффициента произвольны, а остальные четыре коэффициента определяются соотношениями (28), (29).

Из (26) с учетом (23), (29) имеем для схем 3-го порядка

3(b₁₁ + b₂₂) – 1 – 6(b₁₁b₂₂ – b₁₂b₂₁) = 0, (30)

3[(b₁₁ + b₁₁) + (b₂₁ + b₂₁)] – 6(b₂₂ + b₂₁)(b₁₁ + b₁₂) – 2 = 0.

Для полуявных схем (b₁₂= 0) (30) определяет однопараметрическое семейство схем 3-го порядка с коэффициентами

b₂₂ = (1 – 3b₁₁)/[(1 – 2b₁₁)], b₁₁  1/2, (31)

b₂₁ = 1/[3(1 – 2b₁₁)],

а на 4-й порядок уже не хватает свободных коэффициентов.

Полагая ,(σ — число Куранта), получим:

vⁿ⁺¹ = q vⁿ, (32)

.

Действительно, для модельного уравнения в (21) r_k = v_k,

v₁ = vⁿ + (b₁₁v₁ + b₁₂v₂ +…+ b_1Kv_K),

v₂ = vⁿ + (b₂₁v₁ + b₂₂v₂ +…+ b_2Kv_K),

..........................................

v_K = vⁿ + (b_K1v₁ + b_K2v₂ +…+ b_KKv_K).

Обозначая v = {v₁, …, v_K},  vⁿ = {vⁿ, …, vⁿ},  B = , имеем (E – B)v = vⁿ, откуда v = (E – B)^–1vⁿ = Dvⁿ,  D = {d_kj} (т. е. v_k = _kvⁿ =). Подставляя это в (21), и получим (32).

Геометрическая прогрессия (32) будет совпадать с точным решениемvⁿ⁺¹=vⁿe^, еслиq() = e^. Этот способ выбора коэффициентов в (21) называютметодом экспоненциальной подгонки, его обобщение на случай линейной системы (7) ограничен лишь возможностями отыскания спектра матрицыA, а в случае нелинейной системы (18) может быть использована линеаризация (18) на малом отрезке [tⁿ, tⁿ⁺¹].

Так как мы рассматриваем здесь только случай Re  < 0 (затухающие или устойчивые по Ляпунову решения), то при |q()| > 1 приближенное решение (32) не будет иметь ничего общего с точным решением. При действительных значенияхчисленное решение ведет себя, как показано штриховыми линиями на рис. 1.9 (точное решение — сплошная кривая).

Рис. 1.9 Вид численного решения при различных q

|q()| ≤ 1. (33)

Для действительных это эквивалентно–1 < q() < 1.

Определение 1.2. Схема называется абсолютно устойчивой, если (33) выполняется при всех значениях 

Определение 1.3. Схема называется монотонной, если для действительных  выполняется 0 < q() < 1.

Если в комплексной плоскости (Re , Im ) нарисовать кривую |q()| = 1, то она будет ограничивать область устойчивости. Примеры областей устойчивости для явного и неявного методов Эйлера показаны (заштрихованы) на рис. 1.10а и рис. 1.10б соответственно.

Для абсолютно устойчивых схем область устойчивости — вся комплексная плоскость . Таких схем на практике не существует, но для жестких уравнений этого и не нужно, достаточно, чтобы схема былаA-устойчивой (Далквист, 1963 г.), как, например, неявная схема Эйлера, рис. 1.10б.

Определение 1.4. Схема называется A-устойчивой, если кривая |q()| = 1 лежит в правой полуплоскости σ.

Это требование довольно тяжелое, поэтому его еще более ослабляют, требуя, чтобы кривая |q()| = 1 лежалавне заштрихованных на рис. 1.11а, б областей (Гир, 1969г.). Это есть разные определенияжестко-устойчивыхсхем. Действительно, для жесткой системы ОДУ (рис. 1.4г), если все Re _j < 0, то при измененииот 0 довсе точки_j_jна плоскости {Re, Im} будут двигаться по радиальным направлениям внутри некоторого угла 2β (рис. 1.11б), оставаясь внутри области устойчивости.

Величина Re q() определяет скорость затухания решения. При Re << –1 точное решение затухает очень быстро, значит и численный метод при ||(точнее, при Re–) должен обладать этим свойством.

Определение 1.5. Схема называется L-устойчивой (Ламберт, 1973 г.), если

|q()|  0 при ||  (34)

(илиRe q()  0приRe  –).

Обобщение всех этих определений на случай линейной системы (7) очевидно, надо, чтобы соответствующие условия выполнялись для всех _j= _j(_j— собственные значения матрицыA). Для общей нелинейной системы (18) все это должно выполнятся локально, при всехtⁿ(принцип линеаризации и замораживания коэффициентов).