УДК 519.61:621.3

В.П. ВОЛОБОЕВ*, В.П. КЛИМЕНКО*

ОБ ОДНОМ ПОДХОДЕ К РЕШЕНИЮ ПЛОХО ОБУСЛОВЛЕННОЙ СИСТЕМЫ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ, ОПИСЫВАЮЩИХ ФИЗИЧЕСКИЙ ОБЪЕКТ

Институт проблем математических машин и систем НАН Украины, Киев, Украина

Анотація. Показано, що вірогідність результатів моделювання фізичних об"єктів, дискретна модель яких описується системою лінійних алгебраїчних рівнянь (СЛАР), залежить не від поганої обумовленості матриці, а від некоректного вибору змінних СЛАР на етапі складання рівнянь методом вузлових потенціалів або його аналогів, а сам метод є окремий випадок методу коректної постановки завдання. Запропоновано методику перевірки на коректність СЛАР, складеної методом вузлових потенціалів, що має невироджену й симетричну матрицю, і якщо необхідно перетворення її до коректного виду.

Ключові слова: система, моделювання, некоректне завдання, погана обумовленість, система лінійних алгебраїчних рівнянь, метод вузлових потенціалів, метод коректної постановки завдання, перевірка на коректність.

Аннотация. Показано, что достоверность результатов моделирования физических объектов, дискретная модель которых описывается системой линейных алгебраических уравнений (СЛАУ), зависит не от плохой обусловленности матрицы, а от некорректного выбора переменных СЛАУ на этапе составления уравнений методом узловых потенциалов или его аналогами, а сам метод есть частный случай метода корректной постановки задачи. Предложена методика проверки на корректность СЛАУ, составленной методом узловых потенциалов, имеющей невырожденную и симметричную матрицу, и если необходимо преобразование её к корректному виду.

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

Abstract. The paper shows that reliability of results of simulation of the physical objects, which discrete model is described by a system of the linear algebraic equations (SLAE) depends not on poor-conditioned matrix but on an incorrect choice of variable SLAE at generation of equations stage by a node potential method or its analogues, and the method is a special case of a method of correct statement of a problem. It was suggested the check-out method on a correctness of SLAE, made by a node potential method, having nonsingular and a symmetric matrix and if it is necessary its transformation to a correct form.

Keywords: system, simulation, incorrect problem, poor-conditioned, system of the linear algebraic equations, node potential method, method of correct statement of a problem, check-out on a correctness.

1. Введение

Многие задачи моделирования физических (технических) объектов сводятся к решению систем линейных алгебраических уравнений (СЛАУ). Поскольку все вычисления при решении таких систем выполняются с конечным числом значащих цифр, то точность может значительно теряться из-за ошибок округления. Плохо обусловленной (неустойчивой) системой или в более общей формулировке - некорректно поставленной задачей принято считать ту задачу, которая при фиксированном уровне ошибок входных данных и точности вычислений не гарантирует в решении никакой точности. В качестве априорной наихудшей оценки возможных ошибок решения СЛАУ используется число обусловленности . Как следует из литературы, разработка методов решения некорректно поставленных задач рассматривается как чисто математическая задача, в которой не учитываются особенности физических (технических) объектов, несмотря на то, что численное решение многих задач математической физики и математического моделирования сложных физических процес-

© Волобоев В.П., Клименко В.П., 2014

сов и технических систем является неиссякаемым источником задач линейной алгебры. Для перечисленного класса задач при разработке методов решения не рассматривается этап составления СЛАУ, на котором тем или иным способом можно учесть особенности конкретной задачи. В том, что этот этап необходимо учитывать, подтверждают результаты следующих работ.

Прежде всего, следует отметить работу , где приведены примеры матриц, для которых потеря точности при решении СЛАУ невелика, а величина числа обусловленности огромна, то есть показано, что общепринятый критерий априорной оценки точности решения СЛАУ по числу обусловленности есть необходимый, но недостаточный. Совершенно новый подход к решению некорректно поставленной задачи предложен в работах . Он заключается в том, что с целью повышения точности решения СЛАУ даже при большой величине числа обусловленности на этапе описания дискретной модели физического объекта предлагается корректно составлять СЛАУ. Это означает не только то, что такие матрицы существуют, как об этом сообщалось в работе , но и то, что предложен метод корректного составления матрицы СЛАУ, описывающей дискретную модель объекта. Метод составления матрицы СЛАУ рассмотрен применительно к задачам моделирования поведения электрических цепей , энергосистем , стержневых систем механики и эллиптических уравнений матфизики .

Суть данного метода заключается в том, что, в отличие от существующих методов, при формировании СЛАУ целенаправленным выбором переменных учитываются параметры дискретной модели физического объекта. Следует заметить, что метод применим только к тем объектам, топология дискретной модели которых представлена графом.

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

В данной работе будет предложен метод корректировки некорректно поставленной задачи для случая, когда топология дискретной модели не представлена в виде графа. При разработке метода учитывается тот факт, что общепринятый метод описания дискретных моделей задач математической физики и сложных физических процессов и технических систем (метод узловых потенциалов ) есть частный случай метода корректного составления матрицы СЛАУ.

2. Связь между точностью решения СЛАУ, описывающей дискретную модель объекта, и методом составления уравнений

Академик Воеводин В.В. в работе показал, что наибольшая точность результатов решения СЛАУ методом Гаусса достигается в случае применения метода с выбором главного элемента. На основе этой идеи было опубликовано огромное число работ. Однако решение практических задач показало, что точность решения СЛАУ, особенно в случае плохо обусловленных матриц, значительно теряется из-за ошибок округления, то есть для повышения точности результатов на этапе решения недостаточно одного применения метода Гаусса с выбором главных элементов.

