Введение в вычислительную физику - Федоренко Р.П.
ISBN 5-7417-0002-0
Скачать (прямая ссылка):
Мягкий спектр. Собственные значения и векторы обозначим ХДх), и Фу(х) (J = 1, 2, ..., У). Для мягкого спектра выполняются условия
I Я,-(X) H/«L. (4)
Время интегрирования T является средним относительно I и очень большим относительно L: IT равно, например, 10, 20, 30, а LT может быть порядка IO3, IO6 и больше. Отношение LU называют показателем жесткости системы. В приложениях встречаются ситуации, когда LU равно IO6, IO9, IO15. Будем считать, что I= 1.
Стандартные методы интегрирования требуют шага т*, малого в том смысле, что TtIIZjcII«. 1. Так как HZjcII — IAiI « L, то t*<s< 1/L и расчет требует LT шагов численного интегрирования, что иногда просто неприемлемо. С точки зрения гладкости квазистационарной части траектории приемлем большой шаг, например т « 10_3Т. Основной проблемой в теории жестких систем является разработка ал-
§17] ЖЕСТКИЕ СИСТЕМЫ ОДУ 211
горитмов численного интегрирования с таким большим шагом. Она была решена на основе неявных схем. Ниже мы объясним, почему оказалось достаточно такого относительно несложного вычислительного аппарата.
Заметим еще, что нелинейная система может быть жесткой в одной части фазового пространства и не быть таковой — в другой. Числа I, J следовало бы обозначать 1(х), J(x) (/ + / = N).
Линейные жесткие системы. Этот удобный объект позволяет оценить возможности неявных схем, сравнивая точное решение с приближенным. Рассмотрим систему х = Ах, х(0) = X0, считая постоянную матрицу А жесткой (если f(x) =Ax, то Jx = А). Точное решение дается явной формулой
*(0=2 CieA‘t ф/ + 2 cJe^it Vr ‘ J
Первое слагаемое убывает, как е~и, и становится пренебрежимо
малым вне пограничного слоя — интервала времени [о, O(^p)];
второе слагаемое представляет «квазистационарное» движение x(t) (см. рис. 21).
Попытаемся интегрировать систему по явной схеме Эйлера (см. § 5):
(*п+1-XaVx =Axa, или хп+1 = (Е+ хА)хп. (6)
Решение (6) легко получить в терминах спектра:
= 2 cA1 + хАі)п ф1 + 2 cA1 + x^n Vj- (7)
Для мягкой компоненты при шаге 1 (на практике это означает, например, хI « 0.1) имеем 1 + = exXt(\ + 0(т2/2)) и
(1 + xkj)n = enxXi (I + 0(т/)) при п як Т/х.
Таким образом, мягкая компонента (7) аппроксимирует соответствующую компоненту точного решения (5) в обычном смысле слова. Для жесткой компоненты (1 + тА)" « (—xL)n, и при tL»1 (а на практике это величины порядка IO3 IO6 и т.д.) за несколько шагов числа хп просто выходят из разрядной сетки ЭВМ. Рассмотрим неявную схему:
(хп+1~ хп)/х = Ахп+1, или xn+l = (E-XA)-lXn. (8)
Точное решение (8) имеет вид
= 2 c^1 - хАіУП ф1 + 2 cA1 - ХЬУП Vj- (9)
212
ЧИСЛЕННЫЕ МЕТОДЫ МАТЕМАТИЧЕСКОЙ ФИЗИКИ
[Ч. и
Мягкая компонента так же аппроксимирует соответствующую часть (5), а жесткая быстро стремится к нулю (как (1/тL)") и, таким образом, качественно правильно описывает поведение «погранслойно-го» слагаемого (5). Обычно в прикладных задачах подобного рода не представляют особого интереса ни структура пограничного слоя, ни его длительность (лишь бы она была много меньше Т).
Наиболее интересным содержательным результатом является квазистационарный режим. Если это так, решение (9) нас устраивает. На рис. 21 крестиками показано удовлетворительное приближенное решение: толщина пограничного слоя (т или 2т, Зт, ...) существенно больше его реальной толщины Структура слоя
совершенно не описывается, но квазистационарный режим описан достаточно аккуратно. Кстати, если по каким-то причинам нас интересует и пограничный слой, его расчет малым шагом t'«1/L не представляет труда и может быть выполнен по любой явной схеме.
При численном интегрировании линейной системы х = A(t)x с Медленно меняющейся матрицей (в том смысле, ЧТО ||Л||т<3< ||Л||) иногда используют квазианалитические методы численного интегрирования:
*„+1 = еАпгхп. (10)
Для их эффективной реализации нужно вычислять матричную экспоненту. В жесткой системе (при ||Л||т55>1) это не просто. Использовать ряд Тейлора практически нельзя как по «экономическим» причинам (предоставим читателю оценить необходимое число членов ряда, оно порядка TL), так и по причине вычислительной неустойчивости. В силу (3), (4) имеем еАх = 0(1), но эта величина получается суммированием тех же величин, что и е~Ах = 0(eLx).
Метод удвоения аргумента. Для вычисления матричной экспоненты разработан достаточно эффективный алгоритм. Выберем некоторое целое р, такое, что \\А\\х/2р^1. Тогда еАтПР ^B = = E 4- (т/2р)А, а еАх « в2\ последняя же матрица может быть вычислена за р умножений матриц (вычисляются значения Bl = BxB, B2= B1X B1 и т.д.). Очевидно, что р= 0(ln (xL)), и с точки зрения числа операций этот метод вполне эффективен. Однако его применение требует известной осторожности. В самом деле, мягкое собственное значение В есть I -I- (x/2p)Xj + є, где є — относительная погрешность представления чисел в ЭВМ. Пусть, например, xL/2p = 1. Тогда в В информация о Xj передается с относительной погрешностью, не меньшей I txL/Xj I. При некоторых вполне реальных соотношениях между жесткостью и длиной машинного слова величина |ЗуР= (1 + (т/2р)Ay)2” может не иметь ничего общего с
