книги / Метод продолжения решения по параметру и наилучшая параметризация в прикладной математике и механике

3.1. Особенности численного интегрирования жестких систем ОАУ

81

алгоритмов, успешно решающих жесткие системы, относится к категории неявных. Однако и здесь имеются исключения. Так, в [37, 106] показано, как можно решать жесткие системы при помощи явного метода Эйлера. Предложенный подход основан на достижениях в теории численной устойчивости разностных схем и итерационных процессов чебышевского типа. На основании этого разработана программа DUMKA.

Позднее было доказано, что для многошаговых методов значительно полезнее некоторые менее ограничительные требования к устойчивости, и были введены Понятия Л(аг)-устойчивости [118] и жесткой устойчи­ вости [93].

Численный метод называется Л(а)-устойчивым, а е (0, тг/2), если его область устойчивости включает бесконечный клин | arg(—А)| < а.

Если в предыдущих определениях речь шла только об устойчивости, то в [93] было введено более сложное понятие жесткой устойчивости,

включающее		в себя как устойчивость, так и точность приближения
к экспоненциальному фундаментальному решению уравнения (3.9).
Метод называется жестко устойчивым, если он абсолютно устойчив
в области	J?i	и точен в области J?2, где J?i : {Re(Aft) ^ D}, R i : {D <
Re(Aft) ^	a,	| Im(Aft)| < в}, D, а, в — числа.

В работе [86] была описана формальная методология решения жестких обыкновенных дифференциальных уравнений, использующая формулы дифференцирования назад. Первый пакет программ DIFSUB, использующий эти формулы, которые до этого не пользовались популяр­ ностью, был разработан Гиром и описан в [93] Этот комплекс программ позволяет также применять модифицированные методы Адамса для ре­ шения нежестких систем уравнений. Следующий комплекс назывался STIFF (см. [93]). Рассмотрим идеи, лежащие в основе данных программ.

Формулы дифференцирования назад (BDF) Л-того порядка для ре­ шения на т + 1 шаге системы уравнений (3.1) имеют вид

*-1
Уга+1 " W(jfm+li ^m+l) —^ 1 ®iVm-i ~~	(3.10)

При к > 6 методы, описываемые этими формулами, неустойчивы. Для решения системы нелинейных уравнений (3.10) используется метод Ньютона—Рафсона

Р(Ут\.\> *m+1)

(3.11)

/р(Ут+1) ^т+|)

82	Глава 3. Жесткие системы дифференциальных уравнений

где j — номер итерации,
„0>	*	л = „0) . _
ЧУп+и	W	I) = Ут+\ -		h b f ( V m + V
	7 /,,0>		*	\ -	дР(У' ^
	•Ы Ута+1>		*m + l)—		Q y

W l) - Е “ij/m-i. »=0

w_„0 )

у y«+l

Начальное приближение задается полиномиальной аппроксимацией неизвестных функций. Из формул (ЗЛО), (3.11) следует, что рассматрива­ емые методы являются методами прогноза-коррекции. Для оптимального применения методов был разработан алгоритм управления порядком ме­ тода и шагом интегрирования.

Дальнейшие усовершенствования рассматриваемых программ связа­ ны с исследованиями Хиндмарша. Разработанный под его руководством комплекс программ GEAR [113] является модификацией DIFSUB, ко­ торая связана с изменением способа определения точностных характе­ ристик решения (осуществляется контроль только локальных ошибок), способом вычисления якобиана и рядом программных изменений.

Пакет программ EPISODE [79] предназначен для решения как жест­ ких (методы дифференцирования назад), так и нежестких (методы типа Адамса) систем обыкновенных дифференциальных уравнений. Основное отличие этого вычислительного комплекса от пакета программ STIFF заключается в различных способах изменения шага интегрирования, частоте вычислений якобиана, оценке локальной ошибки.

Самая последняя, известная нам версия программ серии Гир называ­ ется LSODE [100, 99]. В ней для решения нежестких дифференциальных уравнений используется неявная формула Адамса—Мултона. Жесткие уравнения решаются при помощи формул дифференцирования назад. В итерационном процессе (3.11) применяется модифицированный метод Ньютона—Рафсона. Начальное приближение определяется по формулам, аналогичным (3.10), но явным. Якобиан вычисляется только при плохой скорости сходимости итерационного процесса.

