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

2.3. Алгоритмы, программы, примеры

71

Если искать решение этого уравнения, удовлетворяющее начальному условию

у(0) = О,

(2.36)

то при * = 0 правая часть уравнения (2.35) обращается в бесконечность, и большинство численных методов интегрирования ОДУ потерпело бы неудачу. Будем искать решение задачи (2.35), (2.36) для у € [О, L]. После применения к этой задачи A-преобразования получаем задачу

Для оценки точности численного решения найдем аналитическое решение задачи, которое может быть представлено в параметрическом виде

х = R(t - sin t), у = R(l - cos t), Д = у .

Исключая параметр t, получаем

Поэтому для оценки точности воспользуемся выражением

у

Д= — - 1 + cos

R

где х* и у* — численные значения.

Задача (2.37) с одинаковым успехом интегрировалась как при по­ мощи программы РС1, так и с использованием метода Рунге—Кутга 4-го порядка при L = 2. При шаге интегрирования Ад = 0,5 • 10-2 методом Рунге—Кутта наибольшее значение Д = 0,8 • 10~4 достигалось при А = 1,25. Конечной точке интегрирования соответствовало значение аргумента А = 3,75.

Программа DLSODE при решении задачи (2.35) на первом же шаге попадает в аварийную ситуацию, что естественно. Интегрирование же при помощи этой программы задачи (2.37) происходит успешно.

Пример 2.3. Улитка Паскаля. Рассмотрим задачу Коши

= р (х, У) dx Q(x,y)'

у(о ) = h

(2.38)

72	Глава 2. Задача Коши для дифференциальных уравнений
где
	Р(х, у) = lx2 - S(a, у)(2х -	a), Q(x, у) = y(2S(x, у) - 12),
	S(x, у) =	х1 + у2 - ах,

а и I — заданные числа. Эта задача имеет замкнутую интегральную кри­ вую, которая называется улиткой Паскаля. На рис. 2.10 показан вид этой кривой в случае, когда параметры а

и I удовлетворяют неравенствам

а < I < 2а.

Улитка Паскаля описывается ура­ внением

S2( x , y ) - i 2(x2 + y2) = 0,

поэтому точность вычислений оцени­ валась по формуле

Д = |S2(®*, у*) - 12(х*2 + у*2)|, где х*, у* — численные величины.

Обычные численные методы, в том числе и программа DLSODE, не в состоянии построить всю интегральную кривую задачи (2.38). Но задача может быть решена как при помощи программ DLSODE, PCI, метода Рунге—Кутга, так и других численных методов, если к ней применить Л-преобразование, после которого задача примет вид

dy Р(х,у)	dx Q(x, у)		./ЛЧ_ .,	_/лЧ я
<2Л				(2.39)
где	Z(x, у) =	yJp1{x,y) + Q1(x,y).
Задача (2.39) решалась при		а = 1,	I = 1,5 при помощи программ
DLSODE, PCI, метода Рунге—Кутта 4-го порядка. Для получения всей
замкнутой кривой параметр Л изменялся от 0 до				10,5. При интегриро­

вании задачи методом Рунге—Кутга с шагом h\ = 0,5 • 10-2 величина Д

не превосходила значения 0,2 • 10-4.

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

нения Эйлера [38]
sin в sinip	0	cos ip
= | sin в cos <p	0	- sin ip	(2.40)
cos в	1	0

2.3. Алгоритмы, программы, примеры

73

определяющие связь между компонентами вектора угловой скорости ш ( ш \ о/з(£)) в подвижной системе координат, связанной с объек­ том, и производными по времени от углов Эйлера ip, <р, 0, где 1р — это угол прецессии, <р — угол собственного вращения и 0 — угол нутации.

Для определения закона изменения углов Эйлера разрешим систему уравнений (2.40) относительно производных, тогда получим [38]

dip S(ip)	dip	d$
1>= Ш ’	« “ “ з - Э Д 'в * .	(24l)

S(ip) = o/j sin у/ + о$ cos <p.

Очевидно, что при численном интегрировании данной системы уравнений возникнут трудности при угле нутации 0 = 0, так как в этом случае правые части первых двух уравнений обращаются в бесконечность. Однако проблема исчезает, если к системе (2.41) применить Л-преобра- зование. В этом случае эта система может быть представлена в виде

' dip	S(<p)	sign (sin в),
dX	Q
dip	u>3 sin 0 —S{ip) cos 0
dX =		sign (sin 0),
		Q
dj0		(2.42)
	o/j cos ip - W2 sin ip
dX ~		| sin0|,
		Q
dt	| sin 0\
dX ~	Q	’
		1/2

