Научная литература
booksshare.net -> Добавить материал -> Физика -> Федоренко Р.П. -> "Введение в вычислительную физику" -> 81

Введение в вычислительную физику - Федоренко Р.П.

Федоренко Р.П. Введение в вычислительную физику — М.: Физ-тех, 1994. — 528 c.
ISBN 5-7417-0002-0
Скачать (прямая ссылка): vvedenievvichesleniyah1994.djvu
Предыдущая << 1 .. 75 76 77 78 79 80 < 81 > 82 83 84 85 86 87 .. 210 >> Следующая

ЖЕСТКИЕ СИСТЕМЫ ОДУ

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)-окрестности Г поле направлений почти горизонтально, фазовая скорость очень велика (порядка
Предыдущая << 1 .. 75 76 77 78 79 80 < 81 > 82 83 84 85 86 87 .. 210 >> Следующая

Реклама

c1c0fc952cf0704ad12d6af2ad3bf47e03017fed

Есть, чем поделиться? Отправьте
материал
нам
Авторские права © 2009 BooksShare.
Все права защищены.
Rambler's Top100

c1c0fc952cf0704ad12d6af2ad3bf47e03017fed