Ульяновский государственный университет




Государственное образовательное бюджетное учреждение

Высшего профессионального образования

Ульяновский государственный университет

 

Факультет Математики и информационных технологий

Кафедра Информационных технологий

 

 

Отчет по лабораторному проекту №1

“Стандартные алгоритмы LU-разложения”

Вариант №2:

- разложение на основе гауссова исключения по столбцам

с выбором главного элемента по активной подматрице

 

 

Выполнил: Багаутдинов Н.Р.

гр. ИТСС-О-12/1

 

_______________________

 

 

Проверила: доц. Цыганова Ю.В.

 

_______________________

 

 

Ульяновск – 2014 г.


 

Цель проекта: приобретение практических навыков реализации алгоритмов LU-разложения в виде программного проекта на языке MATLAB.

 

Задания:

1. Система меню для взаимодействия пользователя с программой и генерация исходных данных задачи.

2. Функция факторизации матрицы.

3. Функция решения системы линейных алгебраических уравнений.

4. Функция вычисления определителя матрицы.

5. Функция обращения матрицы через решение системы AX = E.

6. Функция обращения матрицы через элементарные преобразования разложения.

7. Эксперимент 1 “Решение СЛАУ для случайных матриц”.

8. Эксперимент 2 “Решение СЛАУ с плохо обусловленными матрицами”.

9. Эксперимент 3 “Обращение случайных матриц”.

 

Задание №1 “Система меню для взаимодействия пользователя с программой и генерация исходных данных задачи”

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

 

Программный код меню:

clc;

disp('Лабораторный проект №1 "Стандартные алгоритмы LU-разложения"');

disp('Вариант №3');

disp('Выполнил студент гр. ИТСС-О-12/1 Багаутдинов Н. Р.');

while(1)

param = input('Нажмите "1" для ввода данных;\nНажмите "2" для обработки дан-ных;\nНажмите "3" для выхода из программы\n', 's');

switch param

case '1',

A = input('Введите матрицу: \n A = \n');

n_ar =size(A);

n = n_ar(1);

b = input('Введите вектор: \n b = \n');

n_ar_b =size(A);

n_b = n_ar_b(1);

if (n_b ~= n)

disp('Размерности A и b не совпадают!');

return;

end

case '2',

if (exist('A', 'var') == 0)

disp('Введите данные для обработки!');

else

A_copy=A;

%____________________________________________________________

disp('Факторизация:');

[A, p, q, znak] = factorization(A, n);

disp('Факторизация матрицы A с учетом перестановок');

% PA - матрица с учётом перестановок

PA=printPA(A,n,p,q)

disp('Нижняя треугольная часть L');

L=tril(PA)

disp('Верхняя треугольная часть U');

U=triu(PA)-eye(n).*PA+eye(n)

B=L*U

%__________________________________________________________

% x=zeros(n,1);

% x1=zeros(n,1);

% disp('Решение СЛАУ');

% x=SLAE(A, b, n, p)

% disp('Проверка решения');

% x1=A_copy\b

% %___________________________________________________________

% disp('Детерминант матрицы A');

% d=determinant(A, n, p, znak)

% disp('Проверка детерминанта');

% d1=det(A_copy)

% %__________________________________________________________

% A11=zeros(n,n);

% A12=zeros(n,n);

% disp('Обращение матрицы через решение системы AX=E');

% A11=inversion_1(A, n, p)

% disp('Обращение матрицы через элементарные преобразования разложения');

% A12=inversion_2(A, n, p)

% disp('Проверка обращения');

% V=inv(A_copy)

end

case '3', return;

otherwise, disp('Введен некорректный символ. Введите один из предложенных.');

end

end

 

 

 

Рис 1. Пример работы меню

 

Задание №2 “Функция факторизации матрицы”

 

Алгоритм 1. - разложение на основе гауссова исключения по активной подматрице:

Для k =1 до n

Выбираем главный элемент в A ( k -1).

Нормируем первую строку матрицы A ( k -1).

Для i = k +1 до n

Вычитаем первую строку матрицы A ( k -1), умноженную на aik ( k -1),
из i -й строки.

 

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

 

Здесь в качестве

главного выбирается максимальный по модулю элемент активной подмат-

рицы. В общем случае, чтобы поставить этот элемент на место диагональ-

ного, требуется обменивать столбцы и строки матрицы A(k−1), что связано с

введением двух дополнительных векторов: в одном хранятся перестановки

столбцов, а в другом перестановки строк матрицы A.

.

 

Программный код функции факторизации ( - разложение ):

function [A, p, q, znak, op_numb] = factorization(A, n)

p = zeros(n,1);

q = zeros(n,1);

%вспомогательные переменные

max = 0;

imax = 0;

jmax = 0;

buf = 0;

znak=1;

op_numb=0;

 

% заполняем массивы перестановок начальными значениями перед нача-лом факторизации

for i=1:n

p(i)=i;

q(i)=i;

end

% этот фрагмент вставляем в алгоритм факторизации

% n – номер задачи, k – номер текущего шага алгоритма

for k=1:n

%_____________________________________________________

max=abs(A(p(k),q(k)));

