Полак_и_др__Вычислительные_методы_в_химической_кинетике (972296), страница 22
Текст из файла (страница 22)
Е занимает примерно то же время, что и расчет основной системы. Это дает возможность существенно увеличить количество вычисляемых траекторий. Рассмотрим схему алгоритма расчета траектории и производных по начальным условиям. 1. Ввод начальных данных Вводятся начальные значения импульсов, координат и времени. Если рассчитываются производныа по начальным условиям, то вводится матрица базиса подпространства, в котором происходит вариация начальных условий.
— 1.. 2. Вычисления правой части уравнения для яр и ра щений координат и импульсов. В точке (р,ц) вычисляются якобиан и правая часть уравнения (3.103). Если этот пункт выполняется первый раз, то шаг интегрирования выбирается из малости нормы матрицы ((,)И ((, в противном случае шаг интегрирования на изменяется. л Ь 3. Вычисление интегралов )ехр()т)г(т 1, и )ехр()т)дт 1т. о о Интегралы вычисляются с помощью ряда (3.113) с использованием блочной структуры матрицы д и векторов 1~ и 1т. Затем ищется решение уравнения (3.108) методом Ньютона с начальным приближением корня и = 1. Если ньютоновские итерации сходятся, т.е.
получено решение с шагом И, то необходимо оценивать локальную погрешность решения — переход к и. 4. В случае, когда ньютоновские итерации расходятся, шаг интегрирования делится пополам. Если решение 'для этого шага уже было найдено, осуществляется переход к и. 5, если же нет то и. 3 выполнябтся с новым шагом интегрирования. 84 4.Оценка локальной погрешности решения. Проверяется выполнение условия (3.109).
Если условие (3.109) на оценку локальной погрешности выполнено и с данной правой частью уравнения (3.103) не было деления шага, то шаг интегрирования увеличивается вдвое. после чего происходит переход к и. 3. Если же условие (3.109) выполнено, но деление шага интегрирования происходило, т.е. при удвоенном шаге либо не сходились итерации уравнения [3.108), либо не выполнялось условие (3.109), то осуществляется переход к и. Б. В случае, если условие [3.109) не выполнено, шаг интегрирования делится пополаы и осуществляется вычисление решения с новым шагом интегрирования. т.е.
осуществляется переход к и. 3. Б.Пересчет векторов координат и импульсов в конце шага интегрирования г: г+л; ц:=ц+Ьц(Ь); р:=р+ +Ьр(И). Если рассчитываются производные по начальным условиям. то вычисляются их значения в конце шага интегрирования: у = ехр ()л) у. (3.123) Вычисление производится с помощью ряда (3.121) аналогично вычислению Ь интеграла )' ехр Ут) г(т Е Проверяется критерий окончания счета. Если е следует дальше продолжить счет, то необходимо перейти к и.
3. В противном случае счет прекращается. Метод опробован на ряде задач. Рассчитаны траектории движения трех- атомной молекулы Й,О и двух пятиатомных молекул СНРэ и СНэРЕ Сходимость и точность метода проверялись с уменьшением погрешности локальной аппроксимации и просчетом траектории из конечной точки в начальную с обратным временем. Проверка показала устойчивость и надежность работы программы.
С помощью численных экспериментов установлено, что сохранение полной энергии системы атомов до третьего знака (е "10Гз) дает воэможность удовлетворительно сохранять полный момент и выдерживать точность решения при расчете 10 медленных колебаний наиболее трудной для вычисления модели молекулы СНэР(, для которой разброс характеристических времен колебаний составляет около трех порядков. В работе (66) проводилось сравнение различных процедур численного решения обыкновенных дифференциальных уравнений. На основании сравнения программ по разныы критериям эффективности показано, что процедура Кутта — Мерсона (см.
(324) ) имеет оптимальные характеристики при расчете колебательных систеы. Поэтому было проведено сравнение с программой, реализующей метод Кутта — Мерсона. Критерием сравнения являлось время, необходимое для расчета траектории при сохранении алдитивных интегралов движения в пределах заданной точности. Сравнения показали, что при расчетах траекторий трехатоыных молекул и небольших траекторий пятиатомных молекул без существенного разброса характеристических времен программы обладают сравнимым быстродействием.
Так, для расчета траектории движения молекулы Нз О на временах порядка 60 колебаний молекулы этим программам требуется около 20 с ыашинного времени ЭВМ БЭСМ-6, для расчета молекулы СНРз примерно на таком же числе колебаний — порядка 3 мин, а для расчета траектории совместно с оператором производных по начальным условиям — 8-12 мин. Для вычисления траектории молекулы СНРз на временах порядка 10 медленных колебаний по предложенной программе требовалось примерно 20 мин машинного времени на ЭВМ БЭСМ-6. Пои одинаковом шаге интегрирова- 65 ния программа Кутга-Маревна на тех же временах жизни молекулы дает накопление ошибок, приводящее к существенному несохранению полной энергии молекулы.
Предложенный алгоритм оказывается наиболее эффективным при расчетах длинных по времени траекторий движения много- атомных молекул. Процедура поиска точки перевала поверхности потенциалыюй энергии Рассмотрим численную процедуру поиска точки перевала ППЭ. Пусть дана поверхность У(ц), где ц = (сы сз, -., 4„) -л-мерный вектор. Области, в которых Э И Эщ = О, определяют точки равновесия на поверхности У(ц) . Этим точкам соответствуют минимумы функционала ф ~ 2 [3.124) Задача состоит в том, чтобы найти точку перевала поверхности У(п), являющуюся точкой неустойчивого равновесия. Для поиска точки перевала воспользуемся геометрическими сообра.
жениями. Структура ППЭ, перевал которой мы ищем, характеризуется наличием двух локальных минимумов, соответствующих конфигурации стабильной молекулы и продуктам реакции. Перевал разделяет эти два минимума. Надо найти такую точку на поверхности У(ц).
которая лежит в области притяжения минимума функционала (3.124), определяющего перевал поверхности У(ц). Проведем луч, соединяющий области двух минимумов. Точку на луче, в которой поверхность достигает максимума, выбираем в качестве начальной точки спуска, если квадратичная форма, аппроксимирующая У(ц) в этой точке, знакопеременная. Знакопеременность квадратичной формы, по-видимому, должна указывать на то, что точка находится в области притяжения перевала. Минимизация функционала (3.124) должна вестись методом, приводящим к ближайшему" минимуму, т.е.
траектория спуска не должна уходить из области "притяжения'* минимума Ф, связанного с перевалом поверхности У(ц) . Таким свойством обладают процедуры поиска минимуьза функционала на траекториях решения градиентной системы дифференциальных уравнений (ц —, ц(О) це, (3.126) <(г Оп Эти траектории ведут к минимуму, поскольку (3.126) т.е. продвижение по траектории ведет к уменьшению функционала до тех пор, пока не будет достигнута точка минимума, в которой ВФ вЂ” = О. (3.127) Рц Градиентная система дифференциального спуска для функционала (3.124) имеет вид дц аи Ои с(г (79' Оц (3.128) Система дифференциальных уравнений (3.128) решается численно с использованием метода локальной линеаризации [140] по .процедуре, предложенной в работе [21[.
Очевидно, что на каждом шаге линеаризации движение происходит не по траектории наискорейшего спуска, а по близкой к ней траектории. Поэтому при решении линейной системы дифференциальных уравнений можно следить не за аппроксимацией правой части исходного уравнения дифференциального спуска, а лишь за убыванием функционала. Вычислительная практика показала, что предложенная процедура выборе начальных точак спуска и последующее решение градиентной системы обеспечивают быстрое и надежное определение точек перевала ППЭ. Требавалось -3-5 мин машинного времени на ЭВМ БЭСМ-8 для поиска перевала 15-мерной поверхности, т.е.
поверхности потенциальной энергии пятиатомной молекулы. Генерация случайных чисел При исследованиях, проводимых с помощью метода Монте-Карло, определенное значение имеет вопрос о способе генерации случайных чисел. В подавляющем большинстве расчетов, выполненных методом классических траекторий, для генерации равномерно распределенных случайных чисел используется линейный конгруэнтный метод [78].
В этом методе последовательность случайных чисел получается из соотношения Х„+, = (аХ„+ с) пкх) 1т), (3.129) 87 гдв п>0; с =Х„lгп — л.е равномерно распределенное на интервале (0,1) и н я случайное число; т обычно определяется из соотношения 2 ) гл; )Ч вЂ” количество разрядов в ячейке используемой ЭВМ; а и с — некоторые констан. ты. Длина периода такого датчика случайных чисел не больше гл. За несколько последних лет было показано, что многие используемые датчики случайных чисел, основанные на линейном конгруэнтном методе, не являются вполне удовлетворительными, т.е. не выдерживают проверки, производимой с помощью специальных тестов (критерий )(', критерий Колмогорова — Смирнова, тест частот, тест серий, покертест, спектральный тест и т.д.! [78].
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
