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

a1	0	0	…	0	0
a2	b21	0	…	0	0
a3	b31	b32	…	0	0
…	…	…	…	…	…
ak	bK1	bK2	…	bKK–1	0
	c1	c2	…	cK–1	cK

В этом случае для расчёта vn+1 по vn в соответствии с (21) имеем простые рекуррентные соотношения:

— предиктор,

(22)

— первый корректор и т. д.

«Полуявные» или «диагонально неявные» схемы Рунге–Кутты отличаются от явных наличием ненулевых элементов на главной диагонали таблицы Бутчера (L≤k):

a1	b11	0	…	0
a2	b21	b22	…	0
…	…	…	…	…
ak	bK1	bK2	…	bKK
	c1	c2	…	cK

В этом случае на каждом шаге по времени последовательно решаются нелинейные системы:

,

и т. д.; например, методом Ньютона:

  

с неким начальным значением , например, .

В общем случае ( L = K ), обозначая r = {r1, …, rk}, имеем еще более громоздкую нелинейную систему

R(r) = 0, где R = {R1, R2, …, RK}

и итерационный процесс для ее решения

.

1.3.2. Аппроксимация

Параметры схемы (неопределенные пока коэффициенты a1, …, aK; c1, …, cK; b11, …, b1K, …, bK1, …, bKK), конечно, не произвольны. Их (все или часть) находят, прежде всего, из условий аппроксимации, получаемых из разложения (21) в ряд в точке t = tn или t = tn+1 с учетом (18) и его следствий

vt = f(t, v),

(т. е. рассматривая аппроксимацию на решениях (18)).

Поскольку в (21) при rk стоит множитель τ, то в разложении правой части (21) нужно удерживать на один член ряда меньше, чем в левой. Индекс n будем опускать. Итак, левая часть (21):

Далее:

Если коэффициенты ak выбирать таким образом, чтобы

, (23)

то Подставляя эти разложения в (21) и группируя члены при одинаковых степенях τ, получим (после деления на τ )

Очевидно, что при выполнении условия (вместе с (23))

(24)

обеспечивается 1-й порядок аппроксимации, при

(25)

(вместе с (23), (24)) — 2-й порядок аппроксимации. Условия 3-го порядка точности (вместе с (23) – (25)) есть

, (26)

а условия 4-го порядка (вместе с (23) – (26)):

,

(27)

.

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

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

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

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

(28)

Отсюда

, (29)

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

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

3(b11 + b22) – 1 – 6(b11b22 – b12b21) = 0, (30)

3[(b11 + b11) + (b21 + b21)] – 6(b22 + b21)(b11 + b12) – 2 = 0.

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

b22 = (1 – 3b11)/[(1 – 2b11)], b11  1/2, (31)

b21 = 1/[3(1 – 2b11)],

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

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

vn+1 = q vn, (32)

.

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

v1 = vn + s(b11v1 + b12v2 +…+ b1KvK),

v2 = vn + s(b21v1 + b22v2 +…+ b2KvK),

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

vK = vn + s(bK1v1 + bK2v2 +…+ bKKvK).

Обозначая v = {v1, …, vK},  vn = {vn, …, vn},  B = , имеем (E – sB)v = vn, откуда v = (E – sB)–1vn = Dvn,  D = {dkj} (т. е. vk = akvn =). Подставляя это в (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 (точное решение — сплошная кривая).

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

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

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

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

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

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

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

а	б
Рис. 1.10. Области устойчивости схем Эйлера
а	б
Рис. 1.11. Области устойчивости жёстко-устойчивых схем

Для абсолютно устойчивых схем область устойчивости — вся комплексная плоскость 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б), оставаясь внутри области устойчивости.

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

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

|q(s)| ® 0 при |s| ®¥ (34)

(или Re q(s) ® 0 при Re s ® –¥).

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

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

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

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

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

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

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

b11b12 – (b11 + b12) + 1/2 = 0. (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) называется схемой Лобатто, и для неё таблица Бутчера имеет вид

0	0	0		0	0	0
½	0	½	или	b21+½	b21	½	(36)
	0	1

(более общий вариант, b21 ≠ –1/2).

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

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

1	1	0		1	1	0
0	–½	½	или	b21+½	b21	½	(37)
	½	½					(b21 ≠ 1/2).

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

1	1	0		1	1	0
½	–½	1	или	b21+½	b21–½	1	(38)
	0	1					(b21 ≠ 1/2).

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