MATLAB.Exponenta
MATLAB и Simulink на русском
Технологии разработки и отладки
		сложных технических систем
 

Математика\SNAE Toolbox

В.Д.Борисевич, В.Г.Потемкин, С.П.Струнков "Пакет прикладных программ для решения систем нелинейных алгебраических уравнений от двух переменных"
Решение СНАУ

  В оглавление книги \ К следующему разделу \ К предыдущему разделу

2.2. Алгоритм решения СНАУ

Рассмотрим подробно алгоритм решения систем нелинейных алгебраических уравнений (СНАУ) от двух переменных вида:

Ap*N [x image003.gif (186 bytes) y] = 0,                                   (2.15)

где Ap*N - матрица коэффициентов системы уравнений;
x, y - векторы алгебраической структуры вида x = [xs...xj ...x 1]T, s = max j, ri = si +1;
XN - кронекеровское произведение векторов xA y.

Алгебраическая структура вектора Xk позволяет свести задачу определения решения относительно переменных
x1 .... xk к спектральной задаче для системы линейных сингулярных пучков вида [1, 3]

D(x1) r = ( a1 - x1 b1 ) r = 0;

. . . . . . . . . . . .                                        (2.16)

D(xk) r = ( ak - xk bk ) r = 0,

где матрицы ai и bi имеют размеры (N si /ri) * (N - p), когда ранг матрицы Ap*N равен p.

Обратимся к анализу алгоритма решения СНАУ; основные блоки алгоритма показаны на рис. 1 ниже.

Блок 1 выводит на экран параметры исходной СНАУ: матрицу системы, структуру вектора xA y, количество уравнений системы и максимальные степени переменных. В этом же блоке вычисляется нуль-пространства
Q0x = Ker(Axy) и Q0y = Ker(Ayx) размера N x (N-p).

Блок 1. Описание исходной СНАУ. Построение системы спектральных задач.
1.1. Параметры исходной СНАУ
1.2. Построение пучков D(x), D(y)
 
Блок 2. Вычисление структуры Кронекера
2.1. Характеристики структуры Кронекера
2.2. Структура Кронекера
2.3 Возможные источники решений
 
Блок 3. Вычисление решений системы спектральных задач (ССЗ)
3.1. Формирование Mtx, Mvy
3.2. Формирование вспомогательной СНАУ (ВСНАУ)
3.3. Редукция числа уравнений ВСНАУ
3.4. Формирование Qolx, Dlx, Qoly, Dly
3.5. Выделение регулярной части Df lx, Df ly
 
Блок 4. Формирование решений ССЗ
4.1. Вычисление спектра регулярных пучков
4.2. Решения ССЗ
 
Блок 5. Определение допустимых пар решений ССЗ
5.1. Редукция к СНАУ от одной переменной
5.2. Формирование допустимых конечных пар решений
 
Блок 6. Формирование решений СНАУ
6.1. Подстановка допустимых решений в исходную СНАУ
и оценка точности
6.2. Построение графика решений

Рис. 1. Схема алгоритма решения СНАУ

Для вычисления Q0y по известной матрице Q0x используется матрица перестановок P, связывающая векторы y image003.gif (186 bytes) x и
x image003.gif (186 bytes) y, а следовательно, и соответствующие матрицы Q0y = PQ0x. Используя структуру вектора x image003.gif (186 bytes) y, формируются пучки D(x) = ax - bx и D(y) = = ay - by.

Они позволяют составить систему уравнений

(ax - x bx ) Rx = 0;

(ay - y by) Ry = 0.                                                                   (2.17)

Базисным алгоритмом вычисления нуль-пространств является алгоритм SVD-разложения прямоугольной матрицы [10, 11].

Блок 2 вычисляет структуру Кронекера пучков D(x) и D(y). Структура Кронекера сингулярного пучка D(l ) = l  b определяется следующим образом [5-7]. Сингулярный пучок D(l ) размера m x n строго эквивалентен канонической форме Кронекера, если существуют постоянные невырожденные матрицы image004.gif (176 bytes)  Сmx m и image004.gif (176 bytes)  Сnx n такие, что справедливо [6]

P(a - l b)Q = diag{Dr , Di , Df , Dl },                                          (2.18)

где Dr - сингулярная часть пучка, объединяющая блоки LtT размера t x  (t+1); Di - регулярная часть пучка, соответствующая бесконечным собственным значениям; Df  - регулярная часть пучка, соответствующая конечным собственным значениям; Dl - сингулярная часть пучка, объединяющая блоки Ls размера (s+1) x s.

Для того чтобы охарактеризовать структуры сингулярных блоков Dl и Dr , а также регулярных блоков Di и Df   введем в рассмотрение следующие мультииндексы

a l = { l0 ... ls ... lu };                                                                   (2.19а)

a r = { r0 ... rt ... rw};                                                          (2.19б)

a i = { ji1 ... jik ... jiv };                                                               (2.19в)

a f = { jf1 ... jfk ... jfv },                                                         (2.19г)

т.е. a l - мультииндекс, указывающий число блоков типа L0...Ls...Lu, соответствующих левым полиномиальным решениям; a r - мультииндекс, указывающий число блоков типа L0T...LtT...LwT, соответствующих правым полиномиальным решениям; a i  - мультииндекс, указывающий число клеток Жордана Ji1...Jik...Jiv, соответствующих бесконечным собственным значениям; a f  - мультииндекс, указывающий число клеток Жордана Jf1...Jfk...Jfv, соответствующих конечным собственным значениям с номерами f = {0, 1, ..., p}; u,v,w - индексы, задающие максимальные размеры сингулярных блоков и клеток Жордана.

Cледуя работе [6], введем мультииндекс

a = { a0 a1 ... af ... ap ai al ar },                                        (2.20)

который определяет структуру Кронекера сингулярного пучка прямоугольных матриц.

Число блоков Ls и LtT равно, соответственно

image4619.gif (1033 bytes), image4620.gif (1036 bytes),                                                 (2.21а)

а число клеток Жордана, соответствующих нулевым и бесконечным собственным значениям, равно, соответственно

image4621.gif (1061 bytes), image4622.gif (1049 bytes),                                                 (2.21б)

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

Nl = m - rank( [a b] );                                                            (2.22а)

Nr = n - rank( [aH bH] );                                                         (2.22б)

N0 = m - rank( a );                                                                 (2.22в)

Ni = m - rank( b ).                                                                 (2.22г)

Эти соотношения не определяют состава блоков каждого типа, но дают общую интегральную оценку числа клеток Жордана с нулевыми (N0) и бесконечными собственными значениями (Ni ), общее количество правых (Nr) полиномиальных блоков LtT и левых (Nl) блоков Ls.

Выявление структуры мультииндексов ar (2.19б) и ai (2.19в) реализовано с помощью алгоритма RGSVD-разложения, который включен в ППП SNAE в виде М-функции sstruct на языке Matlab и обращение к которому имеет вид

[bk, ak, n(k), m(k)] = sstruct(b, a).                                        (2.23)

Входные параметры: a, b - матрицы пучка a - l b .

Выходные параметры: n(k), m(k) - информационные векторы, описывающие структуру мультииндексов ar и ai; ak, bk - матрицы, полученные после одного шага RGSVD-алгоритма.

На основе информационных векторов n(k), m(k) можно сформировать следующую таблицу

Таблица

k

nk

mk

r0

r1

r2

r3

ji1

ji2

ji3

1

n1

m1

0

     

m1-n2

   

2

n2

m2

 

n2-m2

     

m2-n3

 

3

n3

m3

   

n3-m3

     

m3-n4

4

n4

0

     

n4

     

Она справедлива при следующих условиях, практически всегда имеющих место при анализе пучков, возникающих при решении СНАУ от двух переменных

n1= m1= Ni + Nr .                                                                  (2.24)