Другим известным пакетом программ решения жестких систем, основанном на формулах дифференцирования назад является комплекс FACSIMILE [85], который предназначен для решения жестких задач, возникающих в химической кинетике.

Следует также отметить программы для решения жестких за­ дач, содержащиеся в библиотеке программ ВЦ МГУ [12] и програм­ мы RKF4RW [ПО], RAI4 [78], позволяющие автоматически определять, является ли задача жесткой, основанные на формулах Рунге—Кутта, при­ чем последняя программа использует адаптивные формулы. Упомянем здесь также программу RADAU5 [97].

3.1. Особенности численного интегрирования жестких систем ОДУ

83

Изучим особенности применения Л-преобразования к жестким си­ стемам уравнений (3.1) при использовании линейных многошаговых методов. Общая форма записи этих fc-шаговых методов такова

ft	к
агУт+г = Л	Pifm+i> Ш = 0, 1, 2, . . . ,	(3.12)
1=0	1=0

где а<, Pi — постоянные, ак ф 0, |ао| + |/30| Ф 0.

Формулы (3.12) определяют линейные соотношения между ут+{ и f m+i, i = 0, к. Для того, чтобы вычислить последовательность

приближенных значений ут, необходимо сначала каким-либо спосо­
бом получить	к начальных значений Уо» У1»- - •»2Лк—1- Если рк = 0,
то Ут+к легко	вычисляется, и метод (3.12) называется явным многоша­

говым методом. Если	рк Ф	0, то правая часть	формулы (3.12) содер­
жит fm+k - /{Ут+к,	tm+k),	и в общем случае	для определения ут+к

необходимо решать нелинейное уравнение. В этом случае метод (3.12) называется неявным многошаговым методом. По сравнению с явны­ ми методами [57] неявные являются более точными и устойчивыми

втом смысле, что при возмущении значений ут шаг интегрирования h

внеявных методах может принимать гораздо большие значения.

Ранее была показана предпочтительность неявных методов при ре­ шении жестких систем уравнений. В этом случае проблема заключается в нахождении решения уравнения

Ут+к —

1(Ут+к< tm+к) + 9т> Р к Ф 0,

(3-13)

<*к

где функция дт содержит известные величины ym+j, f m+j, J = 0, fc - 1, определенные на предыдущих шагах процесса интегрирования уравне­ ний (3.1). В [98] доказано, что если шаг интегрирования

	h<	<*к	]_	(3.14)
		Рк	L'
	где L — постоянная Липшица (||/(зм, <) - /(» 2. ОН <			~ 2/2II для всех

t € [fot^l)» то существует единственное решение ут+к уравнения (3.13), которое может быть получено с помощью итерационного процесса

Ут+к^ — ~~ / (Ут+к’ *т+к) + Ут, J —0, 1, 2, . . . .	(3.15)
а к

Важно заметить, что скорость сходимости зависит от величины h\Pk/ak\L, и тем выше, чем эта величина меньше единицы. Такое требо­ вание накладывает сильные ограничения на величину шага интегрирова­ ния h, поэтому при решении жестких систем уравнений вместо метода простой итерации (3.15) используется метод Ньютона—Рафсона (3.11).

84	Глава 3. Жесткие системы дифференциальных уравнений

Но ситуация улучшается, если применить Л-преобразование. Рас­ смотрим это на примере одного уравнения (случай n = 1 для систе-

мы (3.1)). Можно принять, что L = \-щ\ = \fiV\. После применения

к уравнению (3.1) Л-преобразования, первое уравнение полученной си­ стемы примет вид

_ f{y, t)

y/l + f 2'

Константа Липшица для этого преобразованного уравнения будет равна

Ш _ L

| 1 + / 2|3' 2 |1 + /2|3/2-

Очевидно, что при большом значении правой части f(y, t) выполня­ ется условие Lx С L, и требование (3.14) будет менее ограничительным.

Если положить, что h\ = h( 1 + / 2)1/2, то неравенство (3.14) примет вид

А <	Ок	1 + / 2
	А	L *

Это подтверждается и примерами, приведенными ниже, которые показывают, что в этом случае можно использовать модифицированный метод простой итерации.

Исследования, проведенные в главе 2, показали, что Л-преобра­ зование увеличивает область устойчивости явной схемы метода Эйлера и область неустойчивости неявной схемы, однако и после применения Л-преобразования свойство А -устойчивости метода сохраняется, что, как было установлено, очень важно при решении жестких систем уравнений.

Применение Л-преобразования к тестовой системе уравнений (2.18) в случае, когда собственные значения матрицы этой системы удовлетво­ ряют условию |ai| > |d2|, показало: если предположить, что решения у\ и з/2 одного порядка, то спектральное число обусловленности пре­ образованной системы будет меньше числа обусловленности исходной, а, как мы видели выше, эта характеристика может служить косвен­ ной мерой жесткости системы. Это же касается и размаха спектра, который уменьшается после применения Л-преобразования, что свиде­ тельствует о сближении собственных значений преобразованной системы уравнений.

Таким образом, из всего сказанного можно сделать вывод о том, что предложенное Л-преобразование смягчает жесткость исходной систе­ мы уравнений, а значит, повышает эффективность любого численного метода ее решения.

3.2. Сингулярно возмущенные уравнения

85

3.2. Сингулярно возмущенные уравнения

Уравнения, содержащие малый параметр при старшей производной, называются сингулярно возмущенными. Они образуют класс жестких систем, на котором удобно проводить теоретические исследования с це­ лью определения эффективности численных методов, предназначенных для интегрирования жестких систем уравнений. Это возможно благо­ даря достижениям в асимптотической теории [11], теории разностных схем [24] и простоте анализа качественного поведения траектории задачи.

Рассмотрим сингулярно возмущенную систему уравнений [88]

® = / ,
которую представим в форме
Г V = Y{y, z),	(3.16)
ez = Z{y,z), e
1 i = LZ(y, z),	L

Здесь Y , Z — гладкие вектор-функции, имеющие вместе со своими производными порядок малости 0(1). Размерности векторов у, z рав­ ны к и I соответственно. Введем следующие обозначения для векторов:

	* = {У, г), / =		(Y, Z).
Спектр матрицы	определяется уравнением
	Y y - АЕк	Y z	\Е\ = 0,
	L Z y	LZ^z

где Ек, Ei — единичные матрицы порядка к и I.

Очевидно, что жесткий спектр определяется спектром матрицы LZ<z и соответствующими собственными векторами, имеющими определяю­

щие компоненты z. Их компоненты у имеют порядок О ( jr). Качествен­ ная структура траектории хорошо известна и определяется поверхностью

Г = {х = {у, z} : Z(y, z) = 0},

которая	разделяет фазовое пространство на	две	части Z(y,z) > О
и Z(y,z) < 0. Система (3.16) является жесткой			только в тех точ­
ках x(t),	где спектр матрицы ZtZ устойчив,	т.е.	где действительные

части собственных значений отрицательные. На рис. 3.1 показана ти­ пичная траектория системы (3.16) для скалярных функций у и г, но та же картина будет и в общем случае. Вне малой окрестности поверх­

ности Г	фазовая скорость имеет почти горизонтальное	направление
и очень	большую величину (||®|| = 0(Ь), ||у|| = 0(1),	||£|| = 0(L)).

Точка x(t) быстро перемешается направо, если Z(y,z)> 0, или налево,

8 6	Глава 3. Жесткие системы дифференциальных уравнений

если Z(y, г) < 0. За короткий промежуток времени, О Ц ) , точка x(t)

перемещается из начального положения X Q в О ^^-окрестность по­

верхности Г. Здесь х = 0(1), и точка x(t) медленно перемещается вдоль поверхности Г. Поверхность Г разбивается на две части в соответствии со знаком Zy. устойчивая часть Г (если 2 г < 0 в скалярном случае, или*если спектр матрицы Z<z является устойчивым в общем случае) и не­ устойчивая часть (если ZtZ> 0 в скалярном случае, или если матрица Z z имеет по крайней мере одно собственное значение с положительной действительной частью). На рис. 3.1 устойчивые участки Г — это участ­ ки АВ и CD, а неустойчивый участок — BD. Система (3.16) не будет жесткой в малой окрестности неустойчивой части Г, так как здесь матрица / >х имеет большие собственные значения с положительными действительными частями.

Наиболее интересное явление проявляется около точки В или D, где Г теряет устойчивость, и одно из собственных значений матрицы Zz перемещается из левой части комплексной плоскости в правую по­ луплоскость, а точка x(t) перемещается из устойчивой части Г в точку С на другой устойчивый участок Г в течение короткого промежутка време­

ни, О • Эта часть траектории образует так называемый внутренний

слой. Затем x(t) движется согласно знаку Y вдоль поверхности Г со ско­ ростью 0(1) в точку D и так далее. Та же траектория x(t) показана на рис. 3.2. Интервалы с медленным движением разделяются коротки­

ми О пограничными или внутренними слоями быстрых перемеще­

ний со скоростями 0 (£), и кажется, что траектория становится раз­ рывной. Такого типа ситуации наблюдаются с траекториями уравнений,

3.2. Сингулярно возмущенные уравнения

87

описывающих процессы, происходящие в химической кинетике, см., например, поведение траектории в задаче об орегонаторе [92, 63].

К таким же результатам приводит и анализ простейшего сингулярно возмущенного уравнения вида [53]

£ dt = ^ У’	(3.17)
	£ > °'

Рассмотрим случай, когда вырожденное уравнение, соответствующее уравнению (3.17)

/(ж, t) = О,

имеет единственное решение ж = ж(<), и в окрестности этого решения величина d f/d y отрицательна. Последнее условие является достаточным для устойчивости решения ж = ж(<).

Характер поведения решения уравнения (3.17) следующий. Для дос­ таточно малого s касательные к интегральным кривым даже при не­ большом отклонении от функции x(t) почти параллельны оси у. И чем меньше величина е, тем быстрее осуществляется сближение интеграль­ ной кривой и решения ж(t) вырожденного уравнения.

Эта ситуация может быть описана следующим образом. У любой ин­ тегральной кривой из рассматриваемой области выделяются два участка

ссущественно различным поведением решения, причем продолжитель­ ность первого участка значительно меньше второго. Первый участок

сбыстрым изменением искомой функции отражает стремление инте­ гральной кривой к графику функции x(t) и называется пограничным слоем. На втором участке производные решения значительно меньше,

аинтегральная кривая практически совпадает с графиком ж(f). Погра­ ничный слой всегда будет иметь место, кроме случая, когда начальное условие является корнем вырожденного уравнения, т. е. уо = ж(<о). Раз­ личный характер поведения решения на обоих участках проявляется тем отчетливее, чем меньше величина параметра е.

88	. Глава 3. Жесткие системы дифференциальных уравнений

Таким образом, вне пограничного слоя для описания решения дифференциального уравнения (3.17) может быть использовано реше­ ние вырожденного уравнения. То, что даже при небольшом отклонении начальных условий от графика функции x(t) в любой его точке производ­ ная решения dy/dt резко возрастает по сравнению с производной dx/dt и определяет сложность численной реализации задач рассматриваемого типа.

В качестве первого примера рассмотрим решение задачи Коши

		§	= - Н у -		g(t)) + f t , у(о) = ю,	(зле)