Дальнейшим развитием этой идеи есть метод, предложенный в работе , где предлагается на этапе составления описания дискретной модели объекта формировать диагональные элементы матрицы как главные. Для этого при составлении описания используется дополнительная информация, а именно, параметры дискретной модели. Эффективность данного подхода, а именно, зависимость точности решения СЛАУ, описывающей дискрет-

ISSN 1028-9763. Математичні машини і системи, 2014, № 4

ную модель объекта, от метода составления уравнений будет продемонстрирована на модельном примере. Ниже будет рассмотрено составление описания модельного примера методом, описанным в , с выбором главного элемента и без выбора и его решение.

В качестве модельного примера выбрана электрическая цепь, приведенная на рис. 1.

Рис. 1. Электрическая цепь

Известно , что обусловленность СЛАУ, описывающей электрическую цепь, зависит от диапазона разброса величин проводимостей (сопротивлений) компонент цепи. Выбранный диапазон изменения проводимостей компонент электрической цепи, равный 15 порядкам, обеспечивает плохую обусловленность СЛАУ и тем самым, как принято считать, некорректность задачи. На примере расчета потенциала узла 2 (напряжение на компоненте G2) будет анализироваться зависимость достоверности результатов расчета от способа формирования диагонального элемента при составлении описания электрической цепи.

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

где U i - напряжение, падающее на компоненте, I - ток, протекающий через компоненту, Gt - проводимость компоненты.

Для описания графа электрической цепи и, соответственно, уравнений на основе законов Кирхгофа применяются топологические матрицы контуров и сечений. Граф цепи совпадает с электрической цепью. Составление топологических матриц контуров и сечений включает выбор дерева графа цепи и составление контуров для выбранного дерева. Дерево графа электрической цепи выбирается таким образом, чтобы все источники напряжения включались в дерево, а все источники тока в хорды. Элементы в векторах напряжений U и токов I компонент цепи группируются на входящие в дерево (индекс Д), то есть ветви и хорды (индекс X), таким образом:

Контуры образуются присоединением хорд к дереву графа схемы. В этом случае

топологическая матрица контуров имеет вид

где 1 - единичная подматрица хорд, t

Обозначает транспонирование матрицы, а топологическая матрица сечений - вид |1 -F , где 1 - единичная подматрица ветвей. Как следует из , диагональные члены матрицы

ISSN 1028-9763. Математичні машини і системи, 2014, № 4

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

их =-ґид, (3)

Переменные составляемой системы уравнений выбираются из напряжений и/или токов компонент в результате анализа основной системы уравнений. В случае выбора в качестве переменных напряжений компонент, входящих в ветви дерева, компонентные уравнения (1) и уравнения (3), (4) можно преобразовать к следующему виду:

Gд U д - F(Gx (- FUд)) = 0.

Ниже будет приведено составление уравнений для модельного примера. Вначале составляется описание электрической цепи таким образом, чтобы диагональные члены матрицы были главными. Этому требованию удовлетворяет набор компонент E1, G6, G3, G2, входящих в дерево (на рис. 1 ветви дерева выделены жирной линией). Выбранному дереву соответствуют следующие векторы напряжений, токов компонент:

и топологические матрицы

Уравнение (5), с учетом (6), (7) и компонентных уравнений после преобразований, имеет следующий вид:

- (G4 + G5) (G4 + G5) G1 + G2 + G4 + G5

СЛАУ (8) есть плохо обусловленная, так как собственные числа матрицы \= 1,5857864376253, Я2 = 5,0E +14+j5,0E +14, А, = 5,0E +14 - j5,0E +14. Для того, чтобы определить, как зависит точность результатов решения системы от выбора варианта составления уравнений, расчет потенциала Uq узла 2 будет выполнен в общем виде:

ISSN 1028-9763. Математичні машини і системи, 2014, № 4

(g1+g2 +g4 +g5)-

Из анализа вычислительного процесса (9-11) следует, что, несмотря на большой диапазон изменения величин проводимостей (15 порядков), нет жестких требований к конечной точности представления чисел как при составлении уравнений, так и при их решении. Для получения достоверного результата достаточно выполнить вычислительный процесс составления и решения СЛАУ с точностью представления чисел в две значащие цифры.

Следует заметить, что в СЛАУ (8) диагональный элемент второй строки (столбца) матрицы G+G4+G5I значительно больше (на 15 порядков) суммы остальных членов

строки (столбца) | G4 + 2G51. Это означает, что, приняв UG = 0, можно упростить СЛАУ

(8), сохранив достоверность результатов. В эпоху ручного счета этому приему соответствовало объединение узла 2 с 3 (рис. 1).

Во втором случае (без выбора диагонального элемента главным) достаточно выбрать в дерево компоненты Ex, G6, G4, G2 (на рис. 1 ветви дерева помечены прерывистой

линией). Падения напряжений на этих компонентах соответствуют узловым потенциалам 1, 4, 3, 2, отсчитываемым от нулевого узла. Это означает, что при таком выборе компонент в дерево метод корректного составления матрицы СЛАУ совпадает с методом узловых потенциалов. Выбранному дереву и хордам соответствуют следующие векторы напряжений, токов компонент:

U Д = UG UG G4 , Ux = G1 UG3 UG G Д G ig G4 , Ix = G1 IG3 IG

UG G2 G5 ig G2 G5

и топологических матриц

Уравнение (5), с учетом (12), (13) и компонентных уравнений, примет следующий

ISSN 1028-9763. Математичні машини і системи, 2014, № 4

G5 + G6 -G5 0 UG G6 0

G5 G3 + G4 + G5 -G3 Uo. = 0

0 - G3 G1 + G2 + G3 Uo2 G1E1

