Если у вас возникли сложности с курсовой, контрольной, дипломной, рефератом, отчетом по практике, научно-исследовательской и любой другой работой - мы готовы помочь.
Предоплата всего
Подписываем
8
ЭКСПЕРИМЕНТАЛЬНОЕ ИССЛЕДОВАНИЕ ЧИСЛЕННЫХ МЕТОДОВ РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ
1. Цель работы.
Изучение прямых и итерационных методов решения систем алгебраических уравнений с позиций точности получаемых решений и сходимости итерационных процессов.
2. Краткие сведения из теории
Нахождение решения системы алгебраических уравнений
(1)
где:,
является одной из основных задач линейной алгебры.
Хорошо известно, что теоретическим условием существования решения задачи (1) является отличие от нуля определителя матрицы А. Однако, при реализации любого численного алгоритма выполнение этого условия еще не гарантирует нахождение достоверного решения. В общем случае, численно полученное решение отличается от точного решения . Источниками ошибки численного решения являются:
Проблема оценки величины ошибки численного решения
является едва ли не центральной проблемой реализации любого численного метода. Ниже рассматриваются механизмы, которые, будучи включенными в программу решения системы , позволят судить о приемлемости полученного результата.
2.1. Обусловленность вычислительной задачи и простейшие формулы оценки ошибки
Все методы оценки величины ошибки решения системы (1) используют значение числа обусловленности матрицы системы. Последнее определяется разными способами, но, по сути, всегда является коэффициентом пропорциональности между вариациями исходных данных (значения элементов матрицы и вектор правой части) и соответствующей вариацией решения.
Вообще говоря, в силу наличия ошибок округления и возможных погрешностей в исходных данных всегда решается возмущенная задача
(2)
и полученное численное решение (вектор ) удовлетворяет не уравнению (1), а уравнению
(3)
где - вектор невязки. Учтя, что , получим :
и
(4)
Используя свойства векторной и матричной норм, после простых преобразований получим:
(5)
Полученная оценка погрешности мало пострадает, если в правой части выражения (5) заменить на . Тогда
. (6)
Здесь
- естественное число обусловленности.
Очевидно, что такое число обусловленности может быть вычислено только после решения задачи. Вероятно в силу этого, оно и не получило широкого применения.
Зависимость числа обусловленности от нормы вектора решения исключается в стандартном числе обусловленности . Применение его позволяет получить следующее выражение для оценки погрешности решения
. (7)
Отметим, что оба выражения (6) и (7) являются лишь оценками возможной ошибки решения. Они почти всегда дают завышенную оценку, иногда - сильно завышенную, и редко - заниженную.
Трудность использования каждой из этих оценок заключается в необходимость вычисления нормы обратной матрицы || А-1||. Обращение матрицы А и последующее вычисление ее нормы является слишком «дорогостоящим мероприятием» для вычисления оценки решения. Для практического применения выражения (7) вполне достаточно было бы располагать лишь оценкой величины А-1. Рассматривая эту проблему, будем для определенности считать, что система решена методом исключения Гаусса и известны:
Решение системы (1) другими прямыми методами, например посредством QR-алгоритма, не вносит никаких принципиальных изменений в приведенные ниже соображения.
Итак, имеются три альтернативы.
, (8)
где - решение задачи , и - псевдослучайные векторы. При этом оказывается возможным достаточно достоверно оценить величину при сравнительно небольшом количестве () решений задачи (8). Этот алгоритм потребует примерно дополнительных операций, что составляет лишь малую часть общих трудозатрат на решение задачи.
В численном анализе получили распространение две формы учета возмущения исходной задачи (1).
Первая:
, (9)
где P и D - матрицы возмущения.
Поскольку возмущения исходных данных имеют, как правило, случайный характер, то компоненты матриц возмущения можно представить в виде:
, - случайное число
и , - случайное число ,;
- максимальные относительные ошибки в задании элементов матрицы А и вектора .
Вторая:
, (10)
где: и di,j - определяются по приведенному выше выражению.
Ясно, что модель (9) приводится к виду (10), если т.е. при . Но отсюда следует, что при "плохих" матрицах А ( при "больших" ) результаты исследования чувствительности решения системы по (9) и (10) могут существенно различаться при одинаковых значениях .
Задавая в эксперименте различные значения , где - машинный эпсилон, а и , можно оценить приемлемую, с позиций требуемой точности решения, неопределенность в задании исходных данных. Это позволяет дать разумное обоснование требуемой точности измерения (вычисления) компонент исходных данных задачи (1). Одновременно, зная реальныее величины и , можно оценить достоверность полученного решения. Последнему служит, так называемая, составная формула оценки ошибки.
Исходя из (9), для найденного решения можем записать
А (Е + Р) * = (Е + D) + , (11)
в предположении относительной малости компонент матриц P, D и вектора можно получить формулу оценки ошибки представленную ниже:
x || Р ||+ + (12)
C учетом соображений о способах оценки величины нормы || А||формула (12) дает достаточно экономичный способ вычисления оценки ошибки. Действительно,
А * < ||А-1||* ;
Т.к. по (8а)
|| А-1 || 2 ,
то, взяв в качестве вектор невязки , получим :
А 2 * |||| и
2 * .
Аналогично оценивается величина слагаемого .
При выборе в (8а) = , где: - константа , < 1 , i = 1, 2,…, k , можно получить
1,3 * .
Коэффициент 1,3 взят для повышения надежности оценки при малых значениях k.
Норма матрицы Р вычисляется на основе анализа исходных данных - либо известна их неопределенность, либо предполагается, что pпорядка машинного эпсилона.
Обобщением приведенных выше соображений является следующий алгоритм включения составной формулы оценки ошибки в процессе решения системы линейных уравнений методом исключения Гаусса.
Алгоритм 1.
НАЗНАЧЕНИЕ :
Решение системы А = методом Гаусса с составной формулой оценки ошибки.
0. НАЧАЛО алгоритма.
3.1. Вычислить невязку ;
3.2. Вычислить из систем: L = и U = ;
4. Для k = 1,2 выполнить :
4.1. сгенерировать случайное число
4.2 вычислить из системы :
L , U = ;
5. ERR1 = 2 * / ;
6. ERR2 = 1,3 * ;
7. ERR3 = n *;
9. КОНЕЦ алгоритма.
В случае, если возмущенная задача описывается выражением (10), составная формула оценки модифицируется и имеет вид:
x + + . (13)
Эта формула несколько более трудоемка в сравнении с (12), так как дополнительно требуется оценить и . Это можно сделать посредством решения системы А для нескольких случайных векторов .
Отметим, что данные выше оценки решения допускают следующую интерпретацию:
«ошибка решения = неопределенность модели + неопределенность исходных данных + ошибки вычислений»
Здесь матрица задачи интерпретируется как модель исследуемой системы, а правая часть - как вектор возбуждающих воздействий (исходные данные).
Разумное использование этих оценок требует предварительных суждений об относительных размерах неопределенностей модели и исходных данных. Во многих реальных ситуациях ошибки вычислений пренебрежимо малы в сравнении с этими неопределенностями, а невероятные решения обычно являются только отражением высокой чувствительности модели (плохой обусловленности матрицы) и не свидетельствуют о плохо проведенных вычислениях (низком качестве алгоритма).
Интересно, что существуют ситуации, когда преднамеренное возмущение модели является хорошим способом получить правдоподобное решение. Если обусловленность матрицы задачи > , то надежда получить приемлемое решение крайне мала. Столь плохая обусловленность задачи говорит о значительной степени линейной зависимости строк (столбцов) матрицы системы. Внесение в нее малых возмущений порядка часто позволяет существенно улучшить обусловленность системы и, как правило, не приводит к заметному изменению исходной задачи. Практически, можно рекомендовать генерацию матрицы возмущения вида
. (14)
Задача (1) тем самым преобразуется в (А + Md ) = .
2.3. И т е р а ц и о н н ы е м е т о д ы .
Итерационные методы решения задачи (1) являются теоретически бесконечно-шаговыми и приводят лишь к приближенным решениям, что ставит под сомнение целесообразность их использования. Однако существуют ситуации, когда применение прямых методов становиться неэффективным. Так, для задач с сильно разреженными матрицами, (они возникают при решении краевых задач для дифференциальных уравнений в частных производных, а также при анализе больших электронных цепей), итерационные алгоритмы требуют гораздо меньшей оперативной памяти, чем прямые методы. Кроме этого, для задач большой размерности (1000 и более) трудоемкость прямых алгоритмов становиться слишком большой. Эти обстоятельства делают оправданным использование итерационных алгоритмов.
Каноническая форма записи любого стационарного итерационного метода имеет вид
(15)
Методы Якоби и Гаусса - Зейделя являются наиболее простыми из многочисленного семейства этих алгоритмов. Однако их свойства и особенности применения достаточно типичны.
Главной отличительной чертой всех итерационных процедур является сильная зависимость самой возможности получения решения от свойств итерационной матрицы B, которая интегрирует свойства матрицы задачи и свойства метода. Эта связь формулируется теоремой о сходимости метода, которая для всякого стационарного итерационного метода требует, чтобы все собственные числа матрицы В по модулю были меньше 1, т.е.
. < 1 (16)
Для любой матрицы справедлива следующая оценка сверху
< .
Отсюда вытекает чисто качественное свойство желательности "малости" матрицы В.
Иногда пытаются установить чисто качественные, «визуальные» признаки сходимости того, или иного алгоритма. Так, для метода Якоби матрица В и вектор имеют вид :
B = D * ( A - D ) , = D * b , (17)
где D – диагональ матрицы А. "Малость" матрицы В в этом случае, как видно из (17), предполагает хорошую аппроксимацию матрицы системы диагональной матрицей.
Для метода Гаусса-Зейделя:
B = E - (D + L ) * A , = (D + L ) * , (18)
где L – нижний треугольник матрицы. В этом случае "малость" В реализуется при близости матрицы А к нижне-треугольной матрице (D + L). Отсюда следует, что можно ожидать хорошую сходимость процедуры Гаусса-Зейделя для "почти" нижне-треугольных матриц. Однако такие признаки носят исключительно качественный характер и не являются условиями сходимости. Так, например, процедура Гаусса-Зейделя при решении задачи (1) с матрицей
А =
не сходится, несмотря на ее "почти" нижне-треугольный вид. Такой результат является совершенно закономерным, ибо собственные числа соответствующей матрицы В следующие:
= (-2,2853 , -0,6845 , 0,0)
В заключение отметим, что итерационные методы для плохо обусловленных задач нисколько не лучше прямых.
Исследование проводится с использованием программы LinEqv, которая позволяет:
Необходимо провести численные эксперименты, на основании которых:
5. Содержание отчета.
Литература.
ПРИЛОЖЕНИЕ
Procedure TextMatr (n,k:integer; var a:matr);
{НАЗНАЧЕНИЕ – формировать тестовые матрицы}
{ВХОДНЫЕ ПЕРЕМЕННЫЕ
n – размерность системы
k – порядковый номер матрицы}
{ВЫХОДНЫЕ ПЕРЕМЕННЫЕ
a – двумерный массив в первых n строках и n столбцах которого записана сформированная матрица}
Var
i,j:integer;
alf,xmin,xmax,x:real;
Begin
CASE k OF
1: {хорошо обусловленная матрица}
FOR i:=1 TO n DO
FOR j:=1 TO n DO
IF j=i THEN a[i,j]:=i ELSE
IF j=i+1 THEN a[i,j]:=n ELSE a[i,j]:=0.0;
2: хорошо обусловленная матрица }
FOR i:=1 TO n DO
FOR j:=1 TO n DO
IF j < i THEN a[i,j]:=j ELSE a[i,j]:=i;
3:{ матрица с хорошей и практически постоянной обусловленностью }
FOR i:=1 TO n DO
FOR j:=1 TO n DO
begin
IF i = j THEN a[i,i]:=0.01/((n-i+1)*(i+1));
IF i > j THEN a[i,j]:=i*(n-j);
IF i < j THEN a[i,j]:=j*(n-i);
end;
4: { псевдо-случайная матрица с хорошей и практически постоянной обусловленностью }
begin
RANDOMIZE;
FOR i:=1 TO n DO
FOR j:=1 TO n DO
a[i,j]:=2.0 * (RANDOM - 0.5);
end;
5: { Матрица Гильберта обусловленность быстро растет с ростом порядка}
FOR i:=1 TO n DO
FOR j:=1 TO n DO a[i,j]:=1.0/(i+j-1);
6: { плохо обусловленная матрица }
FOR i:=1 TO n DO
FOR j:=1 TO n DO
begin
IF i = j THEN a[i,i]:=1.0/((n-i+1)*(i+1));
IF i > j THEN a[i,j]:=i*(n-j);
IF i < j THEN a[i,j]:=0.0;
end;
7: {матрица с варьируемой обусловленностью}
begin
WriteLn (' What is a l f a ? ( alfa < 0)');
WriteLn (' Hard case is |alfa| very small or very big');
ReadLn (alf);
IF alf > 1.0E-3 THEN alf:=-alf;
FOR i:=1 TO n DO
FOR j:=1 TO n DO a[i,j]:=EXP(alf*i*j);
end;
8: begin
{матрица с варьируемой обусловленностью}
WriteLn (' What is a l f a ? ');
WriteLn (' Hard case is alfa > 1');
ReadLn (alf);
x:=ln(2.0);
RANDOMIZE;
FOR i:=1 TO n DO
FOR j:=1 TO n DO a[i,j]:=alf + ln(i*j)/x +
1.0E-05*(RANDOM - 0.5);
end;
9: { Pei - matrix }
begin
WriteLn (' What is a l f a ? ');
WriteLn (' Hard case is : alfa=(1.0 + - epsilon) ');
Readln (alf);
FOR i:=1 TO n DO a[i,i]:= alf;
FOR i:=1 TO n DO
FOR j:=1 TO n DO
IF j <>i THEN a[i,j]:=1.0;
end;
10: { матрица задачи полиномиального интерполирования }
begin
Writeln ( ' ':10,'Matrix for polinomial interpolation');
WriteLn (' What are the bounds of interval for interpolation');
ReadLn (xmin,xmax);
FOR i:=1 TO n DO a[i,1]:=1.0;
alf:=(xmax - xmin)/(n-1);
FOR i:=1 TO n DO
begin
x:=xmin + (i-1)*alf;
FOR j:=2 TO n DO a[i,j]:=a[i,j-1]*x;
end;
end;
11: { плохо обусловленная матрица }
FOR i:=1 TO n DO
FOR j:=1 TO n DO a[i,j]:=Sin(i*j/(2*n));
12: { плохо обусловленная матрица }
FOR i:=1 TO n DO
FOR j:=1 TO n DO
begin
alf:=i*j/n;
a[i,j]:=1.0/(1.0 + 66.0 * alf*alf*alf*alf);
end;
13: {матрица 11 с малой случайной добавкой }
begin
RANDOMIZE;
FOR i:=1 TO n DO
FOR j:=1 TO n DO
a[i,j]:=Sin(i*j/(2*n)) + 1.0E-06 * (RANDOM - 0.5);
end;
end { case }
end;