где к >	1, g(t) =	10 -	(10	+ t)e~*.
В	[24] при	к =	200	эта	задача названа «искусственной	тестовой

задачей Лапидуса и Зейнфельда». Она имеет аналитическое решение

2/(0= 9(*) + Юе"*‘.	(3.19)
Здесь мы изучим более жесткий случай: к =	1000. Для численно­

го решения задачи использовался метод прогноза и коррекции первого порядка, в котором прогноз осуществлялся при помощи явной фор­ мулы Эйлера, а коррекция — неявной. Более подробно данная про­ грамма РС1 описана в главе 2. Решение с точностью 10-5 строилось на отрезке 2е[0,1].

После применения Л-преобразования задача (3.18) принимает вид

_		k (g - y )+ g
dX	y/l + (b(9 -У)+ У)2		2/(0) = 10, t(0) = 0.	(3.20)
< Л	V	j
. dX	ф	+ (Н 9 - у ) + 9)2’

Время счета этой задачи было в три раза меньше времени счета за­ дачи (3.18). Это объясняется тем, что для получения результатов, хорошо согласующихся с точным решением (3.19) для малых значений t, т.е. в пограничном слое, при решении задачи (3.18) приходилось принимать

очень малый шаг интегрирования по координате t, равный 10~5, тогда как при решении преобразованной задачи (3.20) допустимый начальный шаг по аргументу А был равен 0,1.