Далее нетрудно убедиться, что выполняются соотношения
Nr = (n2+n3+n4) - (m2+m3), Ni = (m1+m2+m3) - (n2+n3+n4),
что подтверждает условие (2.24). Из таблицы следует, что мультииндексы могут быть записаны следующим образом

ar = {r0, r1, r2, r3} = { 0 n2-m2 n3-m3 n4 };                                  (2.25а)

ai = { ji1, ji2, ji3 } = { m1-n2 m2-n3 m3-n4 }.                                   (2.25б)

Регулярная часть пучка Df = af - l bf определяется квадратными матрицами af и bf размера

size(Df) = [nf mf] = [n -image4623.gif (926 bytes) m - image4624.gif (935 bytes)].                          (2.26)

Размер пучка nf = mf   отличен от нуля, когда процедура завершается по одному из условий:
либо nk = 0, либо mk = 0 и image4624.gif (935 bytes)< m.

Результатами работы блока 2 являются регулярные пучки Df (x) и Df (y), мультииндексы ar(x) и ar(y), а также список возможных источников решений.

Блок 3 определяет решения системы спектральных задач, которые включают регулярные и полиномиальные решения.

Полиномиальное решение для пучка D(l ) = a - l b определяется следующим образом [5]. Векторный полином вида

x(l ) = stl t + st-1l t-1 + ... +s1l + s0 ,                                             (2.27)

где st, st-1, ... s1, s0 - векторы размера n x 1, является правым полиномиальным решением пучка, если при любом l выполняется тождество (a - lb) x(l ) image005.gif (169 bytes) 0, из которого следует справедливость соотношений

bst = 0; bst-1 = ast; ... bs0 = as1; as0 = 0 .                             (2.28)

Рассматривая эту систему равенств как систему линейных однородных уравнений относительно элементов столбцов
st, st-1 ... s1, s0 можно сформировать систему вида

image4625.gif (1823 bytes)                                               (2.29а)

Уравнение (2.29a) может быть записано в следующем виде

Mt St = 0 ,                                                                                  (2.29б)

где Mt - матрица коэффициентов уравнения (2.29а) размера [(t+2)*N*s1/r1 (t+1)*(N-p)]; St - ядро Ker(Mt) размера [(t+1)*(N-p) m ], m = dim Ker(Mt).

Введем в рассмотрение блочно транспонированную матрицу
Sl  = StB вида Sl =[St St-1...Sj...S1 S0], где size(Sj) = [N-p m ] для всех j от 1 до t.

Для полиномиальных решений должно выполняться условие rx = ry, которое с учетом соотношений (2.27 - 2.29) можно записать в следующем виде

S lx lx image003.gif (186 bytes)X = S ly ly image003.gif (186 bytes)Y,                                                                (2.30)

где S lx=(Stx)B, Stx = Ker(Mtx), m t = dim Ker(Mtx), mx = m t/rt;

S ly=(Svy)B, Svy = Ker(Mvy), m v = dim Ker(Mvy), my = m v/rv;

lx=[xt ... 1]T, ly=[yv ... 1]T;

X - вектор размера [mx 1]; Y - вектор размера [my 1].

Число таких уравнений равно произведению величин Nrx и Nr y, вычисленных по соотношению (2.22б).

Cоотношения (2.28 - 2.30) определяют следующую последовательность шагов алгоритма в блоке 3:

Шаг 3.1. Формирование матрицы Mtx и Mvy с учетом структуры мультииндексов arx и ar y по переменным x и y.

Шаг 3.2. Формирование вспомогательной СНАУ. Для этого вычисляются нуль-пространства
St= Ker(Mtx) и Svy = Ker(Mvy) и по каждой переменной строится мультииндекс m , объединяющий размерности ядер
m = dim Ker(Mt), зависящие от t

m = [ m 0 m 1 ... m t ... mw ],                                                                  (2.31)

который можно сопоставить с мультииндексом ar (2.19б). Так например, структура мультииндекса дефектов матрицы Mt, соответствующая мультииндексу (2.25а) структуры правых полиномиальных решений, вычисленных в блоке 2, должна иметь вид

