Сравнительная характеристика методик решения систем нелинейных алгебраических уравнений в средах MathCAD и MathLab

D(x₁) r = ( a₁ - x₁ b₁ ) r = 0;

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

D(x_k) r = ( a_k - x_k b_k ) r = 0,

где матрицы a_i и b_i имеют размеры (N s_i /r_i) * (N - p), когда ранг матрицы A_p*N равен p.

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

Блок 1 выводит на экран параметры исходной СНАУ: матрицу системы, структуру вектора xA y, количество уравнений системы и максимальные степени переменных. В этом же блоке вычисляется нуль-пространства  
Q_0x = Ker(A_xy) и Q_0y = Ker(A_yx) размера 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. Схема алгоритма решения СНАУ

Для вычисления Q_0y по известной матрице Q_0x используется матрица перестановок P, связывающая векторы y x и  
x y, а следовательно, и соответствующие матрицы Q_0y = PQ_0x. Используя структуру вектора x y, формируются пучки D(x) = a_x - x b_x и D(y) = = a_y - y b_y.

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

(a_x - x b_x ) R_x = 0;

(a_y - y b_y) R_y = 0.                                                                   (2.17)

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

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

P(a - l b)Q = diag{D_r , D_i , D_f , D_l },                                          (2.18)

где D_r - сингулярная часть пучка, объединяющая блоки L_t^T размера t x  (t+1); D_i - регулярная часть пучка, соответствующая бесконечным собственным значениям; D_f  - регулярная часть пучка, соответствующая конечным собственным значениям; D_l - сингулярная часть пучка, объединяющая блоки L_s размера (s+1) x s.

Для того чтобы охарактеризовать структуры  сингулярных блоков D_l и D_r , а также регулярных блоков D_i и D_f   введем в рассмотрение следующие мультииндексы

a _l = { l₀ ... l_s ... l_u };                                                                   (2.19а)

a _r = { r₀ ... r_t ... r_w};                                                          (2.19б)

a _i = { j_i1 ... j_ik ... j_iv };                                                               (2.19в)

a _f = { j_f1 ... j_fk ... j_fv },                                                         (2.19г)

         т.е. - мультииндекс, указывающий число блоков типа L₀...L_s...L_u, соответствующих левым полиномиальным решениям; a _r - мультииндекс, указывающий число блоков типа L₀^T...L_t^T...L_w^T, соответствующих правым полиномиальным решениям; a _i  - мультииндекс, указывающий число клеток Жордана J_i1...J_ik...J_iv, соответствующих бесконечным собственным значениям; a _f  - мультииндекс, указывающий число клеток Жордана J_f1...J_fk...J_fv, соответствующих конечным собственным значениям с номерами f = {0, 1, ..., p}; u,v,w - индексы, задающие максимальные размеры сингулярных блоков и клеток Жордана.

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

a = { a₀ a₁ ... a_f ... a_p a_i a_l a_r },                                        (2.20)

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

Число блоков L_s и L_t^T равно, соответственно

, ,                                                 (2.21а)

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

, ,                                                 (2.21б)

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

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

N_r = n - rank( [a^H b^H] );                                                         (2.22б)

N₀ = m - rank( a );                                                                 (2.22в)

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

Эти соотношения не определяют состава  блоков каждого типа, но дают общую  интегральную оценку числа клеток Жордана  с нулевыми (N₀) и бесконечными собственными значениями (N_i ), общее количество правых (N_r) полиномиальных блоков L_t^T и левых (N_l) блоков L_s.

Выявление структуры мультииндексов a_r (2.19б) и a_i (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) - информационные векторы, описывающие структуру мультииндексов a_r и a_i; 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

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

n₁= m₁= N_i + N_r .                                                                  (2.24)

Далее нетрудно убедиться, что выполняются  соотношения  
N_r = (n₂+n₃+n₄) - (m₂+m₃), N_i = (m₁+m₂+m₃) - (n₂+n₃+n₄),  
что подтверждает условие (2.24). Из таблицы следует, что мультииндексы могут быть записаны следующим образом

a_r = {r₀, r₁, r₂, r₃} = { 0 n₂-m₂ n₃-m₃ n₄ };                                  (2.25а)

a_i = { j_i1, j_i2, j_i3 } = { m₁-n₂ m₂-n₃ m₃-n₄ }.                                   (2.25б)

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

size(D_f) = [n_f m_f] = [n - m - ].                          (2.26)

Размер пучка n_f = m_f   отличен от нуля, когда процедура завершается по одному из условий:  
либо n_k = 0, либо m_k = 0 и < m.

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

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

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

x(l ) = s_tl ^t + s_t-1l ^t-1 + ... +s₁l + s₀ ,                                             (2.27)

где s_t, s_t-1, ... s₁, s₀ - векторы размера n x 1, является правым полиномиальным решением пучка, если при любом l выполняется тождество (a - lb) x(l ) 0, из которого следует справедливость соотношений

bs_t = 0; bs_t-1 = as_t; ... bs₀ = as₁; as₀ = 0 .                             (2.28)

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

                                              (2.29а)

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

M_t S_t = 0 ,                                                                                  (2.29б)

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

Введем в рассмотрение блочно транспонированную  матрицу  
S_l  = S_t^B вида S_l =[S_t S_t-1...S_j...S₁ S₀], где size(Sj) = [N-p m ] для всех j от 1 до t.

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

S _lx l_x X = S _ly l_y Y,                                                                (2.30)

где S _lx=(S_t^x)^B, S_t^x = Ker(M_t^x), m _t = dim Ker(M_t^x), m_x = m _t/r_t;

S _ly=(S_v^y)^B, S_v^y = Ker(M_v^y), m _v = dim Ker(M_v^y), m_y = m _v/r_v;

l_x=[x^t ... 1]^T, l_y=[y^v ... 1]^T;

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

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

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

Шаг 3.1. Формирование матрицы M_t^x и M_v^y с учетом структуры мультииндексов a_r^x и a_r ^y по переменным x и y.

Шаг 3.2. Формирование вспомогательной СНАУ. Для этого вычисляются нуль-пространства  
S_t^x = Ker(M_t^x) и S_v^y = Ker(M_v^y) и по каждой переменной строится мультииндекс m , объединяющий размерности ядер  
m _t = dim Ker(M_t), зависящие от t

m = [ m ₀ m ₁ ... m _t ... m_w ],                                                                  (2.31)

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

m = [ 0 m ₁ m ₂ m ₃ ].                                                                       (2.31а)

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

                                                        (2.32)

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

size l _x = [t+1 1]; size l _y = [v+1 1];

size X = [m_x 1]; size Y = [m_y 1];

size S _lx = [N-p m_x* (t+1)]; size S _ly = [N-p m_y* (v+1)];

size[S _lx -S _ly ] = [N-p m_x*(t+1)+m_y *(v+1)].

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

pl* = m_x + m_y* (v+1) .                                                                           (2.33)