Вторая задача, взятая из [24], это нелинейная задача Эдсберга:

“ = “ 2ку2, 2/(0) = Ю,

(3.21)

3.2. Сингулярно возмущенные уравнения

89

имеющая точное аналитическое решение

10

y{t) =

1 + 20JW

После Л-преобразования задача (3.21) принимает вид

_	—2ку2
dX	у/ 1 4. 4k2y* ’
	у(0) = 10, <(0) = 0.
dt_ _	1
dX	у /\ + 4й2у 4 ’

Эта задача интегрировалась с точностью 1(Г10 при к — 103 на ин­ тервале t € [0, 1] при помощи программы РС1. Время счета на PC АТ 286/287 при начальном шаге интегрирования равном 0,1 составило 3 секунд. Решение иепреобразованной задачи потребовало в два раза большего времени счета.

В качестве еще одной задачи такого типа рассмотрим так называемую «специальную задачу» Далквиста [24]

e j = ( ! - % - У2, »(0) = 0,5, е = 1(Гб.

(3.22)

Хотя это нелинейное уравнение является уравнением Бернулли, но оно не имеет аналитического решения, однако вне пограничного слоя, который для данной задачи чрезвычайно мал, численное решение может быть сравнено с решением вырожденного уравнения х = 1 —t.

Решение задачи (3.22) будем отыскивать на отрезке £ € [0,1] с точ­ ностью 10~7. После применения Л-преобразования задача (3.22) прини­ мает вид

dy _	(1 - t ) y - y 2
, dX	^ e 2 + ( ( l - t ) y - y 2)2’	(3.23)
dt _	y(0) = 0,5, <(0) = 0.
	e
dA	y/e2 + ((l ~ t ) y - y 2)2

Программа PCI, описанная в главе 2, потребовала слишком большо­ го времени счета при решении этой задачи, поэтому была предпринята попытка, находясь в рамках процесса простой итерации, модернизиро­ вать его таким образом, чтобы добиться приемлемого времени счета. Модернизация заключалась в том, что итерационный процесс решения нелинейной системы уравнений

Ут+ 1 — Ут ~i~h’f{ym+Ь ^m+l)> ^ — 0» 1» 2, - - - ,

90	Глава 3. Жесткие системы дифференциальных уравнений

к которой сводится интегрирование задачи Коши для системы диффе­ ренциальных уравнений

^ г = /(у,<)> у(*о) = уо, у е м " , / =

неявным методом Эйлера, строится на основе метода Ньютона—Раф- сона, в котором учитываются только элементы, стоящие на главной диагонали матрицы Якоби. В этом случае итерационный процесс будет описываться формулой

		(*)	,	(к)
(*+1) _	(*) _	У«,т+1	Уч» ПМУт+У W l)
2/,>+1 -	У|',т+1		9fi(ym, tm)		(3.24)
		1	— Л ----- ------------
т = 0,1,2,...;			^У«,т	___
		Л = 0,1,2,...;		»=1,п,

а начальное приближение вектор-функции у, вычисленное на т + 1 шаге интегрирования, определяется явным методом Эйлера.

Принят следующий алгоритм вычислений: если итерационный про­ цесс, задаваемый формулой (2.30), не достиг заданной точности за 10 ите­ раций, то включается итерационный процесс, опысываемый форму­ лой (3.24), в котором диагональные элементы матрицы Якоби dfi/dy^m вычисляются не на каждой итерации, а только в случае, если этот итера­ ционный процесс не сходится за 10 итераций, если же и эта мера не дает результата, то только в этом случае происходит деление шага пополам.

Согласно вышеописанному алгоритму была разработана программа

SUBROUTINEРС1М(п, d, h, е, f c t, t , у, k, m) ,

написанная на алгоритмическом языке ФОРТРАН для PC типа IBM, имеющая такие же идентификаторы, что и программа РС1.

При помощи программы РС1М задача (3.23) была проинтегрирована на PC АТ 286/287 за 5 минут.

Программы РС1 и РС1М реализуют одношаговый метод перво­ го порядка точности. Такого типа методы в основном используются при больших расчетах реагирующих течений в двухмерных и трехмерных задачах [48], так как они требуют хранения в памяти минимального объема информации. Однако, многошаговые методы более высокого по­ рядка точности сходятся быстрее и являются более точными. Поэтому были разработаны еще три программы РС2, РСЗ, РС4, использующие многошаговые формулы дифференцирования назад [93] второго, третьего и четвертого порядка точности соответственно.

Согласно рекомендациям [5, 53], прогноз целесообразно вычислять по экстраполяционным формулам, использующим только значения са­ мого вектора решения в предыдущих точках.