Введение в вычислительную физику - Федоренко Р.П.
ISBN 5-7417-0002-0
Скачать (прямая ссылка):
ЖЕСТКИЕ СИСТЕМЫ ОДУ
213
ехіг и квазистационарный режим будет рассчитан неверно. Плохо и то, что погрешность такого рода не имеет явных признаков: численное решение будет гладким и будет выглядеть внешне таким, каким оно могло бы быть. Поэтому численное интегрирование жестких систем на машинах типа IBM (четырехбайтное слово) ведется как минимум с двойной точностью.
Системные методы численного интегрирования. Вычисление матричной экспоненты используется в некоторых алгоритмах. В качестве характерного примера опишем метод, разработанный в Институте химической физики. Стандартный шаг численного интегрирования состоит из следующих операций.
Пусть точка хп уже известна. Вычисляется матрица А = fx(xn). Дальнейшее основано на некоторых преобразованиях. Сделаем замену переменных x(t) — хп + u(t). Для и имеем систему
м — Au(t) = f(xn + и) — Au(t) (11)
с начальными данными u(tn) = 0. Положим, для простоты, tn = 0 и проинтегрируем (11) на интервале [0, т], используя оператор (d/dt — A)'1. Тогда.
T
ы(т) = еАх J e~At [f(xn + u(t)) — Au(t)] dt.
о
Это точная формула.
Приближенную формулу можно получить, взяв в правой части и(х) вместо u(t). Основания для этого следующие. Множитель e~At есть величина, быстро возрастающая вправо. Поэтому основной вклад в интеграл дают значения на правом конце интервала интегрирования. Кроме того, по смыслу замены u(t) есть малая величина:
f(xn + и) Kf(Xn) + fx(xn)u + 0(||м||2),
т.е. выражение в квадратных скобках почти (с точностью до 0(11 mII2)) постоянно. С учетом указанной аппроксимации это выражение выносится из интеграла, после чего он явно вычисляется:
г
и(т) « еАг J e~At dt If (хп + м(т)) — Аи(х)\ =
о
= A'l(eAx - Е) [/(*„ + и(х)) — Ли(т)].
Для определения функции ы(т) получено нелинейное уравнение, которое решается, как показал опыт, простыми итерациями типа ы(1+1) = В [f(xn + u(l)) — AuW]. Ho сначала нужно вычислить
214
ЧИСЛЕННЫЕ МЕТОДЫ МАТЕМАТИЧЕСКОЙ ФИЗИКИ
[Ч. II
матрицу В =* А~1(еЛх — Е). Дело осложняется тем, что часто А~1 не существует. Это естественное свойство тех систем, в которых компоненты X суть относительные концентрации некоторых веществ. Тогда система (1) имеет очевидный первый интеграл — закон сохранения (х(0, е) = 1, где е= {1, 1, ..., 1}. Дифференцируя по t, получаем (х, е) = (/(*), е) =0. Дифференцируя по х, находим f Х(х)е = 0, т.е. нуль есть точка спектра матрицы f Х(х) при всех х, а е — соответствующий собственный вектор.
Определение В корректно (так как (еАх — Е)е = 0), и нужно правильно раскрыть неопределенность <»•(). Для этого используется специфическая форма алгоритма удвоения аргумента. Она основана на следующих преобразованиях. Обозначим А~1(еЛх — Е) = Bx. Тогда
Bi = А~1(еЛх -E) = А~1(еАг12 - Е)(еАг12 + E) =
= А~\еАгП - Е)АА~1(еАх12 -Е+2Е)= Вх/2(АВх12 + 2 Е).
Эта формула позволяет эффективно вычислять Bx.
Вернемся к алгоритму численного интегрирования. Вычисляется матрица
Вх/2р = A~l(eA%nP - Е) » А~1(Е + (х/2р)А -E) = (х/2р)Е.
Далее за р шагов удвоения аргумента находится Bx. Итерациями вида, например,
M('+1) = Bx[f(xn + и(,)) — Au^]
находится с заданной погрешностью ы(т). Скорость сходимости регулируется выбором шага т. Если потребовалось слишком много итераций, следующий шаг интегрирования выполняется с меньшим шагом т. Если итерации сходятся слишком быстро, шаг т увеличивается.
Перейдем к более сложным задачам.
Сингулярно-возмущенные системы. Рассмотрим класс систем, также относящихся к числу «жестких», для которых уже давно создана асимптотическая теория и качественная картина поведения траекторий достаточно ясна. Это позволяет четко формулировать требования к методам приближенного интегрирования и проверять, в какой мере те или иные методы интегрирования таким требованиям удовлетворяют.
Итак, рассмотрим систему
к = Lf (х, у) (или гх = /, е = L-i<k1),
t \ (12) У=<р(х, у).
а
§ 17]
ЖЕСТКИЕ СИСТЕМЫ ОДУ
215
Здесь /, ф и их производные — величины порядка 0(1); х, у — векторы размерности /, / соответственно. Спектр вариационной матрицы (12) определяется уравнением
det
(Lfx-IE1 Vx
Vy
Vy -
= 0.
(13)
Легко показать, что жесткая часть спектра определяется спектром матрицы }х (умноженным на L, если, конечно, f х не имеет собственных значений O(LTi)), а соответствующие собственные векторы Ф(- имеют по существу лишь х-компоненту (их у-компонента
есть 0(LT1)). Хорошо известна и качественная структура траекторий. Она определяется многообразием Г, уравнением которого является f(x, у) = 0. Оно разбивает фазовое пространство на две части: f(x, у) > 0 и f(x, у) <0. Система (12) является жесткой в случае «отрицательности» спектра /Х(х, у).
Теория систем вида (12) хорошо развита. На рис. 23 показана типичная картина в случае, когда х и у — скаляры. (В многомерном случае картина, конечно, более сложная, но примерно того же характера.) Вне малой O(LTi)-окрестности Г поле направлений почти горизонтально, фазовая скорость очень велика (порядка