где

r = f{t + Шг, y + MZ ЬгэЛ Г = 1, 2. ..., /?, (7.51Ь)

ar=tbrs. г = 1, 2. .... /?. (7.51с)

Формулы (7.51) охватывают неявные схемы Рунге-Кутты, неэкономичные в вычислительном плане, так как уравнение (7.51 Ь) является нелинейным и должно решаться с применением итераций для каждого из значений f* и на каждом шаге ло времени. Поэтому больший интерес проявляется к явным <:хемам Рунге -Кутты, в которых верхние пределы суммирования в формулах (7.51Ь) и (7.51с) заменяются на г-1. Так, например, схему Эйлера (7.48) можно интерпретировать в качестве одноэтапной (=1) явной схемы Рунге - Кутты при А=1, Ьи = 0, ai=0.

Существует бесконечное число двухэтапных явных схем Рунге - Кутты второго порядка. Типичной схемой такого рода является улучшенная схема Эйлера

u = u + Atr> ,7 52)

и+1=г/ + 0.5Л/(Г + Г).

Существуют и многочисленные трехэтапные явные схемы Рунге - Кутты третьего порядка [Lambert, 1973]. Нижеследующая четырехэтапная и четырехшаговая схема Рунге - Кутты четвертого порядка получила широкое распространение при решении задач с начальными значениями, определяемых обыкновенными дифференциальными уравнениями

п-ы и-+{Г + 2Г + 2Г + Г), (7.53)

и- = + 0.5Д/Г, r = f (t-\u% = и + 0.5 А/Г, Г = f {t \ tt), = ц + д -, r = f{t\ О.

Очевидно, что схемы Рунге - Кутты повышенного поряда требуют большего числа оценок для f. При решении таких полудискретных систем, как (7.43), серия повторных оценок всех членов, связанных с пространственной дискретизацией, зачастую играет определяющую роль в подсчете суммарного времени исполнения.

Если скалярное уравнение (7.44) рассматривается вместе

с начальным условием u{0)=uo, то соответствующее точное решение имеет форму

a{t) = Uoe = Uss. {7.54}

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

Е(а/-ЯАР/) = 0. (7.55)

При заданном XAt абсолютная устойчивость достигается [Lambert, 1973] тогда, когда все корни Гз полинома (7.55) удовлетворяют условию

г,<1. (7.56)

Это условие в принципе эквивалентно ограничению, выполнение которого требовалось для устойчивости решения дискретизированного дифференциального уравнения в частных производных при применении анализа устойчивости по Нейману (§ 4.3).

В случае схемы Эйлера абсолютная устойчивость достигается, если -2 <! XAt 0. Поскольку Л < О, это соответствует условию Д<2/2. Для трапецоидальной схемы абсолютная устойчивость связана с условием ЯА < 0; следовательно, трапецоидальная схема устойчива при любом выборе А.

Критерий абсолютной устойчивости, эквивалентный (7.56), описан в книге [Lambert, 1973] и приводит к ограничению -2 XAt О, относящемуся как к одноэтапной схеме Эйлера, так и к явной двухэтапной схеме Рунге-Кутты (7.52). Однако четырехэтапная явная схема Рунге-Кутты (7.53) имеет слегка расширенный интервал устойчивости -2.78 XAt 0.

При применении метода прямых мы встречаемся с системами обыкновенных дифференциальных уравнений. Следовательно, уместно проверить абсолютную устойчивость уравнений в форме (7.42). Как и прежде, можно получить полином устойчивости (7.55), если вместо X вводить - собственные значения матрицы А. В общем случае величины Xk комплексные. Однако чтобы задача имела физически устойчивое решение, необходимо, чтобы вещественная часть Xhy т. е. ReXh, была положительной.

Например, применение схемы Эйлера (7.48) к уравнению (7.41) дает следующий критерий абсолютной устойчивости:

-2<НеЯл. Д/<0. (7.57)

Для системы из М уравнений, соответствующей представлению (7.41), матрица А получается трехдиагональной и ее собственные значения определяются формулами (9.48), т. е.

kf, = а (-2+ 2 cos m/jM + D) 7

В результате наиболее строгое ограничение на величину А, вытекающее из условия (7.56), имеет форму

-2< или А/<-. (7.59)

Таким образом, здесь получено то же самое ограничение на величину Д, обусловленное устойчивостью, как и для схемы ВВЦП (7.5). Этого и следовало ожидать, так как метод Эйлера с шагами по времени в применении к уравнениям (7.41) совпадает со схемой ВВЦП. Применение к (7.41) трапецоидаль-ной схемы (7.49) приводит к ограничению ReXk-AtO. Но так как в силу (7.58) ReXk О, то ограничение на отсутствует. Трапецоидальная схема, примененная к (7.41), совпадает со схемой Кранка - Николсона (7.22). Как можно предположить на основании этих примеров, абсолютная устойчивость алгоритма решения системы обыкновенных дифференциальных уравнений эквивалентна устойчивости алгоритма решения родственной задачи о решении дискретизированного дифференциального уравнения в частных производных (§ 4.3).

Если система, полученная при помощи метода прямых, является нелинейной, например система (7.43), то определение устойчивости алгоритма зависит от собственных значений якобиана dF/dV. Однако сами собственные значения зависят от решения и могут нуждаться в корректировке в процессе построения решения.

Точное решение системы уравнений (7.42), образуемой из {7.41), может быть записано в виде

U=S a.A4 + U.s, (7.60)

где собственные значения Xk определяются формулами (7.58), а величины tk - соответствующие собственные векторы. Коэффициенты ak определяются начальными условиями Uo(/). Как ясно из формул (7.58), диапазон изменения собственных значений