imax=k;

jmax=k;

for i=k:n

for j=k:n

if(abs(A(p(i),q(j)))>max)

max=abs(A(p(i),q(j)));

imax=i;

jmax=j;

end

end

end

% обмен

if(imax ~= k)

znak=-znak;

buf=p(k);

p(k)=p(imax);

p(imax)=buf;

end

if(jmax ~= k)

znak=-znak;

buf=q(k);

q(k)=q(jmax);

q(jmax)=buf;

end

%____________________________________________________

 

div = A(p(k),q (k));

for j=(k+1):n

A(p(k),q(j))= A(p(k),q(j))/div;%Нормировка

op_numb = op_numb + 1;

end

for i=(k+1):n

mult = A(p(i),q (k));

 

for j=(k+1):n

A(p(i),q (j))= A(p(i),q (j))- mult *A(p(k),q (j));

op_numb = op_numb + 1;

end

end

end

 

 

Рис 2. Вывод результатов факторизации

 

 

Задание №3 “Функция решения системы линейных алгебраических уравнений”

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

 

Программный код функции решения СЛАУ:

 

function [x, op_numb] = SLAE(A, b, n, p, q)

 

y = zeros(n,1);

x = zeros(n,1);

op_numb = 0;

 

for i=1:n

sum=0;

for j=1:i

sum =sum+A(p(i),q(j))*y(j);

op_numb = op_numb + 1;

end

y(i)=(b(p(i))-sum)/A(p(i),q(i));

op_numb = op_numb + 1;

end

 

for i=n:-1:1

sum=0;

for j=i+1:n

sum =sum+A(p(i),q(j))*x(q(j));

op_numb = op_numb + 1;

end

x(q(i))=y(i)-sum;

end

end

 

 

 

Рис 3. Вывод результатов решения СЛАУ

 

Задание №4 “Функция вычисления определителя матрицы”

Определитель находится как произведение диагональных элементов в разложенной матрице, поскольку

(определитель треугольной матрицы с единичной диагональю равен 1, а определитель треугольной матрицы равен произведению диагональных элементов).

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

 

Программный код функции вычисления определителя:

 

function [det] = determinant(A, n, p,q, znak)

det=1;

for i=1:n

det=det*A(p(i),q(i));

end

det= det*znak;

end

 

 

Рис 4. Вывод результатов вычисления определителя

 

Задание №5 “Функция обращения матрицы через решение системы AX=E”

Решение осуществляется исходя из равенства AX=E, где X – обратная матрица, E – единичная матрица. Разобьём матрицы X и E на столбцы и запишем систему СЛАУ:

, где xkk -ый столбец матрицы X, ekk -ый столбец матрицы E.

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

 

Программный код функции обращения матрицы через решение системы AX=E:

 

function [V, count_ops_1] = inversion_1(A, n, p, q)

 

V = zeros(n,n);

I = eye(n);

count_ops_1 = 0;

count_ops_cur = 0;

 

for count = 1:n

b = I(:,count);

y = zeros(n,1);

x = zeros(n,1);

[x,count_ops_cur]=SLAE(A, b, n, p, q);

V(:,count) = x;

count_ops_1 = count_ops_1 +count_ops_cur;

end

 

end

 

 

 

 

Рис 5. Вывод результатов обращения матрицы через решение системы AX=E

 

Задание № 6 “Функция обращения матрицы через элементарные
преобразования”

Пусть . Тогда . Матрицы L и уже известны. Следовательно, сначала найдём обратные к ним, а затем найдём их произведение и получим матрицу, обратную к A.

 

Программный код функции обращения матрицы через элементарные преобразования разложения:

 

function [A_ob, count_ops_2] = inversion_2(A, n, p,q)

 

count_ops_2 = 0;

% Первый этап - подготовка (обращение элементарных матриц)

 

for i=1:n

for j=(i+1):n

A(p(i),q(j))=-A(p(i),q(j));

end

end

 

for j=1:n

A(p(j),q(j))=1/A(p(j),q(j));

count_ops_2=count_ops_2+1;

for i=(j+1):n

A(p(i),q(j))=-A(p(i),q(j))*A(p(j),q(j));

count_ops_2=count_ops_2+1;

end

end

 

disp('Результат выполнения первого этапа');

PA=printPA(A,n,p,q)

 

% Считаем матрицу U^{-1}

for k=n:(-1):2

for i=1:(k-2)

for j=k:n

A(p(i),q(j))= A(p(i),q(j))+A(p(i),k-1)*A(p(k-1),q(j));

count_ops_2=count_ops_2+1;

end

end

end

 

 

% Считаем матрицу L^{-1}

 

for k=1:(n-1)

for i=(k+2):n

for j=1:k

A(p(i),q(j))=A(p(i),q(j)) + A(p(i),k+1)*A(p(k+1),q(j));

count_ops_2=count_ops_2+1;

end

end

for j=1:k

A(p(k+1),q(j))=A(p(k+1),q(j)) * A(p(k+1),q(k+1));

count_ops_2=count_ops_2+1;

end

end

 

disp('Результат выполнения второго этапа');

