Дифференциальный метод продолжения решений систем конечных нелинейных уравнений, зависящих от параметра Текст научной статьи по специальности «Физика»

УЧЕНЫ Е 3 А ПИ С К И Ц А Г И Том XVII 1986 №5

УДК 629.735.33.015 519.6

ДИФФЕРЕНЦИАЛЬНЫЙ МЕТОД ПРОДОЛЖЕНИЯ РЕШЕНИЙ СИСТЕМ КОНЕЧНЫХ НЕЛИНЕЙНЫХ УРАВНЕНИЙ, ЗАВИСЯЩИХ ОТ ПАРАМЕТРА

М. Г. Гоман

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

Изучение стационарных (установившихся) форм движения во многих прикладных задачах приводит к необходимости расчета и исследования семейства решений систем нелинейных уравнений, зависящих от ряда параметров, и проведения бифуркационного анализа этих решений [1—3]. Хотя существует достаточно большое количество различных методов решения систем нелинейных уравнений, наиболее приемлемым для решения подобных задач является метод продолжения по параметру [3—6].

Рассмотрим векторное уравнение

F(x, Х) = 0, (1)

где F — вектор-функция размерности п, х — вектор фазового состояния размерности п, X-—параметр, без ограничения общности принимаемый скаляром. Функции Fh i — 1, 2, ..., и непрерывны и имеют непрерывные частные производные по всем переменным xt и параметру X. Под вектор-функцией F в общем случае можно понимать некоторое гладкое отображение.

При непрерывном изменении параметра X решение х(к), задаваемое неявно векторным уравнением (1), также будет изменяться

непрерывно. Связь между малым изменением параметра и соответствующим смещением точки решения ёх определяется соотношением

дР йх+-^-(И = 0, (2)

дх

д\

получаемым приравниванием нулю дифференциала левой части (1).

Если матрица Якоби левой части не вырождена, то траектория решения л: (X) может быть получена интегрированием дифференциальной системы уравнений [4]

с1х

ёХ

дР\-1 дх

дР_

д\

при начальных условиях л:(Х0) = х0, где Т7 (х0, Х0) = О;

дР

матрица Якоби, =

д1

дх

дР2

дР„

дх

(3)

т-

Вырождение матрицы

дХ ’ '

в точках траектории решения х(Х)

совпадает с разворотом траектории или появлением новых ее ветвей и делает невозможным использование уравнений (3) [7].

Однопараметрическое семейство решений х(Х) в расширенном пространстве г = (х, I) образует траекторию (рис. 1), которая может быть естественно параметризована ее длиной [5]. Такой подход позволяет устранить искусственно возникающие особенности в предельных точках при выборе в качестве независимого аргумента одной нз фазовых переменных или параметра.

Рис. 1

/ X

1-предельные точхи О

2-точки 5ет5ления йеЬ Ф =0

При движении вдоль траектории решения (х, X) приращения компонент вектора состояния и параметра будут определяться проекциями касательного к траектории единичного вектора 5. Этот касательный вектор, определяемый бесконечно малым смещением (с1х, <А) в пространстве г, перпендикулярен строкам мат-

/ дР дР \ дР , „ п. ,

рицы ------, — =--------, являющимся градиентами функции г Ах, л),

\ дх дХ ) дг

г=1, 2, 3, ..., /г, задающими в пространстве г гиперповерхности единичной коразмерности, пересечение которых и образует искомую траекторию решения.

др•

Если среди векторов г = 1, 2, 3, ..., п нет линейно зави-

/дР\

симых, т. е. если ранг —— = и, то в совокупности с единичным

касательным вектором 5 они образуют в пространстве г базис, а матрица Ф размерности (я+1, п 1), составленная из них, не вырождена

Ф =

дРг ал ал \

дх1 дх2 дхп ах \ / дР\

дР„ дРп дРп 1! Ъ ^1

дхі дх2 дхп ах /

5, 52 . • 5я+1 /

(4)

Компоненты вектора 5 могут быть получены посредством ре-

др

шения линейной системы уравнении -^—5 = 0 (2) при фиксации

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

др

Поскольку вектор 5 ортогонален строкам матрицы , его компоненты определяются соотношениями

пп + \,і <ІЄІ Ф

(5)

где Ап+^1 — алгебраические дополнения элементов нижней строки

Гп+1

матрицы (4), а <1е1Ф=1/ £ Ап+ц.

/=1

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

йг

:о5

или

(ІХі

<и

<и

— і— 2, 3, ..., п,

= а8п+1 = а йеі(~) І йеіФ

(6)

при начальных условиях г — г0, Р(г0) = 0; здесь г1—„время", скалярная величина, параметризующая траекторию, а о—параметр, знак которого определяет направление движения вдоль траектории.

Уравнения (6) позволяют получить непрерывную ветвь решения (1) в пространстве г. При этом, как видно из последнего уравнения в (6), в предельных по параметру точках траектории решения (см. рис. 1) происходит смена знака определителя матрицы

Якоби , поскольку величина сІеіФ положительна. Эти точки

соответствуют точкам бифуркаций динамической системы л: — Г(х, X), в которых, например, происходит слияние устойчивых и седло-вых особых точек.

Численное интегрирование уравнений (6) будет приводить к накоплению ошибки и отклонению от истинной траектории из-за отсутствия коррекции соотношения (1).

Следуя [6], можно построить непрерывный алгоритм, обеспечивающий устойчивость траектории решения (1), задаваемой уравнением (6). Дополним уравнения (6) линейным векторным уравнением относительно Т7:

-^+*/=■ = 0, (7)

аг

решение которого имеет вид Р — И(г^е~м, где положительная константа £ определяет скорость сходимости невязки =

= ~\/ уравнений (1) к нулю при увеличении

1=1

Раскрывая полную производную в левой части (7) и переходя к переменной г, получим уравнение

дР йх , дР йк дР е1г ис /0\

--------!---------=----------= — кР, (8)

дх еН дк еН дг еИ

которое после скалярного умножения его на вектор 5 может быть объединено с уравнением (6) в единую систему уравнений

Ф-Ё*. = '~кР

dt

о

) • (9)

В предельных по параметру точках траектории (см. рис. 1),

дР

где вырождается матрица Якоби > матрица Ф остается невырожденной, и соответственно уравнения (9) не имеют особенностей. В силу этого система (9) может быть преобразована к виду

(10)

Система (10) при о = +1 обеспечивает не только движение вдоль траектории, но и асимптотическую сходимость к ней из некоторой достаточно широкой области значений z, не принадлежащих траектории решения. При а = 0 реализуется только сходимость к траектории решения z(t) по семейству траекторий, ортогональных семейству изолиний Р(х, Х) = С, где С — постоянный вектор размерности п.

Таким образом, получена единая система дифференциальных уравнений (10), которая обеспечивает расчет искомого решения с автоматическим уменьшением невязки flFfl^O, возникающей при интегрировании этой системы с конечным шагом At, а также при выборе начальной точки z0, не являющейся решением (1).

Наличие на траектории решения z{t) критической точки, где происходит ветвление решений (см. рис. 1), обусловлено вырождением матрицы Ф из-за появления линейно-зависимых строк dF

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

Если в точке на траектории решения z(t) ранг то точка

неособая, существует единственное касательное направление S (5).

Если ранг <п, то в критической точке пересекаются две

или более ветвей, что определяется нелинейными членами разложения вектор-функции F(x, Ц в окрестности особой точки.

Для прохождения точки ветвления при интегрировании уравнений (10) может быть применен метод, основанный на использовании представления о структуре решения уравнений в ее окрестности. В случае пересечения четного числа, в частности, двух ветвей, переход через точку ветвления изменяет направление движения на противоположное из-за смены знака алгебраических дополнений (5), так что при конечном шаге интегрирования может возникнуть зацикливание (точки / и 2 на рис. 2, а). Поэтому необходимо осуществить смену знака параметра а в уравнениях (10). В случае пересечения нечетного числа ветвей (рис. 2, б) направление движения при переходе через критическую точку сохраняется. Штриховые линии на рис. 2, а, б показывают траектории приближения к решению в окрестности критических точек при о = 0-

Для определения направления движения по второй ветви в двукратной точке ветвления (рис. 2, а) может быть использована информация об ориентации гиперплоскости, касательной к ее ветвям в точке ветвления. В окрестности траектории решения z(i) скалярное произведение FT-F в первом приближении выражается следующей квадратичной формой FT -F х AzT ЛДг, где

\ dz ) \ дг

Дг — приращения вектора z. Поскольку вдоль траектории решения квадратичная форма должна равняться нулю, то матрица А вырождена и собственный вектор, соответствующий нулевому собственному числу,совпадает с направлением касательного вектора S. В двукратной точке ветвления возникают два нулевых собственных числа, и соответствующая им пара собственных векторов определяет искомую гиперплоскость, касательную к двум ветвям решения z(t). Одномерный поиск направления в этой гиперплоскости по смене знака величин F, позволяет определить направление второй ветви [3].

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

■квадратная матрица размерности (/г+1, п + 1),

-траектория решения Sточке Seml/!enxa приближения к решению

Рис. 2

Рис. 3

Р, задаваемой, в частности, в виде /% = Т7 + £ (г — г*), з — малый параметр. При малой деформации векторного поля пересекающиеся в критической точке ветви распадаются на несколько неособых траекторий [7]. На рис. 3 приведена схема обхода точки ветвления. При подходе к ней введение деформации е=^0 приводит к отклонению от исходной ветви (е = 0) и движению по неособой траектории, гладко аппроксимирующей пересекающиеся ветви. Знак параметра е определяет положение этой неособой траектории. После снятия деформации (е 0) происходит возврат на новую ветвь исходного решения (е = 0).

Уравнения (10) при о=+1 и Р = 0 обеспечивают движение вдоль траектории, при а = 0, ^^0 имеет место сходимость к траектории решения. Рассмотрим вычисление поправки, определяемой

системой (10) при о = 0 и следующем выборе коэффициента

где М — приращение параметра в линейном приближении. Она определяется следующей системой линейных уравнений

■£л* — />; (11)

Дгт 5 = 0. (12)

Видно, что однозначный выбор корректирующей поправки Дг, определяемой уравнением (11), осуществляется с помощью условия об ортогональности поправки Дг к касательному вектору 5 [уравнение (12)].

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

1) запуск процедуры, вычисление начальной точки х0, Х0: Р{х0, Х0) = 0 с помощью известных методов при фиксированном А = Х0 или с помощью уравнений (10) с варьируемым параметром при <з = 0 или по соотношениям (11), (12) до удовлетворения требуемой точности по невязке ||/^Ц. Варьирование параметра в ряде случаев значительно расширяет область начальных значений, при которых реализуется сходимость к решению;

