Разностные схемы для решения жестких обыкновенных дифференциальных уравнений в пространстве неопределенных коэффициентов (стр. 5 )

vn+1 = q vn, (32)

.

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

v1=vn + s (b11 v1 + b12 v2 +...+ b1k vk),

v2=vn + s (b21 v1 + b22 v2 +...+ b2k vk),

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

vk=vn + s (bk1 v1 + bk2 v2 +...+ bkk vk),

Или, обозначая v={v1,...,vk }, vn={v1,...,vn }, B= , имеем: (E – sB) v=vn, откуда v=(E – sB)vn=Dvn, D={d kj} (т. е. vk=ak vn=). Подставляя это выражение в (21), и получим (32).

Геометрическая прогрессия (32) будет совпадать с точным решением vn+1 = vn es, если q(s)=es. Этот способ выбора коэффициентов в (21) называют методом экспоненциальной подгонки, его обобщение на случай линейной системы (7) ограничен лишь возможностями отыскания спектра матрицы A, а в случае нелинейной системы (18) может быть использована линеаризация (18) на малом отрезке [tn,tn+1].

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

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

|q(s)|<1 (33)

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

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

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

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

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

А	Б
Рис. 1.10. Области устойчивости схем Эйлера

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

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

Это требование довольно тяжелое, поэтому его еще более ослабляют, требуя, чтобы кривая |q(s)|=1 лежала вне заштрихованных на рис. 1.11А, Б областей (Гир, 1969г.). Это есть разные определения жестко устойчивых схем. Действительно, для жесткой системы ОДУ (рис. 1.4Г), если все Re(lj)<0, то при изменении t от 0 до ¥ все точки sj = tlj на плоскости {Re s, Im s} будут двигаться по радиальным направлениям внутри некоторого угла 2β (рис. 1.11Б), оставаясь внутри области устойчивости.

А	Б
Рис. 1.11. Области устойчивости жёстко устойчивых схем

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

Определение 1.5. Схема называется L-устойчивой (Ламберт, 1973 г.), если
|q(s)|®0 при |s|®¥ (или Re q(s)®0 при Re s® –¥) (34)

Обобщение всех этих определений на случай линейной системы (7) очевидно, надо, чтобы соответствующие условия выполнялись для всех sj=t lj (l – собственные значения матрицы A). Для общей нелинейной системы (18) все это должно выполнятся локально, при всех tn (принцип линеаризации и замораживания коэффициентов).

2.3.4. Примеры схем Рунге–Кутты

Рис. 1.12. Полуявные схемы Рунге–Кутты в плоскости неопределённых коэффициентов (K=2)

Вернемся теперь к конкретному выражению для q из (32) и для простоты ограничимся случаем действительных l<0. Это особенно не меняет существа дела, но упрощает выкладки. Будем также полагать b12=0 (полуявные схемы), тогда q не зависит также и от b21, а оставшиеся два свободных коэффициента примем за оси координатной плоскости, в которой будем вести рассмотрение (рис. 1.12).

Область устойчивых при всех –¥<s< 0 схем, т. е. область A-устойчивых схем –1<q(s)<1 на рис. 1.12 заштрихована (горизонтальная штриховка: –1<q(–¥)<0, вертикальная штриховка 0<q(–¥)<1).

Множество схем третьего порядка аппроксимации (31) на рис. 1.12 показано сплошной кривой (гипербола с асимптотами b11=1/2, b22=1/2).

Для L-устойчивых схем из (32), (34) имеем

b11b12 – (b11+b12)+1/2=0, (b11+b12)≠1/2 (35)

На рис.1.12 L-устойчивым схемам (35) соответствует штриховая кривая (гипербола с асимптотами b11=1, b22=1). Видно, что в данном случае (среди полуявных схем Рунге–Кутты с K=2) ни L-устойчивых, ни монотонных схем третьего порядка нет, а L-устойчивые схемы второго порядка расположены на кривых A3A4΄A5΄ и A4.

Несколько конкретных примеров. На прямой b11=b22, проходящей через точки O, A4΄, B1, A2, A4, расположены схемы Бутчера. Точка A3 (и точка A3΄, так как точкам, симметричным относительно прямой b11=b22 соответствует формальная замена в (21), (22) индекса 1 на 2 и наоборот) соответствует схеме Лобатто, для которой таблица Бутчера имеет вид:

или (b21≠–1/2) (36)

Это полуявная (с явным предиктором, т. е. достаточно экономная по объему вычислений) схема второго порядка аппроксимации. Она расположена на границе области A-устойчивых схем, не L-устойчива (q(–¥)=–1) и монотонна лишь при –2<s<0.

Точка A1 (и симметричная ей точка A1΄) также расположена на границе области A-устойчивых схем, не L-устойчива (s(–¥)=–1), монотонна лишь при –1–<s< 0 и по своим свойствам мало отличается от схемы Лобатто. Соответствующая ей таблица Бутчера имеет вид:

или (b21≠1/2) (37)

Точка A2 является A-устойчивой, но не L-устойчивой схемой (s(–¥)=–1/2), монотонна при –1–<s<0, как и схема A1, имеет второй порядок аппроксимации, по своим свойствам несколько лучше предыдущих схем. Соответствующая ей таблица Бутчера имеет вид:

или (b21¹1/2) (38)

Точки A4 и A4΄ являются примерами A - и L-устойчивых схем (точка A4 также монотонна при всех отрицательных значениях s) второго порядка аппроксимации. В соответствующей им таблице Бутчера знак ‘+’ относится к схеме A4, а знак ‘–’ – к A4΄:

, b21 ¹ 0 (39)

Схема Хаммера–Холлингсвудта – второго порядка аппроксимации, полностью неявная, не L-устойчивая (q(–¥)=1, пример неудачной схемы):

Явная схема Эйлера – Коши (точка 0):