Бабенко - Основы численного анализа (947491), страница 90
Текст из файла (страница 90)
Численное рмиеиие оодачи Коши ВЕЛИЧИНа рит1 дОЛжНа бЫтЬ МОдИфИцИрОВаНа: ВЫЧИСнястея ВЕЛИЧИНа п1„1 = р„1 + —,.„(р„— й„) и затем вычисляется величина т„т1 112 е =,)(Хот1, Гни.1). ФОРМУЛа КОРРЕКЦИИ ПРИМЕТ ВИД 1го~-1 = (9хо и — з —, М(гнм, 1 -~- 2е -- е 1)~, и окончательный результат находится по формуле 9 ' + =" - — — (Р— 1ы+1). Описанный метод предложен Р.Хемыинго11. Эти, казалось бы, незначительные усовершенствования повышают порядок точности метода и раз ьясняют читателю, почему мы выдвинули требование выполнения неравенства С > 0 ( < О). Контролировать величину шага можно по величине бои1 = рие1— — йи,1, ДЛЯ ЭТОГО ЗаДаЮт НЕКОТОРУЮ ПОГРЕШНОСТЬ Е.
ЯВЛЯЮП1УЮСЯ ВХОД- ным параметром алгоритма, и следят за выполнимостью неравенства до 1 < Е. ЕСЛИ ЗтО НвранвиетВО ВЬШОЛНяЕтея, тО ПОЛуЧЕННЫй рвэуЛЬ- тат считается точным; если бо 1 > е,. то шаг уменьшают вдвое. Если же би» 1 < С 1е, то шаг удванва1от. Конкретное значение константы С определяется эмпирически (С - 50). Многошаговые методы с переменным шагом н способы управления нм см. в [204., 216). Рассмотрим, что изменится в последнем алгоритме, если перейти от одного уравнения к системе о19/11х =- г'(х, у)., где у = (у1, ..., у1)' искомый вектор-столбец, Г(х, й) = (71(х, у),...
..., 71(х, у)) правая часть. Естественно, что формулы (19) — (22) теперь НМЕЮт МЕСТО ДЛЯ КажДОЙ КОМПОНЕНТЫ ВЕКТОРОВ Ри+1, й„+1, Х„1, ПРИЧЕМ, если считается дья компонента, то х'„—.... )' (х„, х„). Величину д„т1 теперь следует вычислять по формуле д -1.=- ~~',пз р,. -1 — ~1, где а весовые коэффициенты, определяемыо, как правило, из эмпирических соображений. Как мы уже отмечали, величина 1пага должна быть такой, чтобы выполнялось условие (11о), которое в случае системы примет впд 11,,1е (х, у)~ << 1. (х, у) й Ю, где, какая-либо из норм матрицы уи(х, р). Естоствепно, что последнее условие не охватывает всех представляющихся случаев, и для других уравнений может быть существенным требование 11 ~ (х, 9)( << 1, Збе Нссл„олька замечаний о численнолл рсшсььи задачи Коши 441 и вообще нельзя высказать никаких суждений о величинах остаточных членов (20), (22) в зависимости от величины 1т~ (ь(х, у) ~.
Поэтому судить о приемлемости величины шага можно на основании сравнения величин рз,ттэл ~з,ь-з-л (1 =1: 2 .,1) Среди неявных разностных формул упомянем еще метод ФДН (см. [204]), основанный на формулах дифференцирования назад. Нередко его называют методом Гира по имени автора, опубликовавшего первую компьютерную програмлзу, в которой были реализованы форлтулы этого метода порядков 1 — 6. Метод ФДН обычно используют для решения жс;сткнх систем (20от1 й' 3. Несколько замечаний о численном решении задачи Коши в экстремальных случаях 1.
Жесткие системы. Некоторые задачи химической кинетикзз приводят к довольно сложным вычислительным проблемам. Уравнения, описывающие поведение концентрации хнллнческих веществ в смеси, где происхо,тнт химическая реакция, имеют вид т1уьтах = зт(у), л = 1, 2,..., 1, у = (ул,..., у~), (1) где ),(у) — чаще всего многочлены от ул, ...., И. Уравнения Д,(у) = 0 = 1, 2, ..., Ц определяют равновесные концентрации компонентов смеси.
Коэффилшенты многочленов Ь(у) связаны с константами скоростей химических реакций и могут в некоторой естественной системе единиц иметь огромные значения. Если в уравнении (1) ул,..., ря т — концент1запни, а 1Š— температура смеси, то функции (,(у) .- многочлены только от ун ..., ут л, коэффициенты которых — функции от уь которые могут изменяться на много порядков. Как следствие, система (1) такова, что матрица Ге(у), где Г(у) = (Гз(у); ., (т(у)), имеет огромную норму. Собственные значения Л (у) могут раЗличаться на много пОрядкОв, причем в Задачах хнмичаской кинетики спектр матрицы 1'„(у) имеет следующий характер: имеется несколько собсэ венных значений Л (у) с очень большими отрнпательными вещественными частями, а остальная часть спектра содержит собственные значения Л„, для которых Ве Лз может иметь любой знак, причем (Л, ьс 1.
Такие системы принято называть жесткими. Решение задачи Коши с шагом Ь, удовлетворшошим установленному ограничению на шаг 6~Ге(у) ~ << 1 нли 6 плах~Л,(у) ~ << 1, абсолютно не реально для существующих ЭВМ. Даже если бы было возможно организовать расчет с многократной точностью (что сняло бы много теоретических проблем), то это потребовало бы огромного времени для решения задачи.
Предлагалось много различных формальных определений понятия жесткости снстелты. но ни одно из них нельзя нрнзнать удовлетворительным. Чтобы лучше понять трудности задачи и сформулировать сущн~ють проблемы, рассмотрим болев простой случай, хорошо теоретически изученный, но имеющий характерные черты, присущие задаче Коши для уравнентий химической кинЕтики. Прилепе.
Рассмотрим задачу Коши для системы т1х)т11 = 1(х, у), ет1у(сХ = у(х, у), (2) Глгша 7. Численное ргласние задачи Коши Рнс, 1. Поведение ннтеграль ной кривой системы (2) где е > Π— малый параметр. Вместо того, чтобы рассматривать систему общего вида, рассмотрим частный скучай, на котором можно увидеть все основные особенности нашей задачи. Пусть 7'(х, у) =- у, д(х,. у) = д — (1/3)у — х. Это известная система Ван-дер-Поля. Обычно вместо системы рассматривают уравнение второго порядка, н его-то н называют уравнением Ван-дер-Палл. Нетрудно подсчитать, что у матрицы В(1, 'д),11)(х, д) одно из собствегпгых значений Лг х ', а второе Лз м 1.
Но Л1 может иметь произвольный знак. Если мы рассмотрим поле направлешпг д для системы (2), то увидим, что, грубо говоря, оно устроено следующим образом. Всюду вне кривой У = (х, у: д(т, у) = 01 оно , д) вертикальнО, на самой кривой гориЗонтально н только в очень узкой переходной области — окрестности кривой У вЂ” происходит перехОд от вертикальноге направления к горизонтальному. Сама кривая У разбивается в на|нем случае па три компоненты (рис. 1): одна из них дуга Уа, лежащая в полосе (х, у: --1 < д < 1), н две дополнительные неограниченные дуги У1 и Ую У1 — —. УС (х, д: у < — 1), У = У О (х, у: у > 1). Нетрудно увидеть, что в окрестности компоненты Уе поле устроено так, что почти вертикальные стрелки направлены от Уе, в то время как на компонентах Уг н Уз стрелки направлены к этим компонентам.
Таким образом, Уо — неустойчивая компонента, а Уц Уз — устойчивые. Заметим, что стационарная точка у системы единственна и лежит на неустойчивой компоненте (точка (О, О)). Точки А( — 2/3, — 1), В(273, 1) — граничные то гки компоненты Уа — играют в дальнейшем существенную роль. Опишем поведение траектории, начинающейся в точке М(хо, уо). Она будет почти вертикальной прямой, пока не встретит одну из компонент У~ или Уз, например компоненту Уь Дойдя до Ум траектория резко изменит направление (пересечет Уг в точке згг); перейдя на другую сторону Ум интегральная кривая пойдет вдоль Уц оставаясь под ней, пока не дойдет до точки А.
Здесь начинается неустойчивая, отталкивающая часть кривой У, и интегральная кривая пройдет от точки А почти вертикально до встречи с компонентой У, которуго она пересечет, резко изменив направление, и далыпе пойдет вдоль Ую оставаясь над ней, пока не дойдет до точки В, где она, сорвавшись почти вертикально, пройдет до встречи с кривой Уг н пересечет ее (далыпе события будут разворачиваться по описанной схеме). В результате установится периодический режим — колебания, известные под наименованием релаксациаииьы. Отметим характерные особегпюсти интегральных кривых.
Каждая из них будет содержать почти вертикальные участки, затем резко изменяет направление, пересекаясь с У, н буквально вне е-окрестностн точки пересечения траектория дальше следует вдаль У. Назовем почти вертикальный участок траектории и ту ее часть, которая лежит в малой окрестности точки пересечения с У, где траектория резко изменяет направление, релаксациаинмм у ~апакам. 3 3.
Несколько замечаний о численном рсшшти задачи Ко1ли 443 Это — хотя и не строгое определение, но достаточно наглядное и недвусмысленное. Для з ого чтобы численно определить релаксационные участки траектории, необходимо интегрировать систему с шагом й « щ нспачьзуя описгяшые выше алгоритмы численного решения задачи Коши. Что касается остальных участков траектории, которые не являются релаксационными, где траектория — гладкая кривая и близка к л, то естественно возникает желание иметь такую процедуру интегрирования на этих участках, когда шаг не зависит от е и может быть много болыпе ". П В уравнениях химической кинетики структура траекторий подобна тсшько что описанной, но возникновение периодических релаксационных колебаний не обязательно.
На траекториях будут релаксационные участки, а сама траектория будет стремиться к устойчивым стационарным точкам, Поэтому для жестких систем задача ях интегрирования ставится так: нужно правильно определять нерелаксапионные участки траекторий, интегрируя с болыпим шагом и допуская возможность того, что релаксационные участки будут определены со значительной погрешностью. Ясно, что не всякий алгоритм, описанный выше, подходит для этой цели. В частности, один из возможных алгоритмов основан на разностном уравнении (2.1б). В нашем случае это уравнение запишется в виде 2 л (3) где в = (в„г,..., сы)' — приближенное значение вектора и =. (у~, ...., р~)' в точке к = к„, Г = (~п, б) — вектор правых частей. Уравнение (3) — нелинейное уравнение относительно в„э ь и его ре|пепие может быть найдено с помощью метода Ньютона — благо мы всегда будем иметь хорошее начальное приближение, Единственное, что нам надлежит проверить, — это то, что дифферецпиал отображения хщю — (л72) Г (вщы ) является невырождепным оператором.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
