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

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

Аz=u, (3; 2,1)

где А -- матрица с элементами a ij , А ={a ij }, z -- искомый вектор с координатами z j , z={z j }, и -- известный вектор с координатами и i ,u = {u i }, i, j =1, 2, ..., п. Система (3; 2,1) называется вырожденной, если определитель системы равен нулю, detA = 0. В этом случае матрица А имеет равные нулю собственные значения. У плохо обусловленных систем такого вида матрица А имеет близкие к нулю собственные значения.

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

В практических задачах часто правая часть и и элементы матрицы А, т. е. коэффициенты системы (3; 2,1), известны приближенно. В этих случаях вместо системы (3;2,1) мы имеем дело с некоторой другой системой Az=и такой, что ||A-A||<=h, ||u-u||<=--d, где смысл норм обычно определяется характером задачи. Имея вместо матрицы А матрицу A, мы тем более не можем высказать определенного суждения о вырожденности или невырожденности системы (3; 2,1).

В этих случаях о точной системе Аz=u , решение которой надо определить, нам известно лишь то, что для матрицы А и правой части и выполняются неравенства ||A-A||<=h, ||u-u||<=--d. Но систем с такими исходными данными (А, и) бесконечно много, и в рамках известного нам уровня погрешности они неразличимы. Поскольку вместо точной системы (3; 2,1) мы имеем приближенную систему Аz= и, то речь может идти лишь о нахождении приближенного решения. Но приближенная система Аz=и может быть неразрешимой. Возникает вопрос:

что надо понимать под приближенным решением системы (3; 2,1) в описанной ситуации?

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

Таким образом, в большом числе случаев мы должны рассматривать целый класс неразличимых между собой (в рамках заданного уровня погрешности) систем уравнений, среди которых могут быть и вырожденные, и неразрешимые. Методы построения приближенных решений систем этого класса должны быть одними и теми же, общими. Эти решения должны быть устойчивыми к малым изменениям исходных данных (3; 2,1).

В основе построения таких методов лежит идея «отбора». Отбор можно осуществлять с помощью специальных, заранее задаваемых функционалов W[ z ] , входящих в постановку задачи.

Неотрицательный функционал W[ z ] , определенный на всюду плотном в F подмножестве F 1 множества F, называется стабилизирующим функционалом, если:

• а) элемент z T принадлежит его области определения;
• б) для всякого числа d>0 множество F 1,d элементов z из F 1 , для которых
• W[ z ]

Итак, рассмотрим произвольную систему линейных алгебраических уравнений (короче -- СЛАУ)

Аz =u , (3; 2,2)

в которой z и u--векторы, z=(z 1 , z 2 , ...,z n)--ОR n , и =(u 1 ,u 2 , ... ,u n)--ОR m , А-- матрица с элементами a ij , А = {a ij }, где j =1, 2, ..., n; i= 1, 2, ..., т, и число п не обязано быть равным числу т.

Эта система может быть однозначно разрешимой, вырожденной (и иметь бесконечно много решений) и неразрешимой.

Псевдорешением системы (3; 2,2) называют вектор z, минимизирующий невязку || Az - u || на всем пространстве R n . Система (3; 2,2) может иметь не одно псевдорешение. Пусть F A есть совокупность всех ее псевдорешений и z 1 -- некоторый фиксированный вектор из R n , определяемый обычно постановкой задачи.

Нормальным относительно вектора z 1 решением системы (3;2,2) будем называть псевдорешение z 0 с минимальной нормой || z - z 1 ||, т. е. такое, что

|| z 0 - z 1 || =

Здесь . В дальнейшем для простоты записи будем считать z 1 = 0 и нормальное относительно вектора z 1 =0 решение называть просто нормальным решением.

Для любой системы вида (3; 2,2) нормальное решение существует и единственно.

Замечание 1. Нормальное решение z° системы (3;2,2) можно определить также как псевдорешение, минимизирующее заданную положительно определенную квадратичную форму относительно координат вектора z--z 1 . Все излагаемые ниже результаты остаются при этом справедливыми.

Замечание 2. Пусть ранг матрицы А вырожденной системы (3; 2,1) равен r < n и z r+1 ,z r+2 , … , z n -- базис линейного пространства N A , состоящего из элементов z, для которых Аz=0, N A = {z; Аz= 0}. Решение z° системы (3; 2,1), удовлетворяющее n--r условиям ортогональности

(z 0 - z 1 , z S)= 0, S= r + 1, r + 2, .. ,n , (3; 2,3)

определяется однозначно и совпадает с нормальным решением.

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

z = Vz*, u = Vu *

ее можно привести к диагональному виду и преобразованная система будет иметь вид

l i z i *=u i * , i= 1, 2,. .., п,

где l i -- собственные значения матрицы А.

Если симметричная матрица А -- невырожденная и имеет ранг r, то n - r ее собственных значений равны нулю. Пусть

l i №0 для i=1, 2, ..., r;

l i =0 для i=r+1,r+2, …, n.

Полагаем, что система (3; 2,2) разрешима. При этом u i *= 0 для i =r + 1, ..., n.

Пусть «исходные данные» системы (А и и) заданы с погрешностью, т. е. вместо А и и заданы их приближения А и u:

|| A - A ||<=h, ||u - u||<=d . При этом

Пусть l i -- собственные значения матрицы А. Известно, что они непрерывно зависят от А в норме (3; 2,4). Следовательно, собственные значения l r+1 , l r+2 , …,l n могут быть сколь угодно малыми при достаточно малом h.

Если они не равны нулю, то

Таким образом, найдутся возмущения системы в пределах любой достаточно малой погрешности А и и, для которых некоторые z i * будут принимать любые наперед заданные значения. Это означает, что задача нахождения нормального решения системы (3; 2,2) является неустойчивой.

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

Вернемся вновь к СЛАУ 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).

УДК 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 с.

Лабораторная работа № 3

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

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

Входные параметры: n-целое положительное число, равное порядку n системы; а - массив из n х n действительных чисел, содержащий матрицу коэффициентов системы; b - массив из n действительных чисел, содержащий столбец свободных членов системы (b(1) = b 1 , b(2)=b 2, …b(n)=b n).

Выходные параметры: х – решение системы; p-количество итераций.

Схема алгоритма приведена на рисунке 18.

Текст программы:

procedure regul(N:Integer;a:Tmatr;b:Tvector;var X:Tvector; var p:integer);

var a1,a2:tmatr; b1,b2,x0:tvector; alfa,s1,s:real; max,eps:real; i,j,k,l:integer;

Out_Slau_T(n,a,b);

For I:=1 To n Do {получение А Т А}

For K:=1 To N Do

For J:=1 To N Do S:=S+A*A;

For I:=1 To N Do {получение А Т В}

For J:=1 To N Do

Begin S:=S+A*B[j];

alfa:=0; {нач-ое знач-е альфа}

k:=0; {кол-во итераций}

alfa:=alfa+0.01; inc(k); a2:=a1;

for i:=1 to N do a2:=a1+alfa; {получение А Т А+alfa}

for i:=1 to N do b2[i]:=b1[i]+alfa*x0[i]; {получение А Т В+alfa}

SIMQ(n,a2,b2,l);

a2:=a1; X:=b2; x0:=X; b2:=b1;

vozm(N,eps,a2,b2);

simq(n,a2,b2,l);

for i:=2 to n do

if abs(b2[i]-X[i])>max then max:=abs(b2[i]-X[i]);

X1 = 1.981 X2 = 0.4735

Рисунок 18 - Схема алгоритма метода регуляризации

Варианты заданий для решения плохо обусловленных систем методом регуляризации приведены в таблице 3.

Метод вращения (Гивенса)

Схема алгоритма приведена на рисунке 19.

Пример. Решить систему уравнений

Текст программы:

PROCEDURE Vrash;

Var I,J,K: Integer; M,L,R: Real; F1:TEXT; Label M1,M2;

Out_Slau_T(nn,aa,b);

for i:=1 to Nn do

For I:=1 To Nn-1 Do Begin

For K:=I+1 To Nn Do Begin

If (Aa0.0) Then Goto M1;If (Aa0.0) Then Goto M1;

1:M:=Sqrt(Aa*Aa+Aa*Aa);

L:=-1.0*Aa/M;

M2:For J:=1 To Nn Do Begin

R:=M*Aa-L*Aa;

Aa:=L*Aa+M*Aa;

R:=M*Aa-L*Aa;

Aa:=L*Aa+M*Aa;

For I:=Nn Downto 1 Do Begin

For K:=0 To Nn-I-1 Do Begin M:=M+Aa*Aa; End;

Aa:=(Aa-M)/Aa; End;

for i:=1 to Nn do x[i]:=Aa;End;

Вычисления по программе привели к следующим результатам:

X1 = 1.981 X2 = 0.4735

Рисунок 19 - Схема алгоритма метода Гивенса (вращения)

Варианты заданий

Таблица 3

	Матрица A				Матрица A

Тема лабораторной работы №3 для контроля знаний проиллюстрирована контрольно – обучающей программой.

Лабораторная работа № 4

Решение нелинейных уравнений и систем нелинейных уравнений

Метод простых итераций

Порядок выполнения лабораторной работы:

Найти нулевое приближение решения;

Преобразовать систему f(x) = 0 к виду х = Ф(х);

Проверить условие сходимости метода.

Схема алгоритма приведена на рисунке 20.

Пример. Решить методом простых итераций систему

В качестве нулевого приближения выберем точку х=1, у = 2.2, z = 2. Преобразуем систему к виду

Текст программы:

PROCEDURE Iteraz;

Var I,J,K,J1: Integer;

X2,X3,Eps: Real;

Eps:=0.01; X2:=0.0; K:=1;

For J:=1 To Nn Do Begin

For I:=1 To Nn Do Begin S:=S+Aa*Xx[i]; End;

For J1:=1 To Nn Do Begin Xx:=R; End; X3:=Xx;

For I:=1 To Nn Do Begin If (Xx[i]>=X3) Then X3:=Xx[i]; Еnd;

For I:=1 To Nn Do Begin Xx[i]:=Xx[i]/X3; End;

X1:=X3; U:=Abs(X2-X1); U1:=U/Abs(X1);

If (U1>=Eps) Then X2:=X1;

Until ((K>=50)or(U1

Вычисления по программе привели к следующим результатам:

X(1)= 1.1132 Х(2)= 2.3718 Х(3)= 2.1365

Количество итераций:5

Рисунок 20 - Схема алгоритма метода простых итераций

Метод Ньютона

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

Входные параметры: n - число уравнений системы (совпадает с числом неизвестных), n£10; х-массив из n действительных чисел, содержащий начальное приближение решения; f- имя внешней процедуры f(n, х, у), вычисляющей по заданным значениям х, расположенным в элементах массива х, текущие значения функции f и размещающей их в элементах массива у; g - имя внешней процедуры g(n, x, d), вычисляющей по заданным значениям х из массива х элементы матрицы
,которая размещается в массиве d размерности n x n; eps - значение условия окончания итерационного процесса.

Выходные параметры: х - массив из n действительных чисел (он же входной) содержит при выходе из подпрограммы приближенное значение решения; k-количество итераций.

Искомый вектор

Если , то система (1) называется плохо обусловленной. В этом случае погрешности коэффициентов матрицы и правых частей или погрешности округления при расчетах могут сильно исказить решение.

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

такую, что

Полагаем, что величины и d известны.

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

Псевдорешением системы (1) называется вектор , минимизирующий невязку на всем пространстве .

Пусть х 1 –некоторый фиксированный вектор из , определяемый обычно постановкой задачи.

Нормальным относительно вектора х 1 решением системы (1) называется псевдорешение х 0 с минимальной нормой , то есть

где F-совокупность всех псевдорешений системы (1).

Причем

где ¾компоненты вектора х.

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

Для нахождения приближенного нормального решения системы (1) воспользуемся методом регуляризации.

Согласно указанному методу построим сглаживающий функционал вида

и найдем вектор , минимизирующий на этот функционал. Причем параметр регуляризации a однозначно определен из условия

где .

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

Компоненты вектора являются решениями системы линейных алгебраических уравнений, которая получается из условия минимума функционала (4)

и имеет вид

где Е-единичная матрица,

¾эрмитово сопряженная матрица.

На практике для выбора вектора нужны дополнительные соображения. Если их нет, то полагают =0.

Для =0 систему (7) запишем в виде

где

Найденный вектор будет являться приближенным нормальным решением системы (1).

Остановимся на выборе параметра a. Если a=0, то система (7) переходит в плохо обусловленную систему. Если a велико, то система (7) будет хорошо обусловлена, но регуляризованное решение не будет близким к искомому решению системы (1). Поэтому слишком большое или слишком малое a не пригодны.

Обычно на практике проводят расчеты с рядом значений параметра a. Например,

Для каждого значения a находят элемент , минимизирующий функционал (4). В качестве искомого значения параметра регуляризации берется такое число a, для которого с требуемой точностью выполняется равенство (5) или (6).

III. ЗАДАНИЕ

1. Построить систему линейных алгебраических уравнений, состоящую из трех уравнений с тремя неизвестными, с определителем, величина которого имеет порядок 10 -6 .

2. Построить вторую систему, аналогичную первой, но имеющую другие свободные члены, отличающиеся от свободных членов первой системы на величину 0,00006.

3. Решить построенные системы методом регуляризации (полагая =0 и d=10 -4) и каким-либо другим методом (например, методом Гаусса).

4. Сравнить полученные результаты и сделать выводы о применимости использованных методов.

IV. ОФОРМЛЕНИЕ ОТЧЕТА

В отчете должны быть представлены:

1. Название работы.

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

3. Описание алгоритма (метода) решения.

4. Текст программы с описанием.

5. Результаты работы программы.

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. - М.: Наука, 1979. 286 с.

2. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. - М.: БИНОМ. Лаборатория Знаний, 2007 636с.

Лабораторная работа № 23