m = [ 0 m 1 m 2 m 3 ].                                                                       (2.31а)

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

image4626.gif (1386 bytes)                                                         (2.32)

где векторы и матрицы имеют следующие размеры

size l x = [t+1 1]; size l y = [v+1 1];

size X = [mx 1]; size Y = [my 1];

size S lx = [N-p mx* (t+1)]; size S ly = [N-p my* (v+1)];

size[S lx -S ly ] = [N-p mx*(t+1)+my *(v+1)].

Шаг 3.3. Редукция числа уравнений. Все отобранные уравнения удовлетворяют системе (2.32), причем число уравнений этой системы может быть уменьшено, если выполняются некоторые условия. Системе (2.32) соответствует пучок Dlx  размера [mx* t mx* (t+1)+my* (v+1)-pl], где pl - число уравнений системы. Условие регулярности пучка определяется соотношением

pl* = mx my* (v+1) .                                                                           (2.33)

Если оказывается, что текущее pl > pl*, то число уравнений, используемых для формирования СНАУ (2.32), следует сократить до величины pl*. Если описанная редукция возможна, то на следующем шаге pl = pl*, а сформированный пучок будет регулярным; если - нет, то число уравнений остается равным pl = N-p, а пучок - сингулярным.

Шаг 3.4. Вычисление

image4627.gif (1405 bytes)                                                            (2.34)

размера [Nl Nl-pl], где Nl = mx*(t+1)+my*(v+1).

Выделить из матрицы Q0l подматрицу Q0l x, соответствующую подблоку l x image003.gif (186 bytes) X и имеющую размер [mx*(t+1) Nl-pl]. Сформировать пучок Dlx alx xblx размера [m *t Nl-pl]. Если редукция на шаге 3.3 была успешной, то сформированный пучок - регулярный и можно переходить к блоку 4, не выполняя шаг 3.3.

Шаг 3.5. Выделение из сингулярных пучков Dlx и Dly, полученных на шаге 3.4, регулярных частей
Dflx aflx-xbflx и
Dflaflxbfl y.

Блок 4 формирует решения ССЗ. Вычислить спектры регулярных пучков Df(x) и Df(y), полученных в блоке 2 и Dflx и Dfl y, полученных в блоке 3, используя процедуру eig(a, b) [12] решения обобщенной проблемы собственных значений. Результатом является множество решений {xi} и {yi} по переменным x и y.

Блок 5 определяет допустимых пар решений ССЗ. Этот блок использует следующую последовательность шагов.

Шаг 5.1. Каждое из вычисленных в блоке 4 решение подставляется в исходную СНАУ и таким образом выполняется редукция к СНАУ от одной переменной. Вычисляются решения этой СНАУ и формируются все допустимые пары решений, включая решения типа Inf и NaN.

Шаг 5.2. Из найденных на шаге 5.2 пар решений извлекаются только конечные пары {xi, yi}.

Блок 6 формирует решения СНАУ. Этот блок включает два шага.

Шаг 6.1. Найденные конечные пары решений подставляются в исходную СНАУ и вычисляются невязки для каждого уравнения системы. Сумма модулей этих невязок служит оценкой точности решений СНАУ.

Шаг 6.2. Построение графика решений, на который выводится графический образ исходной матрицы и множество всех значений {xi } и {yi }, как действительных, так и комплексных.

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

Алгоритм решения СНАУ в полном объеме реализован в среде системы научных расчетов Matlab 5 [13].

В качестве теста рассмотрим пример решения системы нелинейных алгебраических уравнений от двух переменных вида

image4628.gif (1255 bytes)

Протокол решения приведен ниже и не нуждается в дополнительных комментариях. Последовательность решения соответствует всем вышеописанным блокам алгоритма.

   В оглавление книги \ К следующему разделу \ К предыдущему разделу


Поиск по сайту:

Система Orphus

Яндекс.Метрика