Система уравнений (14) есть плохая обусловленная, так как имеет следующие собственные числа матрицы: 1 = 1.0,1 =1015 +у1015,1 =1015-/1015 . Как в первом варианте примера, будет выполнен расчет потенциала UG узла 2 в общем виде:

(G + G + G)------------

V 3 4 У (G + G)

+ (G1 + G2 + G3)

3 4 5" (G5 + G6)

Из анализа вычислительного процесса решения системы уравнений (15-17) следует, что достоверность результатов зависит как при составлении, так и решении уравнений от конечной точности представления чисел. Так, если вычислительный процесс решения системы (15-17) выполняется с точностью меньше 15 значащих цифр, то результат будет

1015 +1015 ~ o,

а в случае, когда с точность больше 15 значащих цифр, будет

1030 + 2*1015 +1030 + %+ 3/1015)

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

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

Диагональные элементы матрицы по модулю больше остальных элементов, как в строках, так и столбцах, независимо от того, матрица составлена с или без выбора максимальных диагональных. Различие заключается только в том, насколько диагональные элементы больше недиагональных. Это означает, что решение такого типа СЛАУ методом Г аусса с выбором главного элемента не повышает точности результатов для данного класса задач.

ISSN 1028-9763. Математичні машини і системи, 2014, № 4

Конечное число значащих цифр, используемых при решении методом Г аусса, существенным образом зависит от того, матрица составлена с или без выбора максимальных диагональных элементов. Отличие одного варианта задачи от другого заключается только в том, что на этапе составления уравнений в одном случае компонента с максимальной проводимостью выбирается в дерево и тем самым напряжение этой компоненты выступает в качестве переменной составляемой СЛАУ. Проводимость этой компоненты участвует только при формировании диагонального элемента матрицы. В другом случае эта компонента попадает в хорды. Как следует из уравнения (3), напряжение компоненты определяется через напряжения компонент дерева. Из уравнения (4) следует, что проводимость компоненты участвует в формировании элементов строк и столбцов и тем самым проводимость хорды определяет величину этих элементов матрицы.

3. Преобразование матрицы СЛАУ, составленной методом узловых потенциалов, к виду, соответствующему корректной постановке

При численном решении задач математической физики и математического моделирования сложных физических процессов и технических систем для составления СЛАУ, описывающих дискретные модели этих задач, в основном применяется метод узловых потенциалов или его аналоги. Отличительной чертой этого метода есть то, что в качестве переменных СЛАУ применяются потенциалы расчетной схемы дискретной модели, отсчитываемые от базового узла к остальным узлам, простой алгоритм составления уравнений, слабо заполненная матрица СЛАУ. Платой за такую эффективность может быть некорректность поставленной задачи. Учитывая, что метод узловых потенциалов есть всего лишь один из вариантов метода корректной постановки задачи, то некорректно поставленную задачу можно откорректировать, применив преобразование матрицы. Ниже будет рассмотрен алгоритм преобразования некорректно составленной методом узловых потенциалов задачи.

Из всего многообразия физических объектов будут рассмотрены только те объекты, линейная дискретная модель которых описывается СЛАУ с невырожденной и симметричной матрицей.

3.1. Алгоритм преобразования матрицы

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

Как следует из исследования, выполненного выше, некорректность поставленной задачи на уровне составленной СЛАУ возникает только в том случае, если хотя бы один из недиагональных элементов строки будет значительно больше параметра дискретной модели, входящего только в диагональный элемент. Ниже приводится методика проверки на корректность составленной СЛАУ.

Пусть СЛАУ имеет вид

где x - вектор узловых потенциалов (узловых воздействий), у - вектор внешних потоков, A - матрица вида

ISSN 1028-9763. Математичні машини і системи, 2014, № 4

а11 а1і a1j a1n

аі1 а,і aj ain , (21)

aJ1 an1 аі aJJ ann

где n - размер матрицы. Элементы матрицы удовлетворяют следующим требованиям:

аи > 0, а. < 0, а. = а]г,1 < i < n, 1 < j < n при j Ф і. (22)

Ниже будет рассмотрена проверка на корректность і -ой строки матрицы и, если необходимо, её корректировка.

Прежде всего определяется параметр дискретной модели ait, входящий только в диагональный элемент і -ой строки матрицы,

Считается, что і -ая строка матрицы корректно составлена, если параметр ait удовлетворяет условию

1 < j < n, при j Ф і.

При невыполнении условия (24) выполняется корректировка і -ой строки. Вначале из недиагональных элементов выбирается наибольший. Пусть это будет j -ый элемент і -ой строки. Нетрудно убедиться, что в силу специфики составления матрицы (условие (22)) параметр дискретной модели, который участвует в формировании элементов о. и а.^ і -ой и j -ой строк, входит как составная часть в элементы aii и а. . Суть корректировки і -ой строки заключается в преобразовании і -ой и j -ой строк матрицы таким образом, чтобы величина элемента а. входила только в элемент aii. Нетрудно убедиться, что, представив переменную xi в виде

X = xj + xj (25)

и выполнив следующее преобразование элементов j -ого столбца матрицы СЛАУ

о = аі. + ai, 1 < 1 < n , (26)

получим новый j -ый столбец матрицы, в котором преобразованные элементы а. и а. не содержат параметр дискретной модели, формировавший элементы а. и а. .

На следующем шаге выполняется преобразование j -ой строки по формуле

aji = а.і + аіі, 1 < l < n . (27)

Элементы a і преобразованной j -строки уже не содержат параметр дискретной модели, соответствующий элементу а і .

ISSN 1028-9763. Математичні машини і системи, 2014, № 4