2) вычисление единичного касательного вектора 5 по формулам (5);

3) определение конечного шага М вдоль траектории из условия, что приращения по всем координатам Дг/ не превышают

некоторых заданных значений Дг*: Д^ = ш1п(Д2^/5*,);

к

4) вычисление следующей точки на траектории по соотношениям (6) и осуществление коррекции с помощью уравнений (11), (12) до требуемой точности. Далее операции, начиная со второй, повторяются.

Уравнения, подобные (9), могут быть использованы для расчета отдельных ветвей корневого годографа р(к) характеристического уравнения, задаваемого комплексным уравнением в(р, Х) = 0, где (/ — комплексная функция, р — комплексное число, X — действительный параметр. Уравнения для изменения корня р и параметра X при условии ортогональности корректирующей поправки к корневой траектории, в отличие от уравнений, приведенных в работе [8], будут иметь следующий вид:

л' л1 I < 1

где = Ох = -д^-; о = + 1; д «■ положительная константа.

Интегрирование корневой траектории с начальным значением в одном из нулей характеристического уравнения в(р, Х) = 0 позволяет получить корневой годограф с постоянным шагом изменения по величине р, что очень существенно при графическом построении годографа. Для прохождения критических точек ветвления корневого годографа могут быть использованы описанные выше приемы.

На базе метода непрерывного продолжения решений по параметру может быть построен алгоритм систематического поиска множества решений системы нелинейных уравнений /'(л;) = 0, где /% х — «-мерные векторы. Подобная задача играет важную роль, в частности, при анализе нелинейных динамических систем.

Рассмотрим подсистему Р^)(х) = 0, где Р(к){х) — вектор-функция размерности п— 1, получаемая из вектор-функции Р (х) отбрасыванием ее &-й компоненты. Расчет траектории решения системы Г(к) (х) = 0, начинаемый из какого-либо известного решения системы Р(х) = 0, позволяет по смене знака функции Рк(х) вдоль рассчитываемой траектории вычислять искомые решения системы Р(х)=0. При этом, если подсистема />)(*) = 0 имеет единственную ветвь решения, то удается найти все действительные решения системы Р(х) — 0, поскольку они должны лежать на этой ветви. Использование достаточных условий для доказательства существования единственной траектории решения у подсистемы Р^х) = О, получаемых на основе теоремы о неявной функции [9], существенно сужает класс подобных систем, поскольку в общем случае траектория решения может содержать [предельные точки и точки ветвления. Поэтому продолжение решения по траекториям, задаваемым последовательно подсистемами Р(/!)(х)=0 для 6 = 1, 2, ..., п, с прохождением особых точек на траекториях значительно расширяет эффективность подобного алгоритма.

Продолжение решения по параметру является существенным элементом гомотопического метода решения систем нелинейных уравнений, суть которого заключается в искусственном введении параметра и замене исходной нелинейной системы новой с известным решением при начальном значении параметра. При непрерывном изменении параметра сама система и ее решение в ряде случаев могут быть приведены к искомым [10]. Примерами гомотопической деформации могут быть следующие замены: й(х, Х) = = Р(х) + (\—\)Р(х0); 0(х, X) = \Р(х)-\- (1 — X) (х — х0), где [0, 1].

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

нениях органов управления нескольких устойчивых установившихся режимов полета [2].

В предположении о линейной зависимости аэродинамических сил и моментов от параметров движения установившиеся режимы полета определяются решениями нелинейной системы уравнений следующей структуры: Р{х, и) = Лх + /(х) + Ви, где х = (а, р, ш^., <«у)«дт—фазовый вектор, а, (3 — углы атаки и скольжения, шд., «>г— проекции вектора угловой скорости на связанные с самолетом оси координат, и = (ЬВ) 8Э, 8Н)Т — вектор управления, 8В, 8Э, 8Н — отклонения руля высоты, элеронов и руля направления соответственно. Матрица А размерности (5x5) содержит производные аэродинамических коэффициентов по фазовым координатам самолета, матрица В размерности (5X3) содержит производные аэродинамических коэффициентов по компонентам вектора управления. Нелинейная функция /(х) содержит члены, описывающие кинематическое и инерционное взаимодействие продольного и бокового движений самолета: /1 =— /2 = аюЛ.; /3 = с1<оушг; /4 = со^.шг;

/5 = с3 «у

На рис. 4 приведен пример расчета установившихся значений скорости крена <ах с помощью метода непрерывного продолжения решения при изменении двух параметров управления 8В, 8Э при 8Н = 0. Расчеты проведены при следующих значениях коэффициентов матриц А и В, функции /:.аи = 0,7; а22 = — 0,17; а31 = — 4; а33 = — 0,8; а42 = -- 3; аи = — 0,6; а45 = 0,055; а52=—15,5; а54= —0,7; #55 = — 2,35; Ьп = 0,002; &28 = —0,0004; 631 = —0,156; 642 == 0,0035; 643 = —0,05; Ь&2 = — 0,5; Ь53 — — 0,07; = — 0,9; с2 = 0,8 (остальные

коэффициенты принимались равными нулю).

