Полак_и_др__Вычислительные_методы_в_химической_кинетике (972296), страница 21
Текст из файла (страница 21)
Сохране ние аддитивных интегралов движения исходной задачи на численных ре шениях обеспечивается специальным выбором аппроксимирующей ППЗ. Рассмотрим алгоритм детальнее. Пусть (7ч (ц, р(г) ) — аппроксимирующая ППЭ; д(г) = (д (г),...,де(г)) — А-мерный вектор свободных парэ метров,' р( д(г) ), ц! д(г) ) — импульсы и координаты, полученные при решении уравнений движения с потенциалом Ц, и, следовательно, гамнль- тонианом зл ди и. (ц) = и(ц„) + а х.— г=~ дц. ч» (3.102) зл зл дэи + — Х Х 2 ) =1 =1 дог да ! где а — свободный параметр, постоянный не шаге интегрирования. Система дифференциальных уравнений для Ьц и Ьрзаписываетсясле дующим образом: ди [ з ди дог ~ (=1 дцгдцу ч» (ЗЛ 03) Р»; аког + — . Ьцг(0) ЬРу(0) Ц )= 1,2,...,3л.
пгг т; векторные и матричные обозначения: Введем (т»Р~ ° т»Рэ. -" тзРзл. дч1 Рцэ - дцзл) — — — — —, О, ..., 0 (3.104) где К вЂ” гамильтониан исходной системы; М вЂ” полный момент импуль са системы; Рт — полный импульс системы. В реализованном алгоритме используется только один свободный пв раметр а(т), обеспечивающий сохранение полной энергии исходной сиота мы вдоль вычисляемой траектории в пределах заданной точности.
Сохре пение импульса обеспечивается правильным вычислением правой части уравнений движения, удовлетворительное сохранение момента импульса подтверждено численными экспериментами. Несохранение полной.энергии исходной системы вдоль вычисленной траектории р(а), ц(а) равно: ЬК = К [ р (а), ц (а) ) — К, [ р [а), ц (а) ) = и — и, (3.101) так как оператор кинетической энергии для обоих гамильтонианов один и тот же. Малая величина т»К вдоль всей траектории определяет ма лесть локальной погрешности точного решения.
Аппроксимирующая поверхность и, — квадратичная форма со свободным параметром, подбором которого можно добиться совпадения и и ил В КОНцаХ ОтрЕЗКа ИитвтрнрОВаиня. ЗтО ОбЕСПЕЧИВавт СОХраНЕНИЕ ПОП ной энергии исходной системы в точках интегрирования. Величина ЬК внутри отрезка интегрирования является критерием локальной погрешности точного решения и определяет выбор величины шага интегрирования. Пусть ц=ц»+Ьц, р р»+ Ьр, где ц» и р» — векто ры координат и импульсов в конце А-го шага интегрирования, а»»ц и ()р — соответствующие приращения на (/с+1)-м шаге. Аппроксимирую. щая поверхность имеет вид Зз =(71/0тт — якобиан системы (3.104) — матрица блочной структуры: р2 (/ О ОЧз (3.105) М О где Π— нулевая матрица (ЗлХ Зл); (7'(//()Ч'„— матрица (Зл Х Зл), в которой злемент (/,/] равен дз(/(Ч,]/ддгдд; М вЂ” диагональная матрица (Зл Х Зл) с диагональными элементами, равными 1/л1;.
Тогда система (3.105) запишется следующим образом: тт = а11 + 1т + З„чт. (3.106) Решение линейной системы дифференциальных уравнений с постоянной матрицей Зз можно выписать в явной форме: и ь ж= а / ехр(Ьзт) г/т 11+ / ехр (Зат)1/г. 1,, о о (3.107) где /1 — шаг интегрирования. Параметра получается при решении уравнения К (тт(а]] — Но = О. (3.108) Полученное решение системы (3.108), (3.106) дифференциально-алгебраи- ческих уравнений считается удовлетворительным, если в середине отрезка интегрирования оценка локальной погрешности не превосходит некоторой заданной величины: /1 /1 (Х (р( — ), Ч( — ),а) — Н (< еН,. 2 2 (3.109) л С (/1] = / ехр (Лзг) т/г, о С (2 Щ = С (Ы + С (/1) (Е + )з С (Ь) ) 6.
Ззк. 1814 (3.110) 81 Очевидно, что достаточно вычислить два интеграла в формуле (3.107), чтобы затем получать решения для любого а практически без дополнительных вычислительных затрат. Таким образом, в процесса вычислений фактически проводятся 'лишь итерации алгебраического уравнения (3.108), которое решается методом Ньютона. 8ыбор зтого метода обусловлен тем, что при а=1 аппроксимирующий потенциал является разложением исходного потенциала в ряд Тейлора до второго члена и при малых приращениях дает удовлетворительное приближение.
Следовательно, а должно быть близко к 1. Известно, что если в начальном приближении мы находимся недалеко от корня, то метод Ньютона обладает квадратичной сходимостью. Использование уравнения (3.107) дает возможность записывать необходимые в методе производные по параметру а в аналитическом виде.
а Рассмотрим процедуру вычисления интеграла / ехр(Ззт]г/г. Зтот ино теграл можно вычислить по рекурренте (6.45), обозначив с начальным приближением )ьдс "ь Ио С(Ио) = ЕИа+ + (3.111) 21 31 В случае жестких систем, т.е. при большом разбросе характеристических времен, это очень эффективная процедура. Если же шаг интегрирования лимитируется в основном величиной нелинейности дифференциального уравнения, то можно использовать формулы, в которых удается избегать большего числа перемножения матриц — наиболее трудоемкой процедуры. Для этого можно использовать рекурренту С(И+ Ио) 1 = С(И) 1+ (Е+ С (Ио)] )ьС (И) 1 с вычислением С (Ис) по формуле (3.111) .
В этой рекурренте перемножение матриц происходит только при вычислении начального приближения С(Иа). Если воспользоваться представлением интеграла в виде ряда, то можно избежать процедуры перемножения матриц, производя вычисления на векторах: Ь Из )и Ии+! ) ехр(.)ьт) г)г 1=И1+ — 1+ ...
+ 1+ .... (3.113) е 21 (3.112) При вычислении ряда (3.113) используется блочный вид матрицы )'„' и векторов 1г и 1э, что дает выигрыш в вычислительных затратах: и) о о (-и — ) )за к )за+ 1 ь (3.114) О В численной процедуре суммирование ряда (3.113) производится до тех пор, пока норма последнего вычисляемого слагаемого не становится существенно меньше нормы вычисленной суммы. Начальный шаг интегрирования выбирается из малости нормы )),)аИ1, что обеспечивает быструю сходи мость ряда. Ржчеты близких траекторий. Производные по начальным усяовиям Рассчитывая траактории с заданными начальными условиями, можно одновременно вычислять и близкие траектории, разложив решение уравнений движения в ряд Тейлора ло приращениям начальных условий. Это дает возможность получать траектории, зависящие от новых начальных 82 условий, не решая уравнения движения, что позволяет существенно увеличивать статистический материал при небольших дополнительных вычислительных затратах и зффективно заполнять начальный фазовый обьем.
Обозначим через х(хе, т) решение задачи Коши [3.97), зависящее от времени и начальных условий. Проварьируем начальные условия и представим новое решение в виде Ох х(хе+ Ьхе) = х(ха)+ — бха. [3.115) Охо Вычисляя х(хс) и оператор У=Ох!Рхе и используя формулу (3.115), можно находить близкие к х (хз, с) траектории. Рассмотрим процедуру вычисления У. Для используемых нами гладких потенциалов правая часть дифференциальных уравнений движения обладает непрерывными производными О1/Ох, позтому оператор У является решением уравнения [188): (3.118) У = 3[т) ° У, У (О) = Е, О1 где 3 (г) = — [ „[„,) — якобиан систеыы (3.97), а Š— едичная матрица. Ох Для того чтобы записать дифференциальное уравнение для У, согласованное с численным решением. необходимо использовать якобиан приближающей системы дифференциальных уравнений.
Тогда получится матрица оператора У, состоящая из производных численного решения по вектору начальных условий. Будем решать (3.116) на той же сетке, что и основную систему уравнений. Пусть хе, хы ...,хь — решение системы (3.97) в концах отрезков интегрирования. Тогда Охз ь- ~ — И У. (3.117) Ока 1=) где Ут=Рх~! Охг Матрица У~ определяется на каждом шаге интегрирования. Текущее значение х(т) на рм шаге равно х;(т) =х; ~ ьтт[т) и х;(т) удовлетворяет уравнению ' дхг — а(~+ 1з+ )(х~,)(х~(т) — х; ~).
(3.118) дт Якобиан системы (3.118) равняется 1, (Оа/Ох)т+)(х~,), тогда уравнение для Уг записывается е виде хоп ~, У1 = )(х~ ~) У~+ 1~ ~ — ~)~ У;; Уг(0) = Е. (3.119) Ох Решение системы (ЗЛ19) может быть записано в виде, удобном для итераций: ат гба~ Ут ехр[)Ь|) + ехр()Ь1))' ехр(-)т)1~ ~ — )ту~(т) гтт (3.120) е Ох где Ь, — величина рго отрезка интегрирования. 63 В расчетах можно ограничиться лишь первой итерацией уравнения (ЗЛ20) у~,') ехр(В1), так как вычисление следующих итераций чрезвычайно трудоемко и неоправданно из-за приближенного характера формулы (3.115) .
Матричную экспоненту можно вычислять как сумму ряда Згбг )яИи ехр(5И) = Е+ 3И+ + ... + + 21 л) (3.121) или с использованием рекурреяты: ехр(12И) = ехр()И) ехр()И). [ЗЛ 22) Ошибки, связанные с приближением решения первым членом ряда Тейлора и численным решением уравнения для производных по начальным условиям, оцениваются величинами несохранения аддитивных интегралов движения, а также выборочно могут быть вычислены экспериментально. В конкретных расчетах на временах порядка 10 собственных колебаний молекулы зти ошибки оказываются незначительными. Рассмотрим вычислительную процедуру получения производных по начальным условиям.
Начальные условия, как правило. варьируются только в некотором подпространства Ч фазового пространства, связанного с возбуждением лишь части нормальных колебаний, поэтому можно вычислять часть оператора У, связанную с Ч. Пусть 1 — матрица размерности (2лХ гл), столбцы которой являются базисом подпространства Ч. Тогда нужно вычислять матрицу У. (., т.е. применять оператор У к гп векторам. Вычисление оператора У на векторах матрицы 1 уменьшает объем вычйслений по сравнению с полным вычислением У. Оценки показывают, что вычисление оператора У.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