PA=printPA(A,n,p,q)

 

% Перемножаем в одной матрице U^{-1} на L^{-1}

 

for i=1:n

for j=1:n

if(i<j)

sum=0;

for k=j:n

sum = sum +A(p(i),k)*A(p(k),q(j));

count_ops_2=count_ops_2+1;

end

end

if(i>=j)

sum =A(p(i),q(j));

for k=i+1:n

sum = sum +A(p(i),k)*A(p(k),q(j));

count_ops_2=count_ops_2+1;

end

end

A(p(i),q(j))=sum;

end

end

 

 

disp('Результат выполнения третьего этапа');

% Учитываем перестановки срок

A_ob=printPA1(A,n,p,q);

 

end

 

Вспомогательная функция для учёта обратной перестановки столбцов:

 

% Вспомогательная функция для учёта обратной перестановки столбцов:

 

function [ PA1 ] = printPA1(A,n,p,q)

% Учитываем перестановки строк при обращении

 

% Вектор обратных перестановок строк

p1=zeros(n,1);

for i=1:n

p1(p(i))=i;

end

 

for i=1:n

for j=1:n

PA1(i,j)=A(p(i),p1(j));

end

end

 

end

 

 

Рис 6. Вывод результатов обращения матрицы через элементарные преобразования разложения

 

Задание № 7 “Эксперимент 1 “Решение СЛАУ”

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

 

Программный код эксперимента №1:

function [] = exp1()

clc();

frmt = get(0,'Format');

frmtV = get(0,'FormatSpacing');

format('shortG')

format('compact')

exp_Data = zeros(1,100);

count = 1;

for n=5:5:100

exp_Data(count) = n;

count = count+1;

A = randi([-100,100],n,n);

A_copy=A;

b = randi([-100,100],n,1);

% Засекаем время

tic

[A, p,q, znak, op_numb] = factorization(A, n);

[x, op_numb_2] = SLAE(A, b, n, p, q);

exp_Data(count) = toc;

count = count+1;

 

x_ex = zeros(n,1);

b_ex = zeros(n,1);

for i=1:n

x_ex(i) = i;

end;

 

% Вычисляем правую вектор b для точного решения

for i=1:n

for j=1:n

b_ex(i) = b_ex(i) + A_copy(i,j) * x_ex(j);

end

end

[x_calc] = SLAE(A, b_ex, n, p, q);

 

% Вычисляем точность найденного решения

for i=1:n

ex = abs(x_ex(i) - x_calc(i));

if (i == 1)

max_ex = ex;

else

if (ex > max_ex)

max_ex = ex;

end

end

end

exp_Data(count) = max_ex;

count = count +1;

 

exp_Data(count) = n^3/3;

count = count +1;

 

exp_Data(count) = op_numb + op_numb_2;

count = count + 1;

end

 

dlmwrite('file.dat',exp_Data,';');

clear A;

clear exp_Data;

 

n_ar = zeros(1,20);

time_ar = zeros(1,20);

ex_ar = zeros(1,20);

th_numb_ar = zeros(1,20);

fact_numb_ar = zeros(1,20);

 

read_data = dlmread('file.dat', ';');

count = 1;

 

% Вывод таблицы

disp(' Порядок Время Точность Теоретическое ЧО Реальное ЧО');

for i=1:5:100

n_ar(count) = read_data(i);

time_ar(count) = read_data(i+1);

ex_ar(count) = read_data(i+2);

th_numb_ar(count) = read_data(i+3);

fact_numb_ar(count) = read_data(i+4);

count = count + 1;

str = [round(read_data(i)), read_data(i+1), read_data(i+2), round(read_data(i+3)), round(read_data(i+4))];

disp(str);

end

 

% Вывод графиков

plot(n_ar, time_ar, 'k');

title('Зависимость времени решения от размера матрицы A (мск)');

grid on;

figure;

 

plot(n_ar, ex_ar, 'k');

title('Зависимость точности решения от размера матрицы A');

grid on;

figure;

 

plot(n_ar, th_numb_ar, n_ar, fact_numb_ar,'--');

title('Реальное и теоретическое количество операций');

grid on;

legend('Теоретич. ЧО','Реальное ЧО',2);

 

format(frmt)

format(frmtV)

end

 

 

Рис 7. Таблица экспериментальных данных для решения СЛАУ

 

 

Рис 8. График зависимости времени решения от порядка матриц для решения СЛАУ

 

Рис 9. График зависимости точности решения от порядка матриц для решения СЛАУ

Рис 10. График зависимости реального и теоретического числа операций от размера матрицы для решения СЛАУ

 

Задание №8 “Эксперимент 2
“Решение СЛАУ с плохо обусловленными матрицами”

Цель данного эксперимента – исследовать скорость и точность решения СЛАУ с плохо обусловленными матрицами.

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

 



Поделиться:




Поиск по сайту

©2015-2024 poisk-ru.ru
Все права принадлежать их авторам. Данный сайт не претендует на авторства, а предоставляет бесплатное использование.
Дата создания страницы: 2016-04-27 Нарушение авторских прав и Нарушение персональных данных


Поиск по сайту: