Полак_и_др__Вычислительные_методы_в_химической_кинетике (972296), страница 34
Текст из файла (страница 34)
Опишем теперь общие принципы построения разностных схем для систем обыкновенных лифлреренциальных уравнений 1-го порядка [174). 1. Разложение правой части системы лифференциальны х уравнений в ряд Тейлора. Запишем разложение У(тч+1) В РЯД ТайЛОРа: лт Ьл У (Гп+~ ) =У(Г~)+Ау (Г„» Ч вЂ” уч(Г„) + + у(Л) (Г ) 21 р( Обозначим через Ь (г. у, л) разность (у(т„ч ~ ) — у(г„))IЬ, тогда й (б у, л) = у (г„) + — у' (г„') + ... + — у (гя) ° г и ю л (л) (б.7) 21 Аппроксймируя Ь(бу,й) той или иной функцией, например ограничиваясь в разложении членами порядка р и заменяя точные функции приближенными значениями у„, получим численный метод, использующий производные р-го порядка.
Если мы ограничимся первым членом по й разложения в ряд Тейлора и пренебрежем членами высших порядков. то получим явную схему Эйлера: У (Гч ~) = У (т„) +Л У' (Г„) = У„+л((Г„. У„). (Е.б) Исходя из такого подхода, можно получить методы Рунге — Кутта. Для етого Ь (г, у, л) аппроксимируется функцией Ф (б у) в разных точках отрезка интегрирования. Это позволяет избежать вычисления производных и обойтись лишь вычислением функции т(г.
У) . Такие ьвтоды называются явными методами Рунге — Купа. Если не каждом шаге интегрирования функция Ф (г. У) вычисляется р раз, то метод называет р-этапным явным методом Рунге — Кутта. Эти методы не требую~ вычислений производных 1(д у) высших порядков, что является, как правило, трудоемкой процедурой. Прн использовании зтих методов на нужно вычислять дополнительные начальные значения, и существует возможность гибко менять шаг интегрирования. 134 2. Л и н ей н ы е м н о га ш а го в ые методы отличаются от методов Рунге-Кутта тем, что для вычисления последующих значений уяг, нужно использовать ранее вычисленные значения у„, у„.
с, у„ Идея получения формул численного решения состоит в том, что задача Коши записывается в интегральной форме: с+» у(т+д) — у(г) = Х у(В, у(И) (в, (5. 9) Затем 1 [т, у (т) ) заменяется интерполяциомным полинамом, проходящим через точки, в которых значения у уже вычислены или будут вычисляться. Если, например, приближенно вычислять интеграл по прямоугольникам с высотой 1(с+гг, у О+И), то получим неявную схему Эйлера: (5.10) В самон общем случае лри использовании такого подхода получаются формулы Адамса: » » ьу+ ° =ах()(п=012 (5.11) 1с, = Ь [Š— (гА[ ' 1(у„), (сг й [Е - ЬА[ '1 (у„— (сг ), (5.12) 135 Из формулы видно, что, для того чтобы получить численное решение в точках с„ь, необходимо иметь значения численного решения в точках ге, гс, ..., т» с.
Линейные многошаговые методы могут быть явными, если (1» = 0 и для вычисления у„+» необходимо иметь значения численного решения в предьщущих точках, и неявными, когда в правую часть уравнения входит величина 1(т„ыи у„+„) и для вычисления у„+» при ходится решать нелинейное по у„г» алгебраическое уравнение. Несмотря на то, что явные схемы представляются в вычислительном плане менее трудоемкими, на практике чаще используются неявные разностные методы, ибо они в большей степени обладают свойством устойчивости, что позволяет выбирать больший по сравнению с явными схемами шаг интегрирования. Нелинейное алгебраическое уравнение, возникающее в неявном методе, реагается методом итераций. Скорость сходимости итераций во многом зависит от того, насколько близко к корню выбрано начальное приближение ус+„. Начальное приближение выбираетея с помощью явного многошагового метода, который.
называется предсказги вающим, а уточнение решения происходит с помощью итерационного процесса для решения нелинеймого уравнения неявной линейной много шаговой схемы. Такой метод носит название предиктор-корректор [174) . Широкое применение для решения прямой кинетической задачи нашли методы Розенброка [389[. Основная идея зтих методов состоит во везде нии якобиана системы в разностную схему Рунге-Кутта. Артемьевым и Демидовым [8 — 10) предложен ряд схем такого типа. Рассмотрим предложенный в работах [8, 9[ метод 4-го порядка.
Численное решение задачи Коши для автономных систем ОДУ осуществляется по формулам: 13 1 2 у = Х + — 1сс + — 1сг — 2)сг + — 1сы "+г " 6 6 3 3 (сэ = Д [Š— ЬА] '1 У + — (т~ + — (гз), н В 3 19 1 )с4 и [Š— дА] 1~ун + (г~ + (гз (тз) В 24 6 Здесь у +„у — векторы решенийна (и+1)-м и л-м шагах интегрирования соответственно; Š— единичная ыатрицэ; А — якобиан системы, вычисленный в точке угг Для того чтобы решать неавтономные системы ОДУ, необходимо добавить еще одно уравнение: г = 1, г (О) = гэ.
В [110] приведена программная реализация этого метода. Выбор шага интегрирования н контроль точности организованы следующим образом: по формулам (5.12) осуществляется решение с обычным и удвоенным шагом (обозначения у и узе соответственно) . Для каждой компонен- А гы вектора решения вычисляется погрешность Ь = (У" — У™[, ( =1,2,...,Д(, (5. 13) И вЂ” число уравнений; если у Ф О, вычисляется относительная погреш ь ) ность л[ (5.14) Шаг интегрирования считается успешным, если выполняется ассимптоти ческая оценка [12] 5 = птах б, ( 1бе, ! 5.15) 1 где е — заданная относительная погрешность на шаге интегрирования.
В противном случае полагается и' = Ы2, н шаг повторяется снова. В случае б С е/2 величина шага удваивается. Для увеличения скорости вычиспе. ний вместо обращения матрицы [Š— пд] осуществляется ее (.О-факториза цня [183], после чего последовательно решаются четыре системы линейных алгебраических уравнений относительно векторов (г,, )гт, (тз, х4.
Кроме описанного метода 4-го порядка, Артемьевым и Демидовым [10] предложен еще ряд методов такого же типа, в том числе метод с переменными порядком и шагом интегрирования. В последние десять лет широкое распространение получил алгоритм численного интегрирования жестких систем ОДУ, предложенный Гиром [263, 264]. Алгоритм основан на использовании линейных многошаговых методов, удовлетворяющих требованиям жесткой устойчивости [263].
При вычислении преднктора применяется алгоритм Норсика [352], исполь зуюцэчй интерполяционный полипом для вычисленных в предыдущих тоэ ках значений вектора решения. За счет этого легко осуществляется пере ход к новому шагу интегрирования, что обычно представляет определенные трудности при традиционной реализации многошаговых методов. Выну ление корректора, как правило, осуществляется методом Ньютона, причем для матрицы [Š— ()ейА] (Š— единичная матрица, й — текущее значение шага, [)э — параметр метода, А — якобнан системы) используется (.(э'.разложение, что, как известно [183], позволяет наиболее эффективно решать возникающие линейные системы алгебраических уравнений.
При решении задачи Коши методом Гира в каждой точке выбирается оптимальный порядок метода„обеспечивающий наибольший возможный шаг интегрирования. 134 л„(т„) = у„в Р(у„,т„), (5.16) л л — з ул-П и вычислении ул, такого, что построение этого полинома возможно. Будем предполагать, что имеется решение в точках и — 1, п — 2, п-З, ... и необходимо построить решение в точке и. Алгоритм Гира можнц представить в виде следующей последовательности шагов. 1. Вычислить значение предиктора 2л" по формуле ~о (5.17! Здесь 2л ~ — вектор Норсика, представляющий собой массив размер ности Ф х (д+1) (д — текущее значение порядка метода, Ф вЂ” число уравнений): 2„, = [у,,бу,, ...,Ьчу'"),1Ч)], (5.16) А(д) — треугольная матрица Паскаля с элементами О, 1<1 П 1,1 = 0,1,...,д.
1>1 /! (1-1')! А(д) = (5.19! 2. Вычислить х = (т — гл)1й, $~ = (т„— г„,)1п, 1= 1,.,ч, и поли ч ном Лл(х) = (! (1+х1Е,.). Коэффициенты попинома обозначаются через (=1 11, т.е. ч Л(х)= Х 1х1, )=о 1 (5.20) показано, что и составляют вектор-строку ! = [1о„1ы ...,1ч]. В [226] ч 1ол1, 1~= Х с,'.
1 =11$ы....~. ч ч' 3. Вычислить значение функции Р (у, т„). (о) 4. Если установлен режим вычисления якобизма, то перейти к шагу 5, в противном случае — к шагу 7. 137 Программная реализация метода Гира несколько раз модифицировалась, Первой реализацией алгоритма является подпрограмма 0(РВОВ (1971 г.! [263]. Затем были разработаны более совершенные пакеты программ 6ЕАЙ (1972 г.) [262] и ЕР!600Е (1975 г.! [226]. В пакетах программ В Г(РР .(1977 г.! [67] и 6ЕРЯВ (1973 г.) [283], являющихся переработками пакета ПЕАЯ, предусмотрена возможность работы с ленточными матрицами, что особенно важно лри решении уравнений в частных произ водных методом прямых.
Пакет программ ВОЕЧЕВ (19ВО г.) [62] ориентирован на работу с системами уравнений, имеющими разреженные матрицы Якоби. При изложении алгоритма Гира будем следовать в основном работе [226], содержащей, на наш взгляд, наиболее удачную методику контроля ошибки вычислений. Основная стратегия состоит в построении интерпопирующего попинома ал степени д ипи меньшей (и — номер шага интегри рования), удовлетворяющего ч + 2 условиям: 5. Вычислить якобиан системы А в точке ул,, гл, вычислить ма трицу Р = Е-(Ы1)А.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