Проверка на корректность матрицы СЛАУ и корректировка некорректных строк выполняются для всей матрицы. В данной работе рассмотрен только подход к построению алгоритма преобразования матрицы к корректному виду. Вопросы, связанные с разработкой эффективно работающего алгоритма преобразования матрицы к корректному виду, в этой работе не рассматриваются. Ниже будет приведен пример преобразования матрицы СЛАУ (14), составленной методом узловых потенциалов.

3.2. Демонстрационный пример

Прежде всего следует отметить, что матрица (14) симметричная и невырожденная. Коэффициенты матрицы удовлетворяют условию (22). Узловые потенциалы соответствуют падению напряжений на компонентах

U4 = UG^, U3 = UG,U2 =UG

Учитывая (28), СЛАУ (14) можно представить следующим образом:

G5 + G6 - G5 0 U 4 0

G5 G3 + G4 + G5 - G3 U3 = 0

0 - G3 G + G2 + G3 U2 GA

Проверка на корректность матрицы включает следующие операции.

Определение по формуле (23) параметра дискретной модели ait, входящего только

в диагональный элемент. Для первой строки матрицы это будет G6, для второй строки G4 и для третьей - (Gl + G2) .

Проверка строк матрицы на корректность выполняется в соответствии с формулой (24). В результате этой проверки оказывается, что вторая строка не удовлетворяет требованию корректности, так как (G4 = 1) ^ (G3 = 1015) . Параметр G3 входит также в третью строку матрицы, поэтому в соответствии с формулой (25) выбирается представление переменной U3 в виде

U3 = U2 + U23 , (30)

В результате преобразования элементов 3 -го столбца, в соответствии с формулой (26), получаем матрицу (29) следующего вида:

G5 + G6 - G5 - G5

G5 g3 + g4 + g5 g4+g5

а после преобразования третьей строки, в соответствии с формулой (27), матрица (31) будет иметь вид

(G5 + G6) - G5 - g5 U 4 0

G5 (G3 + G4 + G) (G4 + G5) U 23 = 0 . (32)

G5 (G4 + g5) (G + G2 + G4+g5) U2 G E

СЛАУ (32) удовлетворяет требованию корректности, поэтому корректировка считается завершенной. Переменные СЛАУ (32) соответствуют переменным СЛАУ (8), то есть в

ISSN 1028-9763. Математичні машини і системи, 2014, № 4

результате преобразования в дерево были выбраны те же компоненты, что и в методе корректной постановки задачи. Из сравнения СЛАУ (8) и (32) следует, что недиагональные элементы матрицы (32) второго столбца и второй строки отличаются по знаку от матрицы (8). Это есть результат того, что при преобразовании матрицы (14) было выбрано направление тока компоненты G3, противоположное направлению, выбранному при составлении СЛАУ (8). Выполнив замену переменной U23 на U23 = -U23 и поменяв во втором уравнении знаки элементов на противоположные, получим матрицу (8).

4. Заключение

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

В отличие от существующего подхода к разработке методов решения некорректных задач, в данной работе предлагается некорректно поставленную задачу (плохо обусловленную) приводить к корректному виду. Показано, что не плохая обусловленность матрицы затрудняет получение достоверных результатов при решении СЛАУ, описывающих дискретные модели физических объектов, а некорректный выбор переменных СЛАУ на этапе составления уравнений, а метод узловых потенциалов и его аналоги, которые применяются для составления СЛАУ, описывающих дискретную модель, есть частный случай метода корректной постановки задачи. Предложена методика проверки на корректность составленной методом узловых потенциалов СЛАУ для случая, когда матрица СЛАУ есть невырожденная и симметричная. Рассмотрен алгоритм преобразования матрицы к корректному виду.

СПИСОК ЛИТЕРАТУРЫ

1. Калиткин Н.Н. Количественный критерий обусловленности систем линейных алгебраических уравнений / Н.Н. Калиткин, Л.Ф. Юхно, Л.В. Кузьмина // Математическое моделирование. - 2011. Т. 23, № 2. - С. 3 - 26.

2. Волобоев В.П. Об одном подходе к моделированию сложных систем / В.П. Волобоев, В.П. Клименко // Математичні машини і системи. - 2008. - № 4. - С. 111 - 122.

3. Волобоев В.П. Об одном подходе к моделированию энергосистем / В.П. Волобоев, В.П. Клименко // Математичні машини і системи. - 2009. - № 4. - С. 106 - 118.

4. Волобоев В.П. Механика стержневых систем и теория графов / В.П. Волобоев, В.П. Клименко // Математичні машини і системи. - 2012. - № 2. - С. 81 - 96.

5. Волобоев В.П. Метод конечных элементов и теория графов / В.П. Волобоев, В.П. Клименко // Математичні машини і системи. - 2013. - № 4. - С. 114 - 126.

6. Пухов Г.Е. Избранные вопросы теории математических машин / Пухов Г.Е. - Киев: Изд-во Академии наук УССР, 1964. - 264 с.

7. Сешу С. Линейные графы и электрические цепи / С. Сешу, М.Б. Рид. - М.: Высшая школа, 1971. - 448 с.

8. Зенкевич О. Конечные элементы и аппроксимация / О. Зенкевич, К. Морган. - М.: Мир, 1986. -318 с.

9. Воеводин В.В. Вычислительные основы линейной алгебры / Воеводин В.В. - М.: Наука, 1977. -304 с.

10. Теоретические основы электротехники: учебник для ВУЗов / К.С. Демирчян, Л.Р. Нейман, Н.В. Коровкин, В.Л. Чечурин. - . - Питер, 2003. - Т. 2. - 572 с.

Выходные данные сборника:

РЕШЕНИЕ ПЛОХО ОБУСЛОВЛЕННЫХ РАЗРЕЖЕННЫХ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ С ПОМОЩЬЮ КРЫЛОВСКОГО ПОДПРОСТРАНСТВА

Гусева Юлия Сергеевна

студент Самарского государственного аэрокосмического университета имени С.П. Королева,

г. Самара

Е-mail:

Гоголева Софья Юрьевна

доцент Самарского государственного аэрокосмического университета имени С.П. Королева,

г. Самара

Е-mail:

Введение

Математические модели многих практических задач приводят к решению СЛАУ с большими и разреженными матрицами коэффициентов, в которых большинство элементов равны нулю. Приписывание матрице свойства разреженности эквивалентно утверждению о существовании алгоритма, использующего её разреженность. Когда большая доля коэффициентов матрицы состоит из нулей, вполне очевидно, что нам не хотелось бы хранить в памяти компьютера все эти нули. Поэтому матричные алгоритмы должны проектироваться таким образом, чтобы обрабатывались только ненулевые элементы и чтобы на основании предварительного знания о расположении ненулевых элементов избегались операции типа сложения с нулем или умножения на нуль. Таким образом число операций, производимых машиной при исполнении алгоритма, пропорционально числу ненулевых элементов, а не числу всех элементов матрицы. Серьезную проблему при работе с разреженными матрицами представляет численная устойчивость .

Когда методы типа гауссова исключения требуют слишком много времени или памяти для решения систем уравнений используются итерационные методы. При решении плохо обусловленных разрежен­ных СЛАУ возникает необходимость выбора метода, который позволит получить при решении точный результат и наименьшее заполнение (возникновение новых ненулевых элементов) . Наиболее эффективными и устойчивыми среди итерационных методов решения таких систем уравнений являются так называемые проекционные методы, и особенно тот их класс, который связан с проектированием на подпространства Крылова. Эти методы обладают целым рядом достоинств: они устойчивы, допускают эффективное распараллеливание, работу с различными строчными (столбцовыми) форматами и предобуславливателями разных типов.

Постановка задачи

Рассмотрим систему линейных алгебраических уравнений

Где: - плохо обусловленная разреженная матрица,

.

В данной работе проводится сравнительный анализ итерацион­ных методов для решения плохо обусловленных разреженных СЛАУ. В качестве исследуемых методов выбраны: метод сопряженных градиентов (CG) , метод минимальных невязок (MinRes), сдвоенный метод сопряженных направлений (CGS), квазиминимальных невязок (QMR) .

В вопросах выбора того или иного способа решения СЛАУ важно учитывать структуру матрицы A . Это связано с тем, что не каждый метод дает возможность получить гарантированный результат для определенной системы линейных уравнений.

Таким образом, критерием сравнения итерационных методов решения СЛАУ будут: погрешность результатов, скорость сходимости, структура матрицы.

Результаты численных исследований показали, что для решения СЛАУ с матрицей A являющейся симметричной/несимметричной и хорошо обусловленной к нормальным уравнениям лучше применять метод CG. Если матрица А- симметричная и плохо обусловленная, то лучшую сходимость показал метод MinRes. Для А- несимметрич­ной, плохо обусловленной - метод квазиминимальных невязок .

Для улучшения скорости сходимости итерационных методов используют предобуславливание матрицы системы. Оно заключается в том, что подбирается такая матрица предообуславливания, что при этом процедура решения СЛАУ является не слишком трудоемкой и численно устойчивой . Правильный выбор предобуславливателя, зависящий от конкретной задачи, может в огромной степени ускорить сходимость . В действитель­ности, хороший предобуславливатель часто необходим для того, чтобы итерационный метод вообще сходился.

В данной работе было рассмотрено несколько видов предобус­лавливания для метода квазиминимальных невязок с разреженными плохо обусловленными матрицами: левое и правое предобус­лавливание с использованием QR - разложения, левое и правое предобуславливание с использованием LU - разложения, а также с использованием модификации LU - разложения .

Таблица 1.

Сравнение относительной погрешности предобуславливателей

Матрица	LU - разложение		LU - разложение(модификация)	QR - разложение
	(левое)	(правое)		(левое)	(правое)

Заключение

В статье был рассмотрен метод квазиминимальных невязок применительно к решению разреженных плохо обусловленных СЛАУ и различные варианты выбора предообуславливателя. Метод квазими­нимальных невязок, основанный на использовании предобуслав­ливателя, полученного с помощью модификации LU- разложения дал наилучший результат по численной устойчивости.

Список литературы:

1.Голуб Дж., Ван Лоун Ч. Матричные вычисления/ Под ред. В.В. Воеводина. - М.: «Мир», 1999. - 548 с.

2.Деммель Дж. Вычислительная линейная алгебра. Теория и приложения / Пер. с англ.Х.Д. Икрамова. - М.: «Мир», 2001. - 430 c.

3.Писсанецки С. Технология разреженных матриц/ Под ред. Х.Д. Икрамова - М.: «Мир», 1988. - 410 с.

4.Станкевич, И.В. Хранение и использование разреженных матриц в конечно- элементной технологии. Журнал «Наука и Образование». - 2005. - 10 октября.

5.Тьюарсон Р. Разреженные матрицы/ Под ред. Х.Д. Икрамова. - М.: «Мир», 1977. - 172 с.

6.Bucker Martin, Basermann Achim. A comparison of QMR, CGS and TFQMR on a distributed memory machine / Bucker Martin //Mathematics of computation. - 1994 - 31 may

7.Harwell-Boeing Collection - [Электронный ресурс] – Режим доступа. - URL: http://math.nist.gov/MatrixMarket/data/Harwell-Boeing/ (дата обращения: 15.12.2012)

8.Roland W. Freund, Noel M. Nachtigal. QMR: a Quasi-Minimal Residual Method for Non-Hermitian Linear Systems / Roland W. Freund, Noel M. Nachtigal // Journal Math. - 1991. - №60. - p. 315-339.

9.Saad, Y. Iterative methods for sparse linear systems / Y. Saad. // SIAM. - 2003. - 447 p.

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

Числом обусловленности линейного оператора A , действующего в нормированном пространстве а также числом обусловленности системы линейных уравнений Ax = у назовем величину

Таким образом, появляется связь числа обусловленности с выбором нормы.

Предположим, что матрица и правая часть системы заданы неточно. При этом погрешность матрицы составляет dA , а правой части - dу . Можно показать, что для погрешности dx имеет место следующая оценка ():

В частности, если dA = 0, то

При этом решение уравнения Ax = у не при всех у одинаково чувствительно к возмущению dу правой части.

Свойства числа обусловленности линейного оператора:

1.

причем максимум и минимум берутся для всех таких x , что Как следствие,

3

где и - соответственно минимальное и максимальное по модулю собственные значения матрицы A . Равенство достигается для самосопряженных матриц в случае использования евклидовой нормы в пространстве

4.

Матрицы с большим числом обусловленности (ориентировочно ) называются плохо обусловленными матрицами . При численном решении систем с плохо обусловленными матрицами возможно сильное накопление погрешностей, что следует из оценки для погрешности dx . Исследуем вопрос о погрешности решения, вызванной ошибками округления в ЭВМ при вычислении правой части. Пусть t - двоичная разрядность чисел в ЭВМ. Каждая компонента вектора правой части округляется с относительной погрешностью Следовательно,

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

Итак – принципиально остаются две проблемы –

1 .не обеспечивается обоснованная сходимость алгоритма к единственной (в случае модельного примера- истинной) структуре и

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

Для АШР даже в случае применения для МНК оценки процедуры Грамма-Шмидта не разрешается вопрос о единственности модели – просто оценки параметров становятся наиболее точными и несмещенными

Т.о. гарантированное нахождение всего множества подходящих решений в реальных задачах (при - количество линейных входных аргументов и степени ПП p >3) получим только после полного

перебора всех подструктур полной структуры как в методе всех регрессий (у Дрейпера и Смита). Тогда мы найдем всете модели , в которых все аргументы входят с уровнем значимости не менее чем заданный. Со всеми выше описанными проблемами – а какая же из них, из этого множества та, которая действительно наша.

Можно еще добавить камень в огород АШР о неиспользуемой возможности вариации уровнем значимости для учета уровня шума в данных

Именно эту проблему предлагает решать МГУА с помощью введения понятия внешних критериев .

Необходимое примечание .

при все типы АШР МВИ, МГУА, другие целесообразные подходы, дают практически одинаково эффективные (или неэффективные) решения. Кривые критериев одинаково асимптотически стремятся к некоторому ненулевому уровню, при подходе к которому и определяется единственная модель.

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

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

И наиболее эффективный подход к решению структурно-параметрического синтеза при данных условиях демонстрирует МГУА

Как видим нарушение уже первого условия порождает необходимомость разрешения проблемы множественности моделей не прибегая к процедуре полного перебора –надо предложить какой-то принцип, позволяющий найти путь к истинной или квазиистинной модели без полного перебора претендентов моделей.

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

Основная проблема – проблема необоснованности выбора структуры модели классическими АШГ многократно обостряется в связи с тем что порог используемый критерием Фишера в виде уровня значимости

на самом деле регулирует не только риск ошибки

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

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

Гораздо сложнее учесть шумы на входе тем более если они проходят нелинейное преобразование модели.

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

Вернемся вновь к СЛАУ Aх=b с квадратной матрицей А размера MхN , которая, в отличие от рассмотренного выше "хорошего" случая (см. разд. 8.Г), требует специального подхода. Обратим внимание на два похожих типа СЛАУ:

• вырожденная система (с нулевым определителем |А|=0 );
• плохо обусловленная система (определитель А не равен нулю, но число обусловленности очень велико).

Несмотря на то, что эти типы систем уравнений существенно отличаются друг от друга (для первого решение отсутствует, а для второго существует и единственно), с практической точки зрения вычислителя, между ними много общего.

Вырожденные СЛАУ

Вырожденная система - это система, описываемая матрицей с нулевым определителем |A|=0 (сингулярной матрицей). Поскольку некоторые уравнения, входящие в такую систему, представляются линейной комбинацией других уравнений, то фактически сама система является недоопределенной. Несложно сообразить, что, в зависимости от конкретного вида вектора правой части ь, существует либо бесконечное множество решений, либо не существует ни одного. Первый из вариантов сводится к построению нормального псевдорешения (т. е. выбора из бесконечного множества решений такого, которое наиболее близко к определенному, например, нулевому, вектору). Данный случай был подробно рассмотрен в разд. 8.2.2 (см. листинги 8.11-8.13).

Рис. 8.7 . Графическое представление несовместной системы двух уравнений с сингулярной матрицей

Рассмотрим второй случай, когда СЛАУ Aх=b с сингулярной квадратной матрицей А не имеет ни одного решения. Пример такой задачи (для системы двух уравнений) проиллюстрирован рис. 8.7, в верхней части которого вводятся матрица А и вектор b , а также предпринимается (неудачная, поскольку матрица А сингулярная) попытка решить систему при помощи функции isolve . График, занимающий основную часть рисунка, показывает, что два уравнения, задающие систему, определяют на плоскости (x0,x1) две параллельные прямые. Прямые не пересекаются ни в одной точке координатной плоскости, и решения системы, соответственно, не существует.