Q = [sin2 0 + S 2(y>)+sin2 в(ш[ cos ip-uii sin <р)2+(шз sin 0-S(ip) cos 0)2j ,

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

Система уравнений (2.42) решалась численно при помощи програм­ мы РС1 для и}[ = —1, W2 = l, 0/3 = 0 и начальных условиях: при А = О, ip = ip = 0 = t = O. Полученное численное решение определяло связь между углами Эйлера вида ip = —ip.

Заметим, что те же проблемы могут быть разрешены и при чи­ сленном интегрировании кинематических уравнений для самолетных углов [36]. Эти уравнения имеют особенность при значении угла танга­ жа, равном тг/2.

Глава 3

Жесткие системы обыкновенных дифференциальных уравнений

Опыт численного интегрирования задач Коши для ОДУ показал, что среди них необходимо выделить так называемые жесткие уравнения, которые требуют специальных численных методов. В этой главе будет рассмотрено влияние Л-преобразования на такого рода системы ОДУ и на методы их интегрирования.

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

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

Такое явление называется жесткостью, а системы обыкновенных дифференциальных уравнений, моделирующие процессы такого типа, называются жесткими системами уравнений.

Отметим, что жесткость задачи — это свойство математической мо­ дели, и она не связана с используемым численным методом. Жесткость задачи является математическим отражением того факта, что в соот­ ветствующем физическом объекте протекают процессы с существенно различными скоростями.

Характерным для всех жестких систем является такое поведение решения задачи Коши, при котором компоненты первого типа претерпе­ вают либо быстрые начальные изменения, либо значительные изменения на некотором участке наблюдения (пограничном слое).

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

75

Необходимость выделения данного вида уравнений в отдельный класс вызвана трудностями их численного интегрирования классичес­ кими методами, например, явными одношаговыми и многошаговыми методами. Выяснилось, что малый шаг интегрирования, используе­ мый для воспроизведения быстропротекающих процессов в пограничном слое, не может быть увеличен вне пограничного слоя, хотя производ­ ные становятся существенно меньше. Даже незначительное превышение некоторой величины шага, определяемой данным методом и решаемым уравнением, приводит к резкому возрастанию погрешности. В самом деле, как показано, например, в работах [5, 53, 57, 55], для того, чтобы обеспечить абсолютную устойчивость численного решения задачи Коши

~ = /(У ,< ) .	У(<0) = У0,	t € [ t 0,T],	(3.1)
У (У1) • • • >Ул) )	f ~ i f 1} •' ‘ } fn)	» уо —(УЮ) • • • I УпО)	>

необходимо использовать такой шаг интегрирования А, при котором каждое из, вообще говоря, комплексных значений Hi = AAj (г = 1,п), где Aj — собственное значение матрицы Якоби df/dy, лежало бы внутри области абсолютной устойчивости. Таким образом, для методов с ограниченной областью устойчивости длина шага лимитируется поряд­ ком величины наименьшей временной постоянной системы. Поскольку интервал интегрирования может во много раз ее превосходить, то не­ обходимое число шагов интегрирования может оказаться чрезвычайно большим.

Для того, чтобы устранить это ограничение, были предложены численные методы, см., например, [93, 80, 112], допускающие значи­ тельное увеличение шага интегрирования вне пограничного слоя, однако и в настоящее время проблема численного решения жестких уравнений является актуальной [4, 97, 37, 60, 106, 88].

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

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

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

Следствием этого, как правило, является с одной стороны относительно высокий порядок системы уравнений, а с другой — ее жесткость.

Следует также отметить и скрытые формы проявления жесткос­ ти. Так, большой класс гладких оптимизационных конечномерных за­ дач [S3, 119] с трудом поддается решению традиционными методами изза овражного рельефа поверхностей уровня. Изучение этого явления [S3] показало, что трудности связаны с жесткостью системы дифференци­ альных уравнений, описывающих траекторию наискорейшего спуска, и поэтому естественно строить овражно-ориентированные алгоритмы оптимизации с учетом этого факта.

Заметим, что жесткими могут оказаться задачи, описываемые урав­ нениями в частных производных, если их решение сводится к системе обыкновенных дифференциальных уравнений. Так, в [80] показано, что жесткой будет задача теплопроводности для стержня, а в [97] указано на жесткость одномерной задачи диффузии и задачи движения упругой балки.

Выделим следующие характерные свойства жестких линейных сис­ тем [53]:

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

2.Собственные числа Л,- матрицы Якоби J = d f/d y полностью опре­ деляют характер решения.

3.Из жесткости неоднородной системы

