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

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

91

При разработке программы РС2 первый прогноз вычислялся по фор­

муле

!/£,+! — 2l/m	У т - 1)
а коррекция по формуле
У т +1 = Г (4Ут - У т -1 +	2 /l / ( j £ +1, *m+l))-	(3.25)

При удвоении шага использовались результаты, вычисленные ранее,

апри делении пополам решение в точке т о - 1/2 вычислялось по формуле

Ут -1 + Ут

Ут -1 /2 = -------2-------'

Впрограмме же РСЗ первый прогноз вычислялся по формуле

У т + l У т - 2 ~ 3j/m-l "I" ^Ут>

а коррекция по формуле

с	18j/m ~ 9Ут-1 + 2Ут-2 +	^m +l)
Ут+I —	п	(3.26)

Недостающие значения у при делении шага пополам подсчитыва­ лись как

У т - 1/2 —	3У т "Ь б У т -1	У т —2
	g	>

апри удвоении шага как

Ут - 3 = Ут —3j/m-l + 3t/m_2-

Впрограмме РС4 первый прогноз вычислялся по формуле

у£,-М = 4У™~ 6Ут-1 + 4Ут-2 - 4ут -3,

(3.27)

а коррекция

г48J/ш “ 36уга_| -1- 16ут _2 —Зут _3 + 12ft/(£ +„ <m+l)

Ут+1	-. (3.28)
	25

Недостающие значения при делении шага пополам подсчитывались согласно выражениям

5Ут+	15yw_l ~ 5ут-2 + У т - 3	(3.29)
У т - 1 /2 =	16
~ У т +	9ут - | + Чут-2 ~ Ут - 3	(3.30)
У т - 3/2	16

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

При удвоении шага новое значение ут- 2 не вычисляется — оно хранится в памяти ЭВМ, а ут-з подсчитывается по формуле

Ут-3 “ $Ут + 16j/m_i/2 15ут _] + 5?/т _2.

(3.31)

Остановимся подробнее на том, как получены формулы построения, например, последней программы РС4.

Формула коррекции (3.28) взята из [93].

Экстраполяционная формула прогноза (3.27) получена как результат решения следующей системы уравнений

(3.32)

Здесь производные вычисляются в точке ут, а значение параметра а принималось равным а = 1.

Формулы (3.29), (3.30) получим, если в системе (3.32) положить а —-0,5, а = -1,5, соответственно.

Наконец, формулу (3.31) можно получить, разрешая соотношение (3.29) относительно ут~з-

Стартовой процедурой для программ РС2, РСЗ, РС4 является усо­ вершенствованный метод Эйлера.

Помимо этих программ были также разработаны программы РС2М, РСЗМ, РС4М, которые построены по принципу программы РС1М, описанной выше, т. е. итерационный процесс организуется на основе метода Ньютона—Рафсона, в котором учитываются только элементы, стоящие на главной диагонали матрицы Якоби.

Так, итерационный процесс в программе РС2М строится на основе формулы (3.25) в виде

У{,т+1 У{,т+1

то = 0, 1,

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

93

Итерационный процесс в программе РСЗМ организуется в соответ­ ствии с формулой (3.26) и имеет вид

„(*+«> =		„(*) _
s f i + 1 “ Т Г (18У«> “	9 V i,m - 1 + 2y,,m_2 +		6	.(*)
				W l )
,	6	dfi (ym, tm)
1 ----- n— ---------
	11	&Уг,т
m = 0,1,2, ... ,	k =	0,1,2, ... ,	* =	l,n.

Наконец, в программе PC4M итерационный процесс в соответствии с формулой (3.28) описывается выражением

J * + i) = u «			_
2/,>+i - 2^(48yi,m - 36y.,m_i +		\ 6 y iim - 2	-	+ ^ h f i ( y ^ + l , t m+ l )
1 _	n	h d fi(y™> ^m)
	25	d y i ^		* = 17».
m = 0,1,2,...,	k ~ 0,1,2,...,

Задача (3.23) была решена при помощи программы РС2М за 2 мин. 20с., однако решение этой задачи при помощи программы РСЗМ при­ вело к увеличению времени счета до 6 минут. Это, по-видимому, объ­ ясняется тем, что линейный прогноз, заложенный в программах РС1М, РС2М, более эффективен при отыскании решения вне пограничного слоя, где, как показано выше, решение хорошо описывается линейной функцией * = 1 —t, являющейся решением соответствующего выро­ жденного уравнения. Это подтверждается и численными результатами, которые согласуются с аналитическими, задаваемые функцией х = 1 —t, до четвертого знака.

В качестве другой проблемы рассмотрим численное решение зада­ чи Коши для уравнения Ван-дер-Поля [97, 63]. Если в классическом уравнении Ван-дер-Поля у" - р( 1 - у2)у' + у = 0 перейти к новой пе­ ременной t = ж//», то уравнение примет вид еу" - (1 - у2)у' + у = 0, где е = I//*2, и задачу Коши можно сформулировать в виде

^ ~ У 2> е ^ = (1- у Ь У 2 - у и V l(0) = 2, у2(0) = - 0,66.

Здесь малый параметр е = 10-6.

Это хорошо известный пример для проверки эффективности работы вычислительных программ решения жестких систем. Для решения задачи

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

использовалась программа РС1. Точность вычислений контролировалась сравнением с «точными» результатами, полученными в [97), и не превос­ ходила 10~3. Решение задачи Ван-дер-Поля при t € [0,0,01] потребовало в десять раз большего времени счета и в двадцать раз большего числа вычислений правой части системы, чем решение той же задачи после применения A-преобразования. При этом шаг интегрирования непреобразованной задачи был в десять раз меньше.

Отметим, что программа РС1 позволила найти решение задачи Ван- дер-Поля после применения A-преобразования при t € [0,2], тогда как решение непреобразованной задачи прекращалось после t = 0,03 из-за переполнения памяти ЭВМ.

3.3. Жесткие системы

Переходя к изучению решения жестких систем обыкновенных диф­ ференциальных уравнений, следует отметить, что существует точка зре­ ния, согласно которой жесткие системы можно рассматривать как систе­ мы, получающиеся из сингулярно возмущенных уравнений при помощи некоторой неизвестной замены переменных. Следовательно, жесткие системы — это системы с многочисленными неявными параметрами, так как не существует явной группы переменных, которые можно было бы разбить на быстрые и медленные (как z и у для системы (3.16)) и нет простого уравнения для поверхности Г, как Z(y, z) = 0 для систе­ мы (3.16). Поэтому не существует асимптотической теории, аналогичной теории для син!улярно возмущенных уравнений.

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

Рассмотрим решение некоторых жестких задач с использованием А-преобразования.

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

ось х направлена по скорости полета V,	ось у	ортогональна оси х
и лежит в плоскости симметрии самолета,	а ось	z направлена вдоль

3.3. Жесткие системы		95
размаха правого крыла, имеем
dv	dO	d?0
m - = P - Q - G ^ e ,	m v - = Y - e cos*,	I , - p = M , - P , „ (JJ3)

a = 0 - 0 ,

где го — масса самолета; ур — плечо силы тяги двигателя относительно центра тяжести самолета; v — величина скорости; Р — сила тяги, приложенная вдоль оси двигателя; Y — подъемная сила, направленная ортогонально к скорости полета вверх; G — вес самолета, направленный вертикально вниз; Q — лобовое сопротивление, направленное по оси скорости набегающего потока; Мг — момент внешних сил относительно оси OZ\ Iz — момент инерции самолета относительно оси OZ; 0 — угол между вектором скорости и горизонтальной плоскостью; О — угол тангажа (угол между хордой крыла и горизонтальной плоскостью; а — угол атаки.

В общем случае силы и моменты, входящие в систему (3.33), зависят от многих параметров движения: угла атаки о , плотности воздуха р, скорости полета, угла отклонения руля высоты, угла тангажа и их производных по времени.

Основываясь на методике малых возмущений и переходя к безраз­ мерной форме, получаем систему линейных уравнений для возмущений

V - VQ	Аа = а - <*о,	„	psvn
Av = ---------,		АО = 0 - OQ,	т =
vo			2го
где s — площадь крыла. Далее знак А перед v, о,			О опускаем.

В условиях горизонтального полета для конкретных значений пара­ метров самолета система уравнений (3.33) принимает вцц:

^ = -0,104v + 0,043а - 0,10, ат

dot

— = -0,57v - 5,12а + wz, ат

(3.34)

dO

— = w z, dr

^ - 12,574v - 43,69a - 9,672a>2. ar

Корни характеристического уравнения этой системы, приведенные

в монографии [53], равны A|i2« - 7,4± 6,2t, A3 и -0,27,	А4«0,16, *2 = —1,
и разделяются по величине модуля на две группы |А|| =	|Аг1 > |Аз| > IA4I,

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

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

движением. Движение, отвечающее двум малым корням — длиннопе­ риодическим или фугоидным движением.

Таким образом, система (3.34) является жесткой на любом отрезке [О, Т], длительность которого значительно превышает длину погранич­ ного слоя (т„с < 0,1).

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

Жесткой является также система уравнений бокового движения са­ молета и полного описания движения самолета, объединяющего боковое и продольное движения.

Отметим, что подобное различие в движениях характерно не только для самолета, но и для ракет.

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

уравнений (3.34), которое следует отыскивать	в виде
v = aext, а = 6еА‘, 0 = сеА‘,	шг =	deM.	(3.35)

Здесь принято, что t = т. Подставляя (3.35) в (3.34), получаем систему линейных алгебраических уравнений относительно а, 6, с, d

'-(0,104 + А)а + 0,0436 - 0,1с = 0, -0,57а - (5,12 + Л)6 + d = 0,

(3.36)

-A c + d = 0,

. - 12,574а - 43,696 - (9,672 + A)d = 0,

Это линейная однородная система уравнений. Для того, чтобы она имела нетривиальное решение, ее определитель должен равняться нулю

	-0,104- А	0,043	- 0,1	0
	-0,57	-5,12- А	0	1	=	0.
	0	0	-А	1
	-12,574	-43,69	0	-9,672

Раскрывая этот определитель, приходим к характеристическому уравнению

А4 + 14.896А3 + 94.774А2 + 9,215А - 3,948 = 0,

корни которого, уточненные при помощи метода Ньютона, имеют зна­ чения

А! =0,1596, Аг = -0,2651, А3,4 =-7,3952 ± *• 6,211.

(3.37)

3.3. Жесткие системы

97

Учитывая систему уравнений (3.36) и принимая в качестве свобод­ ного параметра с — 1, получаем, что собственным значениям (3.37) будут соответствовать следующие значения параметров a, b, с, d

At :

а =-0,36795,

6 = 0,06996,

с =

1,

= 0,1596,

Л2 :

а = 0,6563, Ь= -0,1316,

с = 1,

d = -0,265,

А3 =

-7,3952 -

* • 6,211 : а = 0,57385 • 10~2 -

i •0,5971 • 10_3,

6 =

1,2664-t -0,7274,

с = 1 ,

d = —7,3952 —* • 6,211.

Поэтому общее решение системы (3.34) будет иметь вид:

/

v

(

-0,36795

/

0,6563

а

= С\

0,06996

е°-159б< + С2

-0,1316

e -0,265< +

0

1

1

\ш г

V

0,1596

\

-0,265

/

( 0,57385 - 10_2\

/

-0,5971 • 10-3 \

+

cos 6,211 - Сз

Ч 664

+ с 4

-0,7274

I

0

V

^

-7,395

)

1

- 6,2! 1

)

( 0,5971 • 10_Л

( 0,5738510“Л

\

sin 6,211 • Сз

0,7274

+ с 4

1,2664

0

1

6,211

J

-7,395

У

/

где С |, Ci, С3, С4 — произвольные константы, определяемые из на­ чальных условий. Если взять начальные условия в виде

1 = 0, « = 0 =

= 0, а = 1,

(3.38)

то

Ci =-0,2969, С2 = -0,17106, С3 = 0,46798, С4 =-0,5576.

/

Задача (3.34), (3.38) интегрировалась на отрезке 1 € [0, 5] при по­ мощи программы РС1. Причем при той же точности вычислений время счета задачи после применения Л-преобразования оказалось почти в два раза меньше.

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

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

'	^ =	77,27(02 + 0,(1 -	8,375 • 1<Гб01 -	у2)),
	dyо	1
*	"dT =	77^7^ 3 “ ^ + У	Х	Р-39)

= 0, Ш ( у 1 - у з ) .

Эта система уравнений описывает знаменитую химическую реакцию с предельным циклом в трехмерном случае, так называемый «орегонатор». Это реакция между НВЮг, Вт- и Ce(IV). В [92] система уравнений (3.39) интегрировалась при начальных условиях

2/1(0) = 3, 02(О) = 1 , 0з(О) = 2.

(3.40)

При этом решение задачи обычным методом прогноза-коррекции, в котором стартовой процедурой был метод Рунге—Кутта, не получилось из-за высокой жесткости задачи. Решение удалось получить методом, приведенным в [93].

В монографии [63] приводятся графики решения задачи (3.39), (3.40) (см. рис. 3.3, 3.4) и отмечается, что это пример жесткой си­ стемы уравнений, решение которой быстро изменяется по величине на много порядков. Поэтому данный пример может служить серьезным испытанием для программы численного интегрирования.

Попытка решить задачу (3.39),(3.40) при помощи программы РС1 не привела к успеху, однако после применения к этой задаче Л-преобразования решение было получено на интервале, равном одному циклу орегонатора, т.е. для t € [0, 303]. Но при этом пользовалось очень большое время счета на PC АТ 286/287 (около одного часа). После при­ менения к задаче (3.39), (3.40) Л-преобразования, последняя примет вид:

3.3. Жесткие системы

99

dyi _ Р\_	dy2 _	Рг dy$ _	Pi		dt _	1
~ d \~ Q '	~ d \~	Q ' ~ d \~		Q '	dX ~ Q'		(3.41)
У \ (0) = 3,	У2(0) = 1, !/з(0) =			2,	<(0)=	0,

где Pi — правые части системы уравнений (3.39), Q = \J\+ P } Л-Р^Л-Р%■

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

Было предпринято несколько попыток уменьшить время вычисле­

ний.

Так, ни к чему не привела попытка строить итерационный процесс на ( т + 1)-м шаге согласно формуле

: У,-,т+1 -	[&> + hfiiVm+l’ *m+l) -	y,',m+l] >
в которой параметры	подбирались, исходя из	некоторых оценок

нормы матрицы Якоби правой части системы, обеспечивающей условия сходимости.

Зато существенное улучшение ситуации наступило после примене­ ния итерационной схемы, предложенной в [55], в которой внешние итерации осуществляются по Ньютону, а внутренние — по Зейделю. Применительно к исследуемому случаю с учетом того, что, в виду авто­ номности системы (3.39), последнее уравнение системы (3.41) является следствием интегрирования первых трех уравнений, рассматриваемая

итерационная схема примет вид
(Ь+1)	(к)	У|,м+Ь/|(у»+|А»'и)-у1д+|	’
Vl,m+1 ~ 2^1,т+1 “		**„_[
(fc+l)	(к)		^
У2,т+1 ~ У2,т+1 ~
(fc+l)	(к)	y3.m+ft/3(y,nl|i*mtl)-y]J|.).|~(i>ltm+)|~y!.j,+|)^3l+(y2.M+l|~!>2.1-n)ft^2
Уз,т+	1 ~ Уз,т+ 1 ~		- 1	Лбзз

Данная схема определяет последовательные итерации по Ньютону и только одну итерацию по Зейделю, сохраняя структуру метода Зейделя. Здесь bij — компоненты матрицы Якоби преобразованной системы уравнений (3.41), которые вычисляются по формулам

fij	£	1к ~ /. £ fkfkj
	fc=l	к= I	(3.42)
Здесь /,■ — правые части системы уравнений (3.39), причем пред­
полагается, что /„ = I, fij	—	производная	по переменной у,,

г, j = I, и —1, п = 4.

Но еще более эффективным вычислительный процесс становится, если итерации проводить в соответствии с формулой (3.24), т.е. решать задачу при помощи программы РС1М. В этом случае один цикл пове­ дения орегонатора, соответствующий изменению времени t в пределах от 0 до 303, просчитывается на IBM PC AT 286/287 за 90с. при заданной