Fedorenko-RP-Vvedenie-v-vychislitelnuyu-fiziku (810773), страница 41
Текст из файла (страница 41)
Расчету подлежит только переходной процесс, длительность которого порядка 1 с. Опишем в общих чертах метод, который использовался в расчетах и позволил проинтегрировать уравнение (8) за 20 —: 30 шагов по времени. Правда, шаг (переход от состояния в момент ! к состоянию в момент 1+ т) предполагает выполнение достаточно сложной операции — решения главной спектральной задачи. В целом этот под- 5 161 гллвн»я спскт»Ассов»я злдАчл для кглввых з»длч эоз ход оказался эффективным именно благодаря тому, что был использован очень быстрый способ нахождения главной собственной функции (многосеточный метод; см. 3 !4).
Основу метода составляет определение (разумеется, приближенное) искомого решения системы уравнений (11) из «асимптотиче- ского уравнения» ! —, = А(а)МоехР ~ Лс(С') асС' !Рс(с), а(0) = ао, о где зрс(с) — главная собственная функция спектральной задачи 14а(сН р= Лй'р; Л! — соответствующая точка спектра; Дсо — нормируюший множитель, определяемый «начальными данными» (при С= О множитель дс выбирается таким, чтобы уже началось незначительное, но заметное изменение температуры а). Опуская некоторые технические детали, стандартный шаг численного интегрирования можно представить следующими операциями.
Пусть в некоторый момент времени с уже получены значения а(с) и ~ Л,(с') асс'. Определяются коэффициенты оператора 1.(а), решается о главная спектральная задача, находятся !рс(!) и Л,(с). После этого для вычисления а(С + т) делается один шаг по явной схеме Эйлера: а(С + т) = а(С) + тАСа(С)1М» ЕхР ~ Л, о!С' Р,(С), о с+! ! пересчитывается интеграл ~ = ~ + тЛ,(с), и т.д. о о Естественно возникает вопрос о правомочности такого приближения. Некоторые физически достаточно убедительные аргументы дает следующее рассуждение. Подставим используемую конструкцию в исходные уравнения задачи. Разумеется, они не будут выполнены: возникает некоторая «невязка».
Если она очень мала относительно входящих в уравнение членов, то это в известной мере оправдает вышеизлос женный подход. Обозначая для удобства Л(с) = ~ Л, сй', вычисляем о с »с (а) оо ""с(с) дс елс!'! Л изр .1 ~ ' ца)„р у ехЯ й" пгивлижвнныв ми«оды вычислительной виэики 1ч. п Итак, погрешность метода зависит от соотношения й'д~р~/д~ и членов, входящих в Ь(а) 1рг (В этн выражения входят коэффициенты, погрешности измерения которых часто не так уж малы; кроме того, сама формулировка исходной задачи, которую мы принимаем за истину, в действительности основана на пренебрежении некоторыми относительно малыми величинами.) Заметим, что собственная функция ~р, определена с точностью до нормировки и этим фактором тоже можно разумно распорядиться с целью уменьшения погрешности.
Можно обыграть дза обстоятельства: величины о, и о существенно разные (скорости быстрых и медленных нейтронов связаны соотношением о~ ж 1Озоз). Следовательно, важнее получать медленное изменение второй компоненты 1рг Другое обстоятельство, которое можно было предвидеть на основе опыта расчетов в этой области: компоненты главной собственной функции ~р, (в зависимости от пространственных координат, которые в принятых обозначениях опущены) — просто устроенные функции «колоколообразной» формы. С течением времени (при изменении а(г)) их форма меняется не очень существенно.
Заметно меняется лишь отношение их амплитуд (т.е. соотношение между потоками быстрых и медленных нейтронов в холодном и горячем реакторе). Из сказанного выше следует рецепт нормировки: постоянной полагается гильбертова норма второй компоненты ~рг При численном решении задачи оценивалась и величина идгыр,!дп оиа фактически оказалась достаточно малой по сравнению с другими членами (около 1 —: 2 %). 5 17. Жесткие системы обыкновенных днфференцнапьных уравнений В пятидесятых годах при решении задач Коши для систем, описывающих кинетику реагирующих друг с другом химических веществ, вычислители столкнулись с крайне неприятным явлением. Расчеты проводились с помощью хорошо отработанных программ, в которых использовались методы Рунге — Кутты и надежные алгоритмы автоматического выбора шага.
Эти алгоритмы очень быстро вырабатывали шаг численного интегрирования, столь малый, что часто ие было никакой возможности рассчитать процесс на требуемом для приложений отрезке времени, даже используя наиболее мощные ЭВМ той эпохи. Визуальный анализ правых частей, казалось, не давал оснований для каких-то опасений. Типичная система уравнений химической кинетики выглядела (технические подробности опускаем) примерно так: х'="', А!х~+~Х" А! х'х", й у, й 1, 2, ..., М, (1) » !7! жксткив системы одт где А,', А'.х — константы, характеризующие скорости протекания тех нли иных реакций.
Бросалась в глаза, правда, существенная разница в их значениях: они отличались друг от друга часто на много порядков. В то же время исключить какие-то «малые» члены из (1), оставив только самые большие, было нельзя, ориентируясь лишь на значения А. Существенными были и концентрации разных веществ х'! они могли очень сильно изменяться с течением времени. Затрачивая значительное машинное время, удавалось получать начальные отрезки траекторий и провести анализ ситуации. Он выявил следующую характерную картину. В начале процесса происходит сильное изменение х(() и выбираемый программой шаг численного интегрирования вполне разумен: он очень мал, но так и должно быть для аккуратного инте! рирования столь быстро меняющихся функций.
Через небольшое времв г харак- х((! тер траектории резко меняется, она стано- и х х вится гладкой, медленно меняющейся, но программа этого «не замечает» и выбирает такой же малый шаг, Попытки «подска- ~Ч(»с~х((! зать» программе выбор существенно боль- Я$~ ч(~~~(~х шего шага, согласованного с гладкостью Хс(((7 ~ф~х(П решения, немедленно приводили к вычислительной катастрофе. В з 5, 7 специально указывалось„что прн оценке вычислитель- РиС. 2! ной сложности задачи Коши для системы обыкновенных' дифференциальных уравнений х = У(х) существенны два фактора: строение поля траекторий в окрестности интегрируемой траектории и свойства матрицы /„(х). Анализ поля направлений таких систем, получивших название «жестких», дал характерную картину, качественно представленную рис, 21. Траектория х(() состоит из короткого участка быстрого ее изменения (так называемого «пограничного слоя») и длительного участка очень медленной ее эволюции (иногда его называют «квазистационарным режимом»), Основные трудности связаны именно с расчетом последнего.
Пограничный слой интегрируется очень малым шагом, но он настолько краток, что число шагов интегрирования вполне приемлемо. На рис. 21 при помощи «микроскопа» с последовательно увеличивающимся разрешением показана структура поля направлений в окрестности к((). Сначала видны траектории, отвесно падающие на х(г). При следующем увеличении видно, что, приближаясь к х((), они поворачивают, стремясь двигаться «параллельно» х(Г). И лишь при еще большем увеличении видна стандартная картина практически параллельных линий.
числвнныв мвтолы млтвмлтичвскоа оизлки 1ч и зш Если из точки х(г) траектории сдвинуться по касательной в точку х' = х(Г) + тх(1), то, хотя расстояние х(~ + т) — х' = О(тз) ничтожно, У(х ) не имеет ничего общего с х(» + с), направление г'(х') скорее напоминаег перпендикуляр к л, 1я1 траектории х(г). То же самое поо Лг лУчаетсл и пРи опРеДелении х' отрезком ряда Тейлора из трех- четырех н более членов при том о l о ЯсЛ значения т, которое хотелось бы использовать для численного интегрирования квазистационарною режима.
Рис. 22 Анализ матрицы У„(х) в ок- рестности траектории также привел к специфической картине, которую мы примем за основу при следующем формальном определении жесткой системы. Определение 1. Задачу Коши х = У(х), х(0) = х,, 0 н г я Т, х Е Ял (2) будем называть жесткой„ если спектр матрицы У„(х) достаточно четко делится на две части (рис. 22). Жесткий спектр, Собственные значения и векторы обозначим Лс(х) и Ф,(х) (1 = 1, 2,..., 2).
Для жесткого спектра выполняются условия ИеЛ~(х) я — 1., ~1шЛ,(х)~ < )йе Л,(х)~ (3) Мягкий спекпср. Собственные значения и векторы обозначим Л (х)„и у (х) (у = 1, 2, ..., У). Для мягкого спектра выполняются условия ~Л1(х)/ я Е сЯ.. (4) Время интегрирования Т является средним относительно 1 и очень большим относительно Тл 1Т равно, например, 10, 20, 30, з 1.Т может быть порядка 10з, 10в и больше. Отношение И1 называют показателем жесткости системы. В приложениях встречаются ситуации, когда 2П равно 10~, 10~, 10м. Будем считать, что 1=!. Стандартные методы интегрярования'требуют шага т', малого в том смысле, что т'аУ„И' м1.