Примечание
Во-первых, заметим, что СЛАУ, заданная несингулярной квадратной матрицей размера 2x2, определяет на плоскости пару пересекающихся прямых (см. рис. 8.9 ниже). Во-вторых, стоит сказать, что, если бы система была совместной, то геометрическим изображением уравнений были бы две совпадающие прямые, описывающие бесконечное число решений .

Рис. 8.8 . График сечений функции невязки f (х) = |Ах-b|

Несложно догадаться, что в рассматриваемом сингулярном случае псевдорешений системы, минимизирующих невязку |Ax-b| , будет бесконечно много, и лежать они будут на третьей прямой, параллельной двум показанным на рис. 8.7 и расположенной в середине между ними. Это иллюстрирует рис. 8.8, на котором представлено несколько сечений функции f(x)= | Аx-b | , которые говорят о наличии семейства минимумов одинаковой глубины. Если попытаться использовать для их отыскания встроенную функцию Minimize , ее численный метод будет всегда отыскивать какую-либо одну точку упомянутой прямой (в зависимости от начальных условий). Поэтому для определения единственного решения следует выбрать из всего множества псевдорешений то, которое обладает наименьшей нормой. Можно пытаться оформить данную многомерную задачу минимизации в Mathcad при помощи комбинаций встроенных функций Minimize , однако более эффективным способом будет использование регуляризации (см. ниже) или ортогональных матричных разложений (см. разд. 8.3).

8.2.3. Вырожденные и плохо обусловленные системы

Вернемся вновь к СЛАУ Aх=b с квадратной матрицей А размера MхN , которая, в отличие от рассмотренного выше "хорошего" случая (см. разд. 8. Г), требует специального подхода. Обратим внимание на два похожих типа СЛАУ:

• вырожденная система (с нулевым определителем |А|=0 );
• плохо обусловленная система (определитель А не равен нулю, но число обусловленности очень велико).

Несмотря на то, что эти типы систем уравнений существенно отличаются друг от друга (для первого решение отсутствует, а для второго существует и единственно), с практической точки зрения вычислителя, между ними много общего.

Вырожденные СЛАУ

Вырожденная система - это система, описываемая матрицей с нулевым определителем |A|=0 (сингулярной матрицей). Поскольку некоторые уравнения, входящие в такую систему, представляются линейной комбинацией других уравнений, то фактически сама система является недоопределенной. Несложно сообразить, что, в зависимости от конкретного вида вектора правой части ь, существует либо бесконечное множество решений, либо не существует ни одного. Первый из вариантов сводится к построению нормального псевдорешения (т. е. выбора из бесконечного множества решений такого, которое наиболее близко к определенному, например, нулевому, вектору). Данный случай был подробно рассмотрен в разд. 8.2.2 (см. листинги 8.11-8.13).

Рис. 8.7. Графическое представление несовместной системы двух уравнений с сингулярной матрицей

Рассмотрим второй случай, когда СЛАУ Aх=b с сингулярной квадратной матрицей А не имеет ни одного решения. Пример такой задачи (для системы двух уравнений) проиллюстрирован рис. 8.7, в верхней части которого вводятся матрица А и вектор b , а также предпринимается (неудачная, поскольку матрица А сингулярная) попытка решить систему при помощи функции isolve . График, занимающий основную часть рисунка, показывает, что два уравнения, задающие систему, определяют на плоскости (x0,xi) две параллельные прямые. Прямые не пересекаются ни в одной точке координатной плоскости, и решения системы, соответственно, не существует.

ПРИМЕЧАНИЕ

Во-первых, заметим, что СЛАУ, заданная несингулярной квадратной матрицей размера 2x2, определяет на плоскости пару пересекающихся прямых (см. рис. 8.9 ниже). Во-вторых, стоит сказать, что, если бы система была совместной, то геометрическим изображением уравнений были бы две совпадающие прямые, описывающие бесконечное число решений.

Рис. 8.8. График сечений функции невязки f (х) = |Ах-b|

Несложно догадаться, что в рассматриваемом сингулярном случае псевдорешений системы, минимизирующих невязку |Ax-b| , будет бесконечно много, и лежать они будут на третьей прямой, параллельной двум показанным на рис. 8.7 и расположенной в середине между ними. Это иллюстрирует рис. 8.8, на котором представлено несколько сечений функции f(x)= = | Аx-b | , которые говорят о наличии семейства минимумов одинаковой глубины. Если попытаться использовать для их отыскания встроенную функцию Minimize , ее численный метод будет всегда отыскивать какую-либо одну точку упомянутой прямой (в зависимости от начальных условий). Поэтому для определения единственного решения следует выбрать из всего множества псевдорешений то, которое обладает наименьшей нормой. Можно пытаться оформить данную многомерную задачу минимизации в Mathcad при помощи комбинаций встроенных функций Minimize , однако более эффективным способом будет использование регуляризации (см. ниже) или ортогональных матричных разложений (см. разд. 8.3).

Плохо обусловленные системы

Плохо обусловленная система - это система, у которой определитель А не равен нулю, но число обусловленности |А -1 | |А| очень велико. Несмотря на то, что плохо обусловленные системы имеют единственное решение, на практике искать это решение чаще всего не имеет смысла. Рассмотрим свойства плохо обусловленных СЛАУ на двух конкретных примерах (листинг 8.14).

Листинг 8.14 . Решение двух близких плохо обусловленных СЛАУ

