Полак_и_др__Вычислительные_методы_в_химической_кинетике (972296), страница 35
Текст из файла (страница 35)
Отменить режим вычисления якобизма. 6. Провести Ш-факторизацию матрицы Р. 7, Методом Ньютона решить нелинейную систему уравнений относитель но вектор-столбца в = у — у (О), и л (и — у„) — И|(1) [Р(ц, т ) — ул ] =0 (5.21) на каждой итерации решаатся линейная система относительно о +,1 Рц = Рц — (о — у )- — [Р(ц,т ] — у (О) )' (О) лн.1 м т л ю м' и и )1 ни = у(О), т = О, 1, ..., М. л (5.22) ч Е ((т) = -(~ ' 1+ П (у — у(О) ). л 2, т г и и и-1 л — 1 (5.23) ПуетЬ тт=[И11, ..., Ил]т — ВЕКтср ВЕСОВ ОШИбОК дпя КаждОй ИЗ КОМПОНЕНТ вектора решения, вычислить 0 = П1ак((Ел (д) [)тт().
( ( где Š— Ея компонента вектора Е (д). л. л Проведенный шаг считается успешным в случае, если (5.24) (5.25) 0 < еМз. е Ошибка может контролироваться на интервале времени з или на шаге интегрирования вне зависимости от его величинь1 (в зтом случае з = й): параметры е и з задаются пользователем. Если контроль точности удовлетворителен, перейти к шагу 16, в про тинном случае — к шагу 12. 12. Восстановить массив с„). Перейти к шагу 13. 13.
Вычислить () — козффициент изменения шага по формуле ч = (еЫ1) з) '1(ч"), (526) О е перейти к шагу 14. 14. Вычислить новое значение шага интегрирования ))'= 1) т). Перейти к и1агу 15. 15. Перенормировать массив си 1 на новое знечение шага путем уьь 138 8. Если итерации корректора сошлись за установленное число ша" гов М (как правило, принимается М= 3), то перейти к шагу 11, в яроги(ь ном случае — к шагу 9.
9. Восстановить массив д„). Перейти к шагу 10. 10. Если якобиан был вычислен на текущем шаге интегрирования, то уменьшить значение т) и перейти к шагу 15, в противном случае — уста повить режим вычисления якобиана и перейти к шагу 3. 11. Провести контроль ошибки интегр)фования: вычислить вектор НсжЕНИЯ МатРИЦЫ Уя 1 На ДИаГОНаЛЬНУЮ МатРИЦУ 1 О т)1 перейти к вегу 1. 1В.
Вычислить 2. '= 2(~) + (у — у<о) ) 1. Перейти к шагу 17. 17. Вычислить коэФфициенты изменения шага для случаев метода по. рядка д-1 т) = (елтэо ) метода порядка д (е))(эо ) 1/(а+1) а и ееетода порядка а+1 (еьта11 ) 1/(а+2) а+1 а Выбрать т) =п1ак)(т)~, т(,, т) ~ и очередное значение поРяДка мето. да д, соответствующее максимайьному значению коэффициента измене ния 1вага. Перейти к шагу 18. 18. Если д' д-1, то вычислить коэффициенты полинома а — 2 о(к) = кз П (х+$ ), (5.31) 1=1 х- (т — т„)тл.
$, = (г„— т„,)тл, обозначаемые через <(1. В соответствующих колонках 7 " 2, 3, ..., д-1 марива 2 разместит величины. 1/Юу(а)Ьу1. Если а' а+1, то в колонке а+1 массива 2 разместить нулевой вектор, так как ла+1л(а' ~)1 (д+Ц 1 = О. [5.32) л 19. Перенормировать матрицу 2„на новое значение шага интегрироая ния умножением на диагональную матрицу 1 О т) т) (6.33) 139 ~/г I Рис.
3.2 ех /дат 1 Рис. 3.3 Рис. 3.1. Аппроксимации решвния ыадтиь. наго уравнения х = Л х(х(0) = х„; Л > О) 1 — точное решение! 2 — ненвнан схема; 3 — нвнан схема! 4 — полунвчвнан схема Ет гззг Рис. й2. Оценка величины оцюибки рецтения модельного уравнения по нвявной схема путем сравнения рецгений с магам 6 и 2й Рис. б.3. Оценка величины оыибки паникин модапьнопт уравнения по разности между результатами счете схем резного яорндка точности Рис. бхб Оценка величины оыибки ревенин модельного уравнения па разности между результатами амтв по явной и неявной схемам ной системы обыкновенных дифференциальных уравнений у = 0 'ЛОу, уе =0 'ха, О 1 -1 — 1 1О-4 (5.37) 1О-4 О Зта система получена в результате применения ортогонального преабре.
зования (матрица О) в системе независимых линейных обыкновенных дифференциальных уравнений вида х = Л х, х,(О) = хат, 1' 1,2,...,5. (6.38) !а) С]чевидно, что система (5.37) обладает жесткостью - 1(Г4 и имеет фундаментальное решение, отвечающее положительному корню +1, Вы бором начального значения хе, можно подавить на начальном участке интервала интегрирования компоненту решения, соответствующую положи тельному корню ха, ехр! Т) . При расчете этой задачи по программе ВТ] Рг, являющейся реализацией метода Гира, при соответствующем выборе начальных условий шаг интегрирования выходит за границы жесткой устои чивости метода и, как следствие этого, приводит к неверному численному решению.
Чтобы избвхать этих недостатков, вообще говоря присущих всем ревностным методам, необходимо проводить оценку величин положь тельных корней якобиана системы. Этот вопрос будет обсуждаться в следующем разделе. Во многих задачах химической кинетики процессы взрывного характера происходят с большими временными задержками, математически это и означает медленное появление компоненты решения, соответствующего положительному ' корню якобиана. Для эффективного решения задач, описывающих такого рода процессы, был предложен особый класс числе+ ных методов, рассмотренных в следующем разделе. Э. МЕТОДЫ ЛОКАЛЬНОЙ ЛИНЕАРИЭАЦИИ Проблема жесткости, как следует из определения, заключена уже в линейной задаче. Исследование системы обыкновенных дифференциаль. ных уравнений на жесткость производится с помощью анализа линеаризо ванной задачи.
Поэтому в последнее время широкое распространение получили методы, в которых упор делается на решение или на анализ линеаризо ванной системы [1, 95, 140, 161, 314, 393]. Методы решения систем нелинейных обыкновенных дифференциальных уравнений, использующие идею локальной линеаризации, имеют два аспекта: 1) локальная линеаризация, т.е. способ приближения нелинейной системы ОДУ на шаге интегрирования линейной, и оценка величины возникаю щей при этом ошибки; 2) выбор способа решения линейной системы. Изложим метод локальной линеаризации, следуя работам [44, 140].
Пусть ду/дг 1 (у), у (гэ) = уе — исходная задача Кошм Запишем при ращение о(Ы для у на шаге интегрирования: и(Ы = у(г+Ы вЂ” у(г). (5.39) тогда для приращения ц (Ы можно записать систему диффцоенциапьных уравнений дц — 1(у(г) ) + Ац(Л) + р (у(г), и(Л) ), о (0) = О, (5.40) д)г а1~ где А — якобиан системы; ер(у(г),ц(й)] =1(у(г+Ы)— ау]т=„(г] — 1 (у (г! ) — Аи (Ы вЂ” нелинейная часть. Тогда решение уравнения (5.40) может быть записано в виде л а ц(Л) = ]' ехр(Аз)сй 1(у(г) ) + ехр(АТ) [ ехр( — Аз)ч~ (у, в (й) ) г(з. о о (5.41] Таким образом, возникает проблема вычисления интеграла от матрич. ной экспоненты и интеграла от нелинейной части.
142 Интеграл от матричной экспоненты вычисляется по рекурренте: ь)+ ь ь, ь, ехр(Аз)тй = ( ехр(Аз)тй+ехр(АЬ,) ) ехр(дз)дз. е о о (5.42) Обозначим С(Ы = 1 ехр(Аз)дз,тогда о С(Ьг+Ьт) = С(Ь1) + акр(АЬ!) С(Ьг). (5.43) Используя представление матричной экспоненты и интеграла от матриь ной экспоненты в виде матричных рядов, можно получить простую формулу, связывающую эти две матрицы: ехр[Аг) = Е+ АС(г). (5.44) Тогда, используя соотношение (5.44), можно записать рекурренту (5.42) в виде С(Ь~ +Ьт) = С (Ь~ ) + (Е+ С[Ь~ ) А) С[Ьт).
(5.45) Дпя того чтобы досчитать до времени, соизмеримого с временем медленных процессов [ Хз;яЬ! -1, т.е. Ьт = 10 'а, необходимо сделать всего лишь порядка 30 удвоений [2за 10' ). Рассмотрим просадуру решения нелинейного уравнения (5.41[. Алгеб раическое уравнение (5.41) можно решать по явнойсхеме, положив з = О, тогда ц (Ы = С(Ы1(у(т) ), (5.47) либо по неявной схеме о(Ь) =С(Ь] [1(у(г)) +.р (у, о(Ь))!. (5.48) Величина локальной погрешности решения на шаге интегрирования может быть оценена как разность между решением уравнения (5.41) по явной и неявной схемам.
Уравнение (5.41) можно также решать с помощью апгазоксимации д (у, ц(з) ) полиномом по степеням з. При этом необходима вычислять соответствующие моменты от матричной а экспоненты ) з~ехрАзгй. Эти вычисления можно вести по рекурренте, исгюльзуя формулу интегрирования по частям: а Ьз+! А а ,)з» ехр [Аь ) с(з = — ехр АЬ вЂ” — 1з+ ' ехр (Аз ) т(з. (5.49) е Ь+ 1 Ь+1 о 143 Ясно, что предложенный метод решения линейной системы жестко устойчивый, так как любую линейную систему он решает точно. Пробле ма числа шагов при существенном разбросе характерных времен тоже не возникает.
Положим, что жесткость системы составляет 10 порядков, тогда, выбирая начальный шаг интегрирования из условия [ Х„,ззЬ| ) ( 1. приближая интеграл от матричной экспоненты с помощью ряда с относи тельно небольшим числом членов и воспользовавшись процедурой удвое ния, получаем С(2Ы = С(Ь) + (Е+С(ЫА) С(Ы. (5.51) дрЪ-~3 ар бц Š— Š— С(И) — ) ~С(й1(у ) = С(й — С(й1(у ). (554) и до р' и а и Заметим, что решение уравнения (5.41) по неявной схеме с оценкой локальной погрешности, например по скорости сходимости итераций, имеет тот же порядок точности и дает одинаковую величину выбираемого шага интегрирования по сравнению с методом оценки локальной погрешности по разности между решениями по явной и неявной схемам.
Это связано с тем, что и в том и в другом случае оценка зависит от одной и той же матрицы ((д,р/дц) С(И) ). Следуя изложению работы [44), опишем один нз возможных алгоритмов решения системы обыкновенных дифференциальных уравнений методом линеаризации. 1. Полагаем г = О, бу = О. Переходим к выполнению п. 2. 2. Вычисляем значение якобиана в точке г.'Переходим к выполнению п.З. 3. Вычисляеи начальное приближение для С(И), используя формулу Азлт1 С(й = И Е+ — ь — ~; И))А ()<10 з. 2! 3! (5.55) Устанавливаем режим "выбор шага".
Переходим к выполнению и. 4. 4. Систему (5.41) решаем по неявной схеме (5.48), делая лишь один шаг итераций. Требуем малости добавки к решению по явной схеме. Это позволяет избежать трудоемких итераций и. как показывает вычисли тельная практика, хорошо оценить величину локальной погрешности бц (й.
Если такая оценка величины бц удовлетворяет задаваемому в качестве начальных данных пдэаметру точности, переходим к выполнению и. 5, !44 Попытаемся оценить полную ошибку метода линеаризации. Для зтого запишем приращение решения на л-м шаге интегрирования в оператор ном видв ц(й = 2(И)1(у(г)), (5.50) где 2(И) — оператор, определяющийся тем или иным способом реше. ния нелинейной задачи. В случае решения по явной схеме 2(И) = С(И), а,р 1— а при решении по неявной 2(И) = Š— С(И) — ~ С(й.
Тогда до~ — у„+ ц„— у„+ г„(И1(у ). Если Еи и 2„+ 62„, где 2„— оператор точного уравнения, то дпя погрешности решения, связанной с погрешностью оператора, справедливо при ближенное уравнение бу„, = (Е+ 2„Ат„)бу +62 1(у ). (5.52) Первый член в формуле (5.52) определяет закон накопления ошиб ки, т.е. связан с устойчивостью численного метода, а второй определяет локальную погрешность на шаге интегрирования. Если решать задачу по явной схеме (5.47), а неявную использовать для оценки локальной погрешности на шаге, то уравнание для ошибки (5.52) будет иметь вид бу =ехрАИ бу +бц (5.53) где л Зр(ехр(Аг)) = Х е~н, ~=! (5.56) где 8р(ехр(Аг)) — суыма диагональных элементов матрицы ехр(Аг).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