Получаемые траектории образуют связную поверхность, которая в определенной области параметров управления характеризуется неоднозначной зависимостью <0^.(8,,, 8Э). Существуют области значений параметров управления, где имеются пять, три и одно решение. Штриховой линией на рис. 4 отмечены апериодически

неустойчивые режимы вращения, когда (1е^-^^>0. Исчезновение

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

1. Гоман М. Г., Захаров С. Б., Храброе А. Н. Аэродинамический гистерезис при отрывном обтекании удлиненных тел. — ДАН СССР, 1985, т. 282, № 1.

2. Zagaynov G. 1., О от an М. Q. Bifurcation analysis of aircraft critical flight regimes. — ICAS Proceedings, Toulose, France, 1984, vol. 1.

3. Mehra R. K., Carrol J. V. Bifurcation analysis of aircraft high-of-attack flight dynamics. — AIAA Paper N 80-1599.

4. Давиденко Д. Ф. Об одном новом методе численного решения систем нелинейных уравнений. — ДАН СССР, 1953, т. 88, № 4.

5. Kubicek М. Algorithm 502. Dependence of solution of nonli-ner system's on parameter ACM.— TOMS; 2 (1976).

6. Рыб am OB М. В. Непрерывные алгоритмы продолжения решений конечных уравнений, зависящих от параметров. — Автоматика и телемеханика, 1975, № 4.

7. Арнольд В. И. Дополнительные главы теории обыкновенных дифференциальных уравнений. — М.: Наука, 1978.

8. Pan С. Т., С h а о К. S. A computer-aided root-locus method.— IEEE Transaction on Aut. Control, AC-2, N 5, 1978.

9. Chao K. S.. Liu D. K.. Pan С. T. A systematic search method for obtaining multiple solutions of simultaneous nonlinear equations.-IEEE Transaction, 1975, vol. CAS-22, N 9.

10. Watson L. T. A globaly convergent algorithm for computing fixer points of Cs тарь. —Applied Mathematics and Computation, 1979, N 5.

Рукопись поступила 29jXlI 1984