^ = J y + y ( 0 , y e R " , <€[<0, П

(3.2)

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

4.Вне пограничного слоя между компонентами вектор-функции y(t) можно установить линейные соотношения, число которых равно количеству быстро осциллирующих частных решений системы (3.2), т. е. вне пограничного слоя решение жесткой системы может быть описано решением системы меньшей размерности, уже не являю­ щейся жесткой (см., например, [93]).

Несмотря на большое число публикаций по данной проблеме, до сих пор не существует общепринятой концепции жестких систем [5, 60, 88]. Более того, нет даже общепринятого определения жесткости.

Так, в [112] понятие жесткости связывается с понятием устойчиво­ сти, и под жесткой задачей понимается задача, не имеющая неустой­ чивую компоненту решения (у матрицы Якоби системы уравнений (3.1)

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

77

отсутствуют собственные значения с большими положительными дей­ ствительными частями), но имеющая несколько очень устойчивых ком­ понент (по крайней мере одной компоненте соответствует собственное значение с большой отрицательной величиной действительной части).

В [105, 55] придерживаются следующего определения.

	Задача Коши (3.1)	называется	жесткой на	некотором интервале
	I С [*о. Г], если для t € I
	1)	Re А«(0 < 0, * = 1,т»,
		шах	| Re A,-(i)|	(3.3)
	2) S(t) = l<i<n		»
				1,
		min	|ReAj(t)|

где Л< — собственные значения матрицы Якоби J = df/dy, в кото­ рую подставлено решение задачи (3.1) у = y(t) при значении аргумента равного t.

В [53] проводится подробное критическое обсуждение наиболее известных определений жестких систем уравнений, при этом сами авторы отдают предпочтение следующему определению:

Система обыкновенных дифференциальных уравнений

§ « / ( » , * ) ,

(3 4)

Г С It х К$, 4 = {0 < t < 00}

называется жесткой на отрезке изменения независимой переменной [а,Ь], принадлежащем интервалу существования ее решений, если при любом векторе начальных значений (зд, (о) € Г и на любом отрезке [to, *0+2’] С [а, Ь] найдутся такие числа г, L, N , удовлетворяющие неравенствам

т < Ъ - а	(3.5)
и	(У,*)€ Г,

где p(df/dу) — максимальный модуль собственных чисел матрицы Якоби (спектральный радиус), || • || — принятая норма матрицы, что справедливы неравенства

dyн	L
dt	, max ,lw(*)l. Л = 1,7»,
	N »6|io, to+T]

to + r ^ t ^ t 0 + T, N > \ .

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

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

Авторы считают, что важным моментом введенного определения является неразрывная связь понятия жесткости системы (3.4) с вели­ чиной промежутка наблюдения решения [а,Ь], заложенная в неравен­ стве (3.5). Если жесткую на [а, Ь] систему рассматривать лишь на проме­ жутке [а, с] С [а, 6], включающем только пограничный слой т = с - а, то на [а, с] ее нельзя считать жесткой, так как никакого различия в характере поведения решения не наблюдается.

Понятие жесткой системы уравнений в [37, 106] вводится следую­

щим образом.		J(y,f) = \\dji!dyi\\ (i,j	___
Пусть J(t)	=		= l,n ) — матрица Якоби
функции f(y,t),	где y(t) — решение задачи (3.1).
В окрестности		точки (од, to) системе	уравнений (3.1) можно по­

ставить в соответствие систему линейных дифференциальных уравнений вида

= J {к)(у ~ од) + /(од, *),

(3.6)

которая является линеаризацией по у исходной системы уравнений (3.1) в окрестности (од,to). Предполагается, что решение системы уравне­ ний (З.б) с начальными условиями yfo) = од аппроксимирует с задан­ ной точностью е > 0 при при to ^ t < т решение задачи Коши (3.1), где 0 < т < Т, т — зависящая от to величина. Пусть А,-(<) (* = 1, п) — собственные значения матрицы J(t), которые все являются действи­ тельными. Им соответствует полная система собственных векторов &(t) (t = 1, та). Спектр матрицы J дает верную информацию о качествен­ ном поведении решения системы (3.1). Задача Коши (3.1) называется в [37, 106] жесткой, если существует не зависящая от t постоянная С > 0 такая, что

A ,(t)< C для	te [to ,T ], » = I7 n ,	(3.7)
и при этих условиях имеет место
AT= AT(t)=max(-Aj(f)) > 0,	5 = M ( T - t0) > 1, te [f0,r].	(3.8)

В работах [60, 88] задача Коши (3.1) называется жесткой, если спектр матрицы Якоби системы уравнений (3.1) может быть разделен на две существенно различные части:

1) жесткий спектр, где собственные значения и собственные векторы обозначаются через Лi(t) и Ф*(<) (* = 1, . . . , / ) соответственно. Для жесткого спектра выполняются условия