Каждая строка листинга 8.14 содержит решение двух очень близких плохо обусловленных СЛАУ (с одинаковой правой частью ь и мало отличающимися матрицами А). Несмотря на близость, точные решения этих систем оказываются очень далекими друг от друга. Надо заметить, что для системы двух уравнений точное решение получить легко, однако при решении СЛАУ большой размерности незначительные ошибки округления, неминуемо накапливаемые при расчетах (в том числе и "точным" алгоритмом Гаусса), приводят к огромным погрешностям результата. Возникает вопрос: имеет ли смысл искать численное решение, если заранее известно, что, в силу неустойчивости самой задачи, оно может оказаться совершенно неправильным?

Еще одно соображение, которое вынуждает искать специальные методы решения плохо обусловленных СЛАУ (даже приведенной в качестве примера в листинге 8.14 системы двух уравнений), связано с их физической интерпретацией как результатов эксперимента. Если изначально известно, что входные данные получены с некоторой погрешностью, то решать плохо обусловленные системы не имеет вовсе никакого смысла, поскольку малые ошибки модели (матрицы А и вектора ь) приводят к большим ошибкам решения. Задачи, обладающие такими свойствами, называются некорректными.

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

Рис. 8.9. График хорошо обусловленной системы двух уравнений

Рис. 8.10. График плохо обусловленной системы двух уравнений

Из рис. 8.10 видно, что прямые, соответствующие плохо обусловленной СЛАУ, располагаются в непосредственной близости друг от друга (почти параллельны). В связи с этим малые ошибки в расположении каждой из прямых могут привести к значительным погрешностям локализации точки их пересечения (решения СЛАУ) в противоположность случаю хорошо обусловленной системы, когда малые погрешности в наклоне прямых мало повлияют на место точки их пересечения (рис. 8.9).

ПРИМЕЧАНИЕ

Плохая обусловленность матрицы типична и при реконструкции экспериментальных данных, задаваемых переопределенными (несовместными) СЛАУ (например, в задачах томографии). Именно такой случай приведен в качестве примера в следующем разделе (см. листинг 8.16 ниже).

Метод регуляризации

Для решения некорректных задач, в частности, вырожденных и плохо обусловленных СЛАУ, разработан очень эффективный прием, называемый регуляризацией. В его основе лежит учет дополнительной априорной информации о структуре решения (векторе априорной оценки хо), которая очень часто присутствует в практических случаях. В связи с тем, что регуляризация была подробно рассмотрена в разд. 6.3.3, напомним лишь, что задачу решения СЛАУ Аx=b можно заменить задачей отыскания минимума функционала Тихонова:

Ω (х,λ ) = |Ах-b| 2 +λ |х-х0| 2 . (8.3)

Здесь Я, - малый положительный параметр регуляризации, который необходимо подобрать некоторым оптимальным способом. Можно показать, что задачу минимизации функционала Тихонова можно, в свою очередь, свести к решению другой СЛАУ:

(А T А+λ I)-х=А T В+λ х0, (8.4)

Которая при λ ->0 переходит в исходную плохо обусловленную систему, а при больших x , будучи хорошо обусловленной, имеет решение х 0 . Очевидно, оптимальным будет некоторое промежуточное значение А , устанавливающее определенный компромисс между приемлемой обусловленностью и близостью к исходной задаче. Отметим, что регуляризационный подход сводит некорректную задачу к условно-корректной (по Тихонову) задаче отыскания решения системы (8.4), которое, в силу линейности задачи, является единственным и устойчивым.

Приведем, без излишних комментариев, регуляризованное решение вырожденной системы, которая была представлена на рис. 8.8. Листинг 8.15 демонстрирует отыскание решения задачи (8.4), а полученная зависимость невязки и самого решения от параметра регуляризации Я показана на рис. 8.11 и 8.12 соответственно. Важно подчеркнуть, что решения исходной системы и, следовательно, системы (8.4) при λ =0 не существует.

Листинг 8.15 .Регуляризация вырожденной СЛАУ

Заключительным шагом регуляризации является выбор оптимального λ . Имеется, по крайней мере, два соображения, исходя из которых, можно выбрать параметр регуляризации, если опираться на зависимость от него невязки. В рассматриваемом примере применим критерий определения λ , опирающийся -на подбор нормы невязки, равной априорной оценке погрешностей задания входных данных: матрицы А и вектора b , т. е. величине |δA | + |5λ| . Например, можно выбрать норму невязки и, соответственно, параметр λ и решение х(λ ), которые отмечены на рис. 8.11 и 8.12 пунктирами.

ПРИМЕЧАНИЕ 1

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

ПРИМЕЧАНИЕ 2

Полезно убедиться в том, что формула (8.4) в случае линейной задачи дает тот же результат, что и общая формула (8.3). Для этого достаточно изменить в листинге 8.15 последнюю строку, выражающую решение СЛАУ (8.4), на код, реализующий минимизацию численным методом, как это показано в листинге 8.16. Расчеты (которые требуют значительно большего компьютерного времени) дают тот же самый результат, что был приведен на рис. 8.11 и 8.12.

ПРИМЕЧАНИЕ 3

Попробуйте в расчетах листингов 8.15 и 8.16 взять иную, например, более реалистичную, априорную оценку решения (вместо использованного в них нулевого вектора х0) и посмотреть, как это повлияет на результат.

Рис. 8.11. Зависимость невязки регупяризованного решения вырожденной СЛАУ от параметра А. (продолжение листинга 8.15)

ПРИМЕЧАНИЕ 4

Любопытно также применить вместо формулы (8.3) в качестве функционала Тихонова другую зависимость: Ω (х, λ ) = |Ах-b|+ λ |х-х0 | . Соответствующий пример расчетов вы найдете на компакт-диске.

Рис. 8.12. Регуляризованное решение в зависимости от λ (продолжение листинга 8.15)

Листинг 8.16. Регуляризация СЛАУ при помощи алгоритма минимизации (продолжение листинга 8.15)