Fedorenko-RP-Vvedenie-v-vychislitelnuyu-fiziku (810773), страница 42
Текст из файла (страница 42)
Так как 11УД 1Л,~ Ь, то т*~с1/г. и расчет требует. г,Т шагов численного интегрирования, что иногда просто иепряемлемо. С точки зрения гладкости квазистационарной части траектории приемлем большой шаг, например т 10 зТ. Основной проблемой в теории жестких систем является разработка ал- з 121 211 жесткие системы оде горитмов численного интегрирования с таким большим шагом. Она была решена на основе неявных схем. Ниже мы обьясним, почему оказалось достаточно такого относительно несложного вычислительного аппарата.
Заметим еще, что нелинейная система может быть жесткой в одной части фазового пространства и не быть таковой — в другой. Числа /„Х следовало бы обозначать /(х), У(х) (! + У = /т). Линейные жесткие системы. Этот удобный обьект позволяет оценить возможности неявных схем, сравнивая точное решение с приближенным. Рассмотрим систему х = Ах, х(0) = Х„, считая постоянную матрицу А жесткой (если /(х) = Ах, то /„= А). Точное решение дается явной формулой х(!) = ~К' С,.е».1 Ф, + ~т' с е"/ ~р . (5) ! / Первое слагаемое убывает, как е с', и становится пренебре18имо ппв1 .
малым вне пограничного слоя — интервала времени (О, О! — )~; второе слагаемое представляет «квазистационарное» движение х(!) (см. рис. 21). Попытаемся интегрировать систему по явной схеме Эйлера (см. з 5): (х„»1 — х„)/с=Ах„, или х«м — — (Е+ тА)х„, (б) Решение (б) легко получить в терминах спектра: х„= ~ С1(1 + тЛ1)" Ф1 + ~ сг(1+ тХ!)" р .. (7) Для мягкой компоненты при шаге т! ««1 (иа практике это означает, например, т!ж0 1) имеем 1+ 12/ = е1~(1+ О(т~! )) и (1+ тА/)" = е"'"~ (1+ О(т!)) при л ж Т/т.
(8) (х„1 — х„)/т = Ах„, о илн х„= (Š— тА) 1х„. Точное решение (8) имеет вид х„=~ С,(1 — Л,.) " 01+~' с!(1 — т11!)-" р,. (9) Таким образом, мягкая компонента (7) аппраксимирует соответствующую компоненту точного решения (5) в обычном смысле слова. Для жесткой компоненты (1+ тА)" ( — тЕ)", и при т/.»~1 (а на практике зто величины порядка 102 †: 10» и т.д.) за несколько шагов числа х„просто выходят из разрядной сетки ЭВМ. Рассмотрим неявную схему: г1г числвнныв методы мхткмАтичзской физики 1ч. и Мягкая компонента так же аппроксимирует соответствующую часть (5) „а жесткая быстро стремится к нулю (как (1/т/,)") и, таким образом, качественно правильно описывает поведение «погранслойногоь слагаемого (5). Обычно в прикладных задачах подобного рода не представляют особого интереса ни структура пограничного слоя, ни его длительность (лишь бы она была много меньше Т).
Наиболее интересным содержательным результатом является квазистационарный режим. Если это так, решение (9) нас устраивает. На рис. 21 крестиками показано удовлетворительное приближенное решение: толщина пограничного слоя (т или 2т, Зт, ...) существенно больше его реальной толщины О( — /. Структура слоя гы г.'1 совершенно не описывается, но квазистацнонарный режим описан достаточно аккуратно. Кстати, если по каким-то причинам нас интересует и пограничный слой, его расчет малым шагом т*«м 1Ю не представляет труда и может быть выполнен по любой явной схеме. При численном интегрировании линейной системы х= А(1)х с медленно меняющейся матрнцей (в том смысле, что йАйтмс'йАй) иногда используют квазианалитические методы численного интегрирования: х„« ~ — — е".'х„. (10) Для их эффективной реализации нужно вычислять матричную экспоненту.
В жесткой системе (при йАйт))1) это не просто. Использовать ряд Тейлора практически нельзя как по «экономическим» причинам (предоставим читателю оценить необходимое число членов ряда, оно порядка тЬ), так и по причине вычислительной неустойчивости. В силу (3), (4) имеем е"' = 0(1), но эта величина получается суммированием тех же величин, что н е "' = 0(ес'). Метод удвоения аргумента.
Для вычисления матричной экспоненты разработан достаточно эффективный алгоритм. Выберем некоторое целое р, такое, что йАй т/2е~ 1. Тогда е"мы ж В = = Е+ (т/2")А, а е"' Вг, последняя же матрица может быть вычислена за р умножений матриц (вычисляются значения В, = В х В, В = В, х В, и т д). Очевидно, что р = 0(1п (тЬ)), и с точки зрения числа операций этот метод вполне эффективен. Однако его применение требует известной осторожности. В самом деле, мягкое собственное значение В есть 1+ (т/2е) Х, + «, где е — отно/ сительная погрешность представления чисел в ЭВМ.
Пусть, например, т/./2е = 1. Тогда в В информация о Х передается с относительной погрешностью, не меньшей ) ет//Х ~. Прн некоторых вполне реальных соотношениях между жесткостью и длиной машинного слова величина В/з'=(1+(т/2г)Л )г может не иметь ничего общего с жзсткив снстамы одх е"р и квазистациоиарный режим будет рассчитан неверно. Плохо и то, что погрешность такого репа не имеет явных признаков: численное решение будет гладким и будет вмглядеть внешне таким, каким оно могло бы быть. Поэтому численное интегрирование жестких систем иа машинах типа 1ВМ (четырехбайтное слово) ведется как минимум с двойной точностью.
Системные методы численного интегрирования. Вычисление матричной экспоненты используется в некоторых алгоритмах. В качестве характерного примера опишем метод, разработанный в Институте химической физики. Стандартный шаг численного интегрирования состоит из следующих операций.
Пусть точка х„уже известна, Вычисляется 'матрица А = ~,(х„). Дальнейшее основано на некоторых преобразованиях. Сделаем замену переменных х(8) х„+ и(1). Для и имеем систему й — Аи(1) =У(х + и) — Аи(Г) (11) с начальными данными и(1„) = О. Положим, для простоты, Ф„' = О и проинтегрируем (11) на интервале 10, т), используя оператор (Ы/ой — А) '. Тогда, и(т) = е"' ~ е з' (У(х„+ иЯ) — Аи(1) '1 сй. о Это точная формула. Приближенную формулу можно получить, взяв в правой части и(т) вместо и(г). Основания для этого следующие. Множитель е " есть величина, быстро возрастающая вправо.
Поэтому основной вклад в иатеграл дают значения на правом конце интервала интегрирования. Кроме того, по сммслу замены и(г) есть малая величина: У(х„+ и) ж ~(х„) + /„(х„) и + ОЦ~ и й ~), т.е. выражение в квадратных скобках почти (с точностью до О(Йии~)) постоянно. С учетом указанной аппроксимации это выражение выносится из интеграла, после чего он явно вычисляется: т и(т) е"' ~ е "' Н1 (Дх„+ и(т)) — Аи(т)1= о = А '(е~' — Е) (/(х„+ и(т)) — Аи(т)1.
Для определения функции и(т) получено нелинейное уравнение, которое решается, как показал опыт, простыми итерациями типа ирен=В (У(х„+ ир1) — Аир11. Но сначала нужно вычислить гш численные методы метемлтическоа физики матрицу В = А '(е"' — Е). Дело осложняется тем, что часто А ' не существует. Это естественное свойство тех систем, в которых компоненты х суть относительные концентрации некоторых веществ.
Тогда система (1) имеет очевидный первый интеграл — закон сохранения (х(1), е) = 1, где е = (1, 1, ..., 1). Дифференцируя по 1, получаем (х, е) =(У(х), е) =О. Дифференцируя по х, находим ~„(х)е = О, т.е. нуль есть точка спектра матрицы ~„(х) при всех х, а е — соответствующий собственный вектор. Определение В корректно (так как (е"' — Е) е = О), и нужно правильно раскрыть неопределенность .О.
Для этого используется специфическая форма алгоритма удвоения аргумента. Она основана на следующих преобразованиях. Обозначим А '(ее' — Е) = В,. Тогда В = А >(е»' — Е> = А >(е»™ — Е)(е"Нз+ Е) = = А '(е"нз — Е)АА >(е~*>з — Е+ 2Е) = В, (АВ + 2Е). Эта формула позволяет эффективно вычислять В,. Вернемся к алгоритму численного интегрирования. Вычисляется матрица Вту2» = А '(е"из — Е) ж А '(Е+ (т/2»)А — Е) = (т/2»)Е. Далее за р шагов удвоения аргумента находится В,. Итерациями вида, например, ип+'>= В,Ц(х„+ир>) — Аиц>) находится с заданной погрешностью и(т). Скорость сходимости регулируется выбором шага т. Если потребовалось слишком много итераций, следующий шаг интегрирования выполняется с меньшим шагом т. Если итерации сходятся слишком быстро, шаг т увеличивается. Перейдем к более сложным задачам. Сннгулярно-возмущенные системы.
Рассмотрим класс систем, также относящихся к числу «жестких», для которых уже давно создана асимптотическая теория и качественная картина поведения траекторий достаточно ясна. Это позволяет четко формулировать требования к методам приближенного интегрирования и проверять, в какой мере те или иные методы интегрирования таким требованиям удовлетворяют. Итак, рассмотрим систему х=ег(х, у) (или ех=У, е=Е 'чк'1), у= р(х, у). (12) 21З е 171 ЖЕСТКИЕ СИСТЕМЫ ОДХ Здесь /, Р и их производные — величины порядка О(1); х, у— векторы размерности 1, У соответственно. Спектр вариационной матрицы (12) определяется уравнением 12'.2, — АЕ ЬУ (13) ЄР— 2Е ~ Легко показать, что жесткая часть спектра определяется спектром матрицы 2„(умноженным на Р, если, конечно, 2„не имеет собственных значений О(Е ')), а соответствующие собственные векторы Ф,.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