Re Aj(t) < - L < 0,

| Im Aj(f)| < | Re A,-(f)|, * = T7;

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

79

2)мягкий спектр, где собственные значения и собственные векторы обозначаются через Aj(t), ipj(t) (j = I, J; J = n - I ) . Для этой части спектра справедливы условия

|АД*)|<1, 3 = ^ 7 ,

и если отношение L/1 принимает большое значение.

При этом, так как матрица Якоби df/dy, вообще говоря, зависит от t, то следует писать I(t),J(t).

В[93,53, 105, 57, 80,112,97] приводится обзор работ как по методам

иалгоритмам, применяемым при решении жестких систем [93, 53, 105, 57, 97], так и по вычислительным программам этого направления [80, 112, 97].

Термин «жесткая система», по-видимому, впервые был введен в ра­ боте [86] при анализе решения задачи

у' - -50(у-cost),	у(0) = 0,
методами Адамса и Рунге—Кутга. Была также показана целесообразность
использования неявных методов при	численном решении уравнений

такого типа.

Впервые явление неустойчивости с ограничениями на шаги интегри­ рования для гиперболических систем уравнений в частных производных изучено в знаменитой статье [83].

Однако, применительно к жестким системам обыкновенных диф­ ференциальных уравнений именно в [87] указано на численную не­ устойчивость как на причину затруднений, введены основные понятия и определения.

Для того чтобы обеспечить устойчивость численного решения задачи (3.1), необходимо использовать такой шаг интегрирования Л, при ко­

тором каждое из, вообще говоря, комплексных значений Щ =	ЛА,-
___	d f
(i = 1, те), где А,- — собственные значения матрицы Якоби J(t) =	—-,
	ду

лежало бы внутри области абсолютной устойчивости. Таким образом, для методов с ограниченной областью устойчивости длина шага лими­ тируется величиной наибольшего собственного значения матрицы J. Поскольку интервал интегрирования может значительно превосходить величину 1/ min |А<|, необходимое число шагов интегрирования может

оказаться сравнимым с коэффициентом s(t), задаваемым формулой (3.3). В самом деле, на каждом шаге интегрирования система уравне­ ний (3.1) может бьггь преобразована к линейной системе (3.2) с матри­ цей Якоби. Если предположить, что все собственные значения матрицы Якоби различны, то произведя соответствующую замену переменных,

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

матрицу можно привести к диагональному виду с собственными значе­ ниями Aj, J = 1,п, на главной диагонали, и система уравнений (3.2) примет вид

Щ = XiVi +Ф)> * = М*-

Поэтому изучение свойств методов принято проводить на тестовом уравнении [87, 10S]

(3.9)

где y(t) — скалярная функция. Решение уравнения (3.9) асимптотически устойчиво, если Re Л < 0, неустойчиво, если Re А > 0, и устойчиво, если Re А = 0.

При использовании численных методов дифференциальное уравне­ ние обычно преобразуется в разностное. Разностное уравнение асимпто­ тически устойчиво, если все корни его характеристического уравнения по модулю меньше единицы, неустойчиво, если хотя бы один корень по модулю больше единицы, и устойчиво, если есть некратные корни, по модулю равные единице, а остальные меньше единицы. Очевидно, методы следует использовать при таком шаге интегрирования в зави­ симости от А, когда наблюдается соответствие по всем видам устой­ чивости тестового дифференциального и соответствующего разностного уравнения.

Множество значений ЛА, удовлетворяющих условию асимптотичес­ кой устойчивости решения разностного уравнения численного метода, возникающего при интегрировании тестового урввнения (3.9), называ­ ется областью устойчивости (абсолютной устойчивости) метода в ком­ плексной плоскости ЛА.

В главе 2 было показано, что областью устойчивости явного метода Эйлера является внутренность круга единичного радиуса с центром в точке ft Im А = 0, ft Re А = -1 .

Анализ области устойчивости традиционных алгоритмов Рунге— Кутта и Адамса показывает, что они не годятся для решения жестких систем, так как использование малых значений h требует больших за­ трат на вычисления. В связи с этим в [87] было предложено создавать и применять для жестких уравнений такие методы, у которых область устойчивости при решении тестового уравнения (3.9) содержала бы всю левую полуплоскость плоскости ЛА, т. е. при Re А < 0 решение соответ­ ствующего разностного уравнения было бы асимптотически устойчиво при любом положительном Л. Эти методы называются А-устойчивыми.

В [87] также доказано, что явный линейный многошаговый метод не может быть A-устойчивым и что порядок неявного многошагового метода не может превосходить двух. Поэтому подавляющее большинство