Цель работы. Изучение прямых и итерационных методов решения систем алгебраических уравнений с позиций то
Работа добавлена на сайт samzan.net:
ЭКСПЕРИМЕНТАЛЬНОЕ ИССЛЕДОВАНИЕ ЧИСЛЕННЫХ МЕТОДОВ РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ
1. Цель работы.
Изучение прямых и итерационных методов решения систем алгебраических уравнений с позиций точности получаемых решений и сходимости итерационных процессов.
2. Краткие сведения из теории
Нахождение решения системы алгебраических уравнений
(1)
где:,
является одной из основных задач линейной алгебры.
Хорошо известно, что теоретическим условием существования решения задачи (1) является отличие от нуля определителя матрицы А. Однако, при реализации любого численного алгоритма выполнение этого условия еще не гарантирует нахождение достоверного решения. В общем случае, численно полученное решение отличается от точного решения . Источниками ошибки численного решения являются:
неопределенность в задании исходных данных (матрицы А и вектора )
ошибки представления чисел в ЭВМ
дополнительные возмущения, возникающие при реализации вычислительного алгоритма посредством арифметики конечной разрядности (устойчивость численного алгоритма).
Проблема оценки величины ошибки численного решения
является едва ли не центральной проблемой реализации любого численного метода. Ниже рассматриваются механизмы, которые, будучи включенными в программу решения системы , позволят судить о приемлемости полученного результата.
2.1. Обусловленность вычислительной задачи и простейшие формулы оценки ошибки
Все методы оценки величины ошибки решения системы (1) используют значение числа обусловленности матрицы системы. Последнее определяется разными способами, но, по сути, всегда является коэффициентом пропорциональности между вариациями исходных данных (значения элементов матрицы и вектор правой части) и соответствующей вариацией решения.
Вообще говоря, в силу наличия ошибок округления и возможных погрешностей в исходных данных всегда решается возмущенная задача
(2)
и полученное численное решение (вектор ) удовлетворяет не уравнению (1), а уравнению
(3)
где - вектор невязки. Учтя, что , получим :
и
(4)
Используя свойства векторной и матричной норм, после простых преобразований получим:
(5)
Полученная оценка погрешности мало пострадает, если в правой части выражения (5) заменить на . Тогда
. (6)
Здесь
- естественное число обусловленности.
Очевидно, что такое число обусловленности может быть вычислено только после решения задачи. Вероятно в силу этого, оно и не получило широкого применения.
Зависимость числа обусловленности от нормы вектора решения исключается в стандартном числе обусловленности . Применение его позволяет получить следующее выражение для оценки погрешности решения
. (7)
Отметим, что оба выражения (6) и (7) являются лишь оценками возможной ошибки решения. Они почти всегда дают завышенную оценку, иногда - сильно завышенную, и редко - заниженную.
Трудность использования каждой из этих оценок заключается в необходимость вычисления нормы обратной матрицы || А-1||. Обращение матрицы А и последующее вычисление ее нормы является слишком «дорогостоящим мероприятием» для вычисления оценки решения. Для практического применения выражения (7) вполне достаточно было бы располагать лишь оценкой величины А-1. Рассматривая эту проблему, будем для определенности считать, что система решена методом исключения Гаусса и известны:
вектор численного решения
вектор невязки
LU разложение матрицы А.
Решение системы (1) другими прямыми методами, например посредством QR-алгоритма, не вносит никаких принципиальных изменений в приведенные ниже соображения.
Итак, имеются три альтернативы.
1. Вычислить матрицу и ее норму. При наличии LU разложения матрицы А решение задачи обращения последней потребует около дополнительных операций и примерно в 4 раза увеличит общую трудоемкость решения системы (1).
2. Грубо оценить . Такую оценку можно получить на основании выражения
, (8)
где - решение задачи , и - псевдослучайные векторы. При этом оказывается возможным достаточно достоверно оценить величину при сравнительно небольшом количестве () решений задачи (8). Этот алгоритм потребует примерно дополнительных операций, что составляет лишь малую часть общих трудозатрат на решение задачи.
3. Использовать оценку стандартного числа обусловленности из хорошо зарекомендовавшей себя процедуры DECOMP (пакет Форсайта; ее паскальная версия в программном комплексе LinEqv представлена процедурой PDECOMP). В этой процедуре идея предыдущего алгоритма усовершенствована так, что единственный вектор ищется из системы , где вектор с компонентами 1 строится по принципу максимизации модулей компонентов вектора решения.
2.2 Анализ чувствительности решения к возмущениям исходных данных и составная
формула оценки ошибки
В численном анализе получили распространение две формы учета возмущения исходной задачи (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. НАЧАЛО алгоритма.
1. Выполнить LU-разложение матрицы системы исключением Гаусса с частичным выбором ведущего элемента;
2. Вычислить ;
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 *;
ERR = ERR1 + ERR2 + ERR3
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)
В заключение отметим, что итерационные методы для плохо обусловленных задач нисколько не лучше прямых.
4. Задание на работу.
Исследование проводится с использованием программы LinEqv, которая позволяет:
сгенерировать или ввести с пульта матрицу системы
вычислить по заданному точному решению вектор правой части системы ()
получить решение системы А = различными методами
вычислить фактически полученную ошибку в соответствии с выражением (2)
вычислить обратную матрицу А и оценить качество обращения по величине
вычислить числа обусловленности и и оценку стандартного числа обусловленности .
вычислить вектор невязки и его норму
вычислить оценки величины ошибки по (7) , (12) и (13) ;
внести возмущение в матрицу системы, соответствующее моделям возмущения типа P, M, Md (выражения (9), (10), (14)).
Необходимо провести численные эксперименты, на основании которых:
1 Сравнить между собой естественное и стандартное числа обусловленности матрицы а также - точное значение стандартного числа обусловленности с его оценкой, вычисленной процедурой DECOMP.
2 Оценить точность решений, получаемых методом исключения Гаусса для систем с одинаково хорошо обусловленными матрицами порядка от 3 до 15; провести анализ точности, как функции порядка матрицы; сравнить фактически получаемую ошибку с ее оценками.
3 Выполнить задания п.2 для систем с одинаково плохо обусловленными матрицами.
4 Оценить точность решений, получаемых методом исключения Гаусса для систем одного порядка, но различной обусловленности (от «очень хороших» до «очень плохих»). Результаты анализа представить в виде зависимости относительной точности решения от числа обусловленности. Обратить внимание на величину нормы вектора невязки и проследить ее зависимость от обусловленности системы и связь с фактической ошибкой решения.
5 Исследовать возможность улучшения обусловленности задачи посредством внесения малого случайного возмущения в матрицу системы.
6 Для 2-3 задач с «хорошей» матрицей ( = 102 - 104 ) посредством внесения в матрицу системы возмущений различной величины сделать заключение о приемлемой для получения требуемой (наперед заданной) точности решения степени неопределенности в задании исходных данных.
7 Повторить эксперимент п.6 для 2-3 задач с плохо обусловленной матрицей.
8 Выполняя п.п. 6 и 7 , исследовать работоспособность различных методов оценки ошибок решения ( выражения (7), (12), (13) ) при наличии возмущения левой части системы.
9 Применить для решения нескольких систем из пунктов 2-4 итерационные методы Якоби и Гаусса-Зейделя; проверить реализацию задаваемого критерия точности. Исследованием спектра матрицы В проверить выполнение теоремы сходимости стационарного метода; выявить взаимосвязь скорости сходимости итерационного процесса с величиной спектрального радиуса матрицы В.
10 Провести исследование влияния вида доминирования матрицы задачи на сходимость процедур Якоби и Гаусса-Зейделя.
5. Содержание отчета.
Каждый пункт отчета должен содержать
формулировку цели эксперимента (что исследуется)
условия проведения эксперимента (что не изменяется, что будет изменяться и т.п.)
полученные результаты в форме таблицы и, если возможно, диаграмм (не протокол!!!),
выводы и объяснения результатов.
Протокол вычислительных экспериментов прикладывается к отчету.
Литература.
А.А. Самарский, А.В. Гулин Численные методы. - Наука, 1989.
В.В. Воеводин. Вычислительные основы линейной алгебры. - М.: Наука, 1977
ПРИЛОЖЕНИЕ
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;
2. Курсовая работа- Полномочия органов местного самоуправления.html
3. Реферат- Проблемы урегулирования задолженности предприятий-банкротов
4. Вариант1 ОСД 1
5. постоянных затрат в общем объеме 40
6. неолитическая революция как фактор социального расслоения общества появления классов собственности го
7. ___ ______________ 2009р. Р О Б О Ч А П Р О Г Р А М А Вид і назва практики- Виробнича Комплексна п
8. Оглавление Введение Структурная схема Функциональная схема Принцип действия Основные
9. Дворец детского творчества Из серии Сладкое творчество ПОДСОЛНУХ
10. 1 - 316 - 9301 ЄСИПЕНКО Дмитро Миколайович СЕНС ІСТОРІЇ ЯК СВІТОГЛЯДНИЙ ОРІЄНТИР ОСОБИ
11. лекция для педагогов школы Янушевская Наталья Анатольевна учитель физики высшей категории МКОУ
12. инновационному развитию Республики Казахстан на 2010 ' 2014 годы Астана 2010
13. Его получают из бескарбонатных глин обжигом при температуре 1200 С
14. экономических дисциплин ИНФОРМАЦИОННАЯ ЭВРИСТИКА Учебнометодический комплекс Для
15. Time jobs for schoolchildren is delivering newsppers to people~s doors
16. окклюзивная косметика; истощение
17. Адвокат в современном уголовном процессе- Пособие для адвокатов - Под ред
18. Типы сканирующих устройств
19. Диапазон приемлемости
20. тема управления должны обеспечить свободный поток информации снизу вверх Признак высокого качества систем