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

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

А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)) возмущениям правой части и, основанного на методе регуляризации.

Лабораторная работа № 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-количество итераций.

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

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

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

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

г. Самара

Е-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.

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

Если , то система (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