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

цы Якоби / , которая не меньше максимальной обусловленности матриц

/* . От этой обусловленности зависит погрешность, которая накапливает­

ся в процессе построения ортонормированного баэиса Ух,

, Ут (1.1.18).

Во-вторых, обусловленность метода ортогонапиэаиии зависит от выбора

вектора	и определяет погрешность операции (1.1.23) выделения нз
вектора (} орта искомого решения с]Х/с! X.
В книгах [63, 35] рекомендуется в качестве О выбирать вектор-
строку	[0,. . . , 0,2]. Однако такая рекомендация вызвана использованием

метода ортогоналиэации для решения системы уравнений классического вида (АХ=В ) . Вообще же говоря, обусловленность операции (1.1.23) тем больше (а погрешность тем меньше), чем ближе вектор (? к искомому решению Ут+1. Действительно, если совпадает с Ут+,, то все скаляр­ ные произведения (О, У/), / = 1,. ., т, равны нулю, и Фт+х^т+х) ~ 1- Бели же, наоборот, вектор О ортогонален к Нт+1, то он принадлежит подпространству Р„ ,, т.с. [/,„ +, = 0, и операция нормировки в (1.1.23) становится невозможней хотя бы потому, что она приводит к делению

на нуль. С этой точки зрения использование операции о г ((/, @) в про­ цессе продолжения решения по параметру создает выгодную ситуацию для выбора вектора О, близкого к искомому Ут*-\. Если мы переходим из точки Х&_1 кривой множества решений К (рис. 1.2) в точку X*, то ра­

ционально в точке ХА. б качестве О выбрать вектор ^Х /Л Х !^ ^ ( в точ­

ке _ т), который будет тем ближе к искомому вектору <1Х/4Х\Х/с, чем ближе ХАк ХА_ , . Это обстоятельство может быть использовано в качестве критерия для оценки величины выбранного шага по X при движении вдоль кривой К. Так, малая величина нормы ({/т+1, {/т+1) 1^2 вектора Ит+\

свидетельствует о том, что вектор О = <2Л7<2Х|х^ ^ сильно отличается

шага от ХА-1 к ХА = Х*_, + ДХ сильно изогнулась и, чтобы избежать на­ копления погрешности при интегрировании задачи Коши по X, необходимо уменьшить шаг ДХ.

Рассмотрим теперь конкретные алгоритмы построения решений сис­ темы уравнений (1.1.3) путем численного интегрирования задачи Коши

по параметру X (1.1.24), (1 Л.,25) :

ах

« ( / . а ) . *(о)=дг(0).

~а\

Из всех известных алгоритмов рассмотрим и будем исследовать в дальней­ шем только метод Эйлера, модифицированный метод Эйлера и метод Рунге — Кутта 4-го порядка. Для удобства введем следующие обоз­ начения:

*(Х*) = * С*>,		( 1 .1 .2 9 )
ах/ах-х,	ах/ацКк= * ( х Л) = * (Л ),	( 1 .1 .3 0 )
7 д а Л) ) = 7 ( ^ (А )) = 7 (Л ),		( 1 .1 .3 1 )
( К^к)~0{к)•		( 1 .1 .3 2 )
Простейшим	методом интегрирования задачи Коши	(1.1.28) является

метод Эйлера. Он основан на замене движения по кривой К движением

по касательному вектору	аХ/	аХ = дг. Переход по методу Эйлера от X*,
где известно Х(Хк) = X	, к	Х*+1 = X* + ДХ сводится к	следующим
вычислениям: вычисляется расширенная матрица Якоби 3 (к)			(х (к))‘>

она дополняется снизу до квадратной вектором-строкой ^(к)> в качестве которого удобно взять вектор Х ( * _ 1) = а Х / а Х ^, полученный на

предыдущем шаге; с помощью операции ог1 строится орт касательной

ккривой К множества решений

ах/ах\ч = Х ( * ) = о г1( / л , * ( * _ ! > )

и, наконец, получают приближенное решение АгСл+1) = ЛГ(Х*+1), прод­

вигаясь на ДХ вдоль вектора касательной, как X ^к+1 ) = х (к) +Х(Л) ДХ. В таком алгоритме остается произвол в выборе вектора б на первом шаге, при переходе от Х(0) = О к Х (() . Так как обычно начальная точка X = О,

где известно решение -^(о)> не является предельной, то на первом шаге можно использовать в качестве (2(0) вектор [0,. . . , 0 , 1]. Окончательно алгоритм метода Эйлера для интегрирования задачи Коши (1.1.28) прини­ мает следующий вид:

Хо=0, Х , о)=Х(0), 0(0) * [0,0, - - -, 1];

(1.1.33)

Х*.ц = X* + ДХ*, Х(*) " ОП(/(*), О.(к)),

Х {к+ 1 ) = Х (к) + *(*)

<2(*+1) =

Л = 0 ,1 ,2 ,...

С	помощью этого алгоритма	кривая
К	множества решений уравнения	(1.1.3)

приближенно представляется в виде ло­ манной Эйлера Ь (рис. 1.3).

Как известно (см., например, [35]), на каждом шаге такого алгоритма накап­ ливается ошибка порядка О (АХ3 ).

Более точное приближение множест­ ва решений К (с погрешностью на каж­

дом	шаге порядка О (АХ3) обеспечи­
вает	модифицированный метод Эйлера

(см., например, [35] ). Его алгоритм в нашем случае принимает вид:

Х „ = 0 , * ( 0 ) = * ( 0 ) , ( ? ( о ) = [ 0 , 0 . .	. , 1 ] ;
^ Л + 1 = ХА + Д Х Л , • * ( * ) = о г ! ( У (Х(к)), ( 2 , * ) ) ;
(*н)	(1.1.34)
*(*кМ) = ОГН-МЛ<** + 1),'@<А <1)^’	О =

Х (к + 1) " Х (к) + “ <Л'(к) +Л'(Л + I ) >^ »к

к =0 ,1 ,2 ,.

Еще лучшую точность (погрешность порядка О(АХ4) дают формулы 3-го порядка метода Рунге - Кутта:

Ао=0, * (0)=*(0), 0(о, = [0.......

0,11;

*</*) = ог1(7(^(Л)) , е (Л)),

= 0 ,1 ,2 ,...

И, наконец, погрешность 0(ДХ5) обеспечивают формулы 4-го порядка методаРунге —Кутта, наиболее употребительные среди которых принимают

зэ

в нашем случае следующую форму:

Х0 =0, Г (0)=ДО), <2(0) = [0.........0 , 1];

Х*+1 = Х* + ДХЛ, X ^ = оП(/ (-^(Л)), 0(к))>

*(к) ~ ° « (Ф «)+ ~2*<*> Л^ )’ Х(к))'
		+\ *№ д4	$!)■	< '1-36>
*й,	=	♦4 3>>д**ь	4 3))'
*(* ♦ 1) = •*(») +		<*№ +	+ 2 *(к)	+ Х{к) ^ •
О » .,) =*<?>> ■		* - 0 , 1 , 2 , . . .

Все эти алгоритмы относятся к алгоритмам явного вида и для них характерно накопление погрешности в процессе непрерывного продолже­ ния решения. Возможно применение и других алгоритмов явного вида, как, например, алгоритмы метода Адамса - Штермера и т.п. (см., напри­ мер, [35]).

12. Обобщенные формы дискретного продолжения решения

Напомним	сначала метод Ньютона - Рафсона для решения уравнения
с одним неизвестным
П Х ) =0.	(1.2.1)
Идея метода	состоит в том, что некоторое приближенное решение X®

этого уравнения уточняется с помощью замены в (1.2.1) функции Р(Х) ее линейной аппроксимацией. Для этого Р(Х) раскладывают в ряд Тейлора

в окрестности	и	сохраняют только линейные слагаемые. Геометриче­
ская интерпретация		такой замены	показана на рис. 1.4. Тогда уравнение
(1.2.1) заменяется на следующее:
р(*<'>) + р'(хЮ)(х- х(0)= о,			'=ар/ах.	(1.2.2)
Решение этого уравнения
Х = Х< 1+х) =*<'> -		(РХхЮу-'Р^хМ)		(1.2.3)

и дает новое приближение для решения уравнения (1.2.1). Повторяя эту операцию, мы строим итерационный процесс, который при удачном началь­ ном приближении позволяет найти решение уравнения (1.2.1) с какой угодно точностью (рис. 1.4). При этом, конечно, в процессе итераций Р' = йР/йХ не должна обращаться в нуль. Иногда удобно записать уравне­ ние (1.2.2) как уравнение для приращения

Р(Х&) +Г'(Х(Г>) Д *</+ = 0.

(1.2.4)

Получив отсюда ДХ(,+ *) = - Р { Х ^ )/Р '(Х ^) , найдем

Если У11 т-мерная вектор-функция, а X - вектор в линейном нормирован­ ном пространстве Кт , то метод Ньютона - Рафсона обобщается и требует на каждом шаге итерационного процесса решения следующей системы

линейных уравнений для вектора приращения ДХ*-'* ,) :

/ ( Х (П) Д Х (,*П =	- Г ( Х (П).	О-2-6)
здесь 1(Х) = дР/дХ -	матрица Якоби вектор-функции Р.	Эта матрица
должна быть не вырождена.

При исследовании решений системы нелинейных уравнений, содержащих

параметр Р,
Р(Х,Р) =0,	(1.2.7)

обычно возникает задача определения влияния изменений параметра задачи Р на эти решения. Как отмечалось во Введении (В.1), при выполнении условий теоремы о неявных функциях решения-К уравнений (1.2.7) явля­ ются непрерывными и дифференцируемыми функциями Р:

Х = Х(Р).

(1.2.8)

Простейшей формой представления этой зависимости было бы отыска­ ние решений Х(к) =Х(Рк) для некоторого множества значений параметра Р0 < Л < . . . <Рк< . . . < Р„. М. Лаэй [447, 448] показал, как при исполь­ зовании метода Ньютона —Рафсона можно экономично организовать процесс продолжения решения по параметру. Его предложение, по сущест­ ву, свелось к использованию решения для предыдущего значения параметра задачи Р в качестве начального приближения для текущего значения пара­ метра. Этот процесс можно записать в виде

	= % - • > ;
/(■ $)• р^ л х % °		- - ■ ''0 $ . '*>■	<>•*■«>
у(/+1) = у(/) +ДУ</+1)
А (Л)	Л(*:)Т а Л (*) *
/ = 0 ,1 ,2 ,... при		IIДЛГ^лу > II> е;	* = 1 ,2 .........

здесь использованы следующие обозначения:			= Х(Рк), Х ^ /-е прибли­
жение ДГ(#) в итерационном процессе Ньютона -			Рафсона (1.2.9), который
продолжается, пока	норма	вектора приращения		Д Л Т ^1* превышает
заданную точность	е > 0.
На рис. 1.5 показана геометрия процесса (1.2.9)				для уравнения с одним
неизвестным при переходе от		1 к Рк. Искомым		множеством решений

уравнения Р(Х, Р) = 0 является кривая К, но которой поверхность Р(Х,Р) пересекается с плоскостью X, Р. Итерационный процесс Ньютона - Рафсона проходит в плоскости Р = Рк с начальным приближением = Х ^ _ х у

Этот же рисунок наглядно показывает, что трудности, возникающие в окрестности предельной точки Т, связаны с переходом от Рк к Рк -ц, кото­ рый выводит процесс (1.2.9) из области, где существует решение. Как видно, эти трудности обусловлены тем, что решение разыскивается в плоскости Р=Рк+1, которая не имеет пересечения с кривой множества реше­ ний К. Таких трудностей удастся избежать, если для каждого К организо­

вать итерационный процесс Ньютона —Рафсона поиска	X ^	в плоскости
М(*), которая ортогональна кривой К при *-*<*)■ Но	пока	нс найдено

Х(к), плоскость М(*) неизвестна. Однако можно искать решение в плоско­ сти М (*) , близкой к М(*). Рассмотрим один из способов задания М(Л). Так же, как и в предыдущем параграфе, будем считать неизвестные и параметр задачи равноправными и введем вектор

X = [Хи . . . , Хт, Хт + ,= Р] Е	+ |.
Тогда систему уравнений (1.2.7) можно записать в векторной форме
^ (Г ) = 0.		(1.2.10)
Здесь Р = [Т7! , . . . , Рт ]т - вектор-функция.
Пусть I —величина шага, с которым мы стремимся		двигаться вдоль К.
Тогда к плоскости М(к) будет близка плоскость		проходящая через

з*

точку !) + (• с/ЛГ(*_1)/с/Х<ЕК2 гак, что она ортогональна касательно­

му к кривой К в предыдущей точке Х (к„ ^единичному векторуёХ^к_ х^ ' К (рис. 1.6). Плоскость М(к) поэтому определится уравнением

( ^ 7 Г '

>+, 2 р ) ) - 0 .

(1.2.11)

которое требует равенства нулю скалярного произведения (т.е. ортогональ­ ности) принадлежащего плоскости М(&) вектора

Д* (	>	)	(1. 2. 12)
и вектора6 Х (*_ 1) А*Х, касательного к кривой К в точке X (* _ ,).
Таким	образом, бпределение	следующего после X (*_1)	решения

X сводится к нахождению решения уравнения (1.2.10) в плоскости

М*ку, т.е. к совместному решению уравнений (1.2.10), (1.2.11)

Г(Х) = 0,

(1.2.13)

Здесь второе уравнение получено из (1.2.11) с учетом того, что вектор

/сГХ является единичным, т.е.

(<**(»-1)А*Х, <^(*-1)7</Х) - 1.

(1.2.14)

Итерационный процесс метода Ньютона - Рафсона для системы (1.2.13) имеет вид

-	(О)	-	-	(1.2.15)
X (к) -
—	/«	—и+и	—м \	(1.2.16)
У	( А ^	) АЛ- (*)	= - У (Л Г (* ))’

(4	=“-		~л
- 0 + 0	-(»)	- 0 + 0	(1.2.17)
X (*)	= X (*) +	А Х 00 ’
	при	1> е;	к = 1,2,
Здесь, как и ранее в		§1.1, введена расширенная матрица Якоби вектор-
функции Р
1 = дР/дХ.			(1.2.18)

Для уравнения с одним неизвестным Х х = Х и параметром Хг = Р гео­ метрия такого итерационного процесса показана на рис. 1.6.

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

нения (1.2.17) представляет собой длину проекции вектора (Х^ку

— Х^к_ 1 у) на направление орта йХ^к_ 1уД/Х. Начала этих векторов сов­

падают и лежат в точке Х^ к_ ху, а конец вектора (.Х^ку - Х (^к_ 1 у) ле­

жит в плоскости М*. Но так как плоскость Мк ортогональна вектору

ДА, то в процессе итераций длина этой проекции остается постоян­ ной и такой же, как на начальном шаге (1.2.15), т.е. равной г. Поэтому правая часть в уравнении (1.2.17) обращается в нуль, а само это уравнение приобретает более простую форму:

’ Л^ {(^)1)) = 0‘

(1'2'19)

Геометрически это уравнение означает требование ортогональности

поправочного вектора АХ*1^ * к касательному орту <1 Х ( к_ 1 у/(1 \

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

то по аналогии с (1.2.19) условие (1.2.17) можно заменить следующим

(О	-(«+1)) = 0.		(1 2 21)
(*)’	ТАЛ }	'

Геометрия итерационного процесса с таким условием показана на рис. 1.7. в плоскости К2. В этом процессе на каждой итерации корректируется по­

ложение плоскости М(к) >в которой ищется решение.

Если к тому же вектор 5 на каждой итерации нормировать так, чтобы он имел длину г, то получится итерационный процесс, проиллюст­ рированный в плоскости К2 рис. 1.8. Этот процесс обеспечивает величину шага I по хорде кривой К с заранее заданной точностью. Его алгоритм имеет вид

					"		Г	+
								(1.2.22)
/2(0	“	(/+1)	) = 0,					(1.2.23)
^ (* )’		(к)
7(^+1)	= *еГ(0		+ д у		(к)	|/	II	ч
5 (А)	'^(А )		+ ДЛ			>'	" Е(А)	Л Л (*>
-(/+1)	_	у		4Т0+1)
(А)		Л (А-1 >		*	(к)
			при	II Д * (/(У)) II >				е;

Заметим, что все рассмотренные здесь алгоритмы допускают обыч­ ную для метода Ньютона - Рафсона модификацию с заменой матрицы

1{Х ^ ) , которая меняется от итерации к итерации, на матрицу / (ЛГ(° ^ )

первого приближения.

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

уравнений вида (1.2.10) с т + 1 неизвестными компонентами вектора X совместно с некоторым дополнительным условием.

Так, в процессе М. Лазя (1.2.9) используется простейшее дополнитель­ ное условие Х т + 1 =Рк, ав процессе (1.2.15) —(1.2.17) -условие (1.2.19) и т.д. Все эти дополнительные условия определяют в К ,„+1 некоторую плоскость М(#), которая может изменяться от итерации к итерации. Обобщая этот подход, можно попытаться сформулировать дополнитель­ ное условие, определяющее некоторую поверхность, и искать решение уравнения (1.2.10) как точку пересечения множества решений К с этой поверхностью. Так, если в качестве такой поверхности выбрать сферу ра­

диуса I с центром в точке						то придем к			совместному реше­
нию следующих уравнений:
Г(Х) =0,									(1.2.24)
(X - Х {к_ ,), X -			Х {к_! >) - Г2 - 0.						(1.2.25)
Алгоритм метода			Ньютона		— Рафсона		для	такой	системы уравнений
имеет вид
- ' Ч	Х ^ / а к х Щ				• * < » - „		+ ?,*>:
7 ( Х ^ ) й Х < ‘;к\> —				я * * » ) .					(1.2.26)
			=			О	) .		(1.2.27)
у<1+1>	у (о	+	а х (,+1)		ьО+1) -		у(/+ 1)
А (*)	л (*)	+		(к)	■ *	(*) '	Л (к)	‘
		при		II А Х ^ У		II > е ;

Истолкование рассмотренных выше итерационных процессов как про­ цессов совместного решения основной системы уравнений с дополнитель­ ным уравнением позволяет рассматривать их с общей точки зрения на метод Ньютона - Рафсона, которая подробно развивается во многих мо­ нографиях ([366, 35, 481,212] и др.). В них детально обсуждены вопросы сходимости метода. Мы только отметим, что для сходимости итерацион­ ного процесса метода Ньютона - Рафсона начальное приближение обычно не должно слишком сильно отличаться от искомого решения. В построен­ ных выше итерационных алгоритмах по самому смыслу метода продол­ жения решения это требование удовлетворяется при достаточно малых величинах шага I по параметру продолжения X.