Полезное:
Как сделать разговор полезным и приятным
Как сделать объемную звезду своими руками
Как сделать то, что делать не хочется?
Как сделать погремушку
Как сделать так чтобы женщины сами знакомились с вами
Как сделать идею коммерческой
Как сделать хорошую растяжку ног?
Как сделать наш разум здоровым?
Как сделать, чтобы люди обманывали меньше
Вопрос 4. Как сделать так, чтобы вас уважали и ценили?
Как сделать лучше себе и другим людям
Как сделать свидание интересным?
Категории:
АрхитектураАстрономияБиологияГеографияГеологияИнформатикаИскусствоИсторияКулинарияКультураМаркетингМатематикаМедицинаМенеджментОхрана трудаПравоПроизводствоПсихологияРелигияСоциологияСпортТехникаФизикаФилософияХимияЭкологияЭкономикаЭлектроника
|
Ульяновский государственный университетСтр 1 из 2Следующая ⇒ Государственное образовательное бюджетное учреждение Высшего профессионального образования Ульяновский государственный университет
Факультет Математики и информационных технологий Кафедра Информационных технологий
Отчет по лабораторному проекту №1 “Стандартные алгоритмы LU-разложения” Вариант №8: - разложение по компактной схеме Краута с выбором главного элемента по столбцу.
Выполнили: Антипкина Ю., Пятаков О. гр. МОАИС-О-14/1
_______________________
Проверила: доц. Цыганова Ю.В.
_______________________
Ульяновск – 2016 г.
Цель проекта: приобретение практических навыков реализации алгоритмов LU-разложения в виде программного проекта на языке MATLAB.
Задания: 1. Система меню для взаимодействия пользователя с программой и генерация исходных данных задачи. 2. Функция факторизации матрицы. 3. Функция решения системы линейных алгебраических уравнений. 4. Функция вычисления определителя матрицы. 5. Функция обращения матрицы через решение системы AX = E. 6. Функция обращения матрицы через элементарные преобразования разложения. 7. Эксперимент 1 “Решение СЛАУ для случайных матриц”. 8. Эксперимент 2 “Решение СЛАУ с плохо обусловленными матрицами”. 9. Эксперимент 3 “Обращение случайных матриц”.
Задание №1 “Система меню для взаимодействия пользователя с программой и генерация исходных данных задачи” Взаимодействие с пользователем осуществляется через командную строку МАТЛАВ. Пользователю на выбор предлагается выбрать одно из действий: ввести данные, обработать данные (выполнить факторизацию матрицы, решить СЛАУ, вычислить определитель, вычислить обратную матрицу) или выйти из программы. Если пользователь выбирает обработку данных, но данные не введены, то программа выводит сообщение “Введите данные для обработки!” и предложит выбрать действие заново.
Программный код меню: clc; disp('Лабораторный проект №1 "Стандартные алгоритмы LU-разложения"'); disp('Вариант №8 "LU - разложение по компактной схеме Краута с выбором главного элемента по столбцу." '); disp('Выполнили студенты гр. МОАИС-О-14/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, znak] = factorization(A, n); disp('Факторизация матрицы A с учетом перестановок'); % PA - матрица с учётом перестановок PA=printPA(A,n,p) disp('Нижняя треугольная часть L'); L=tril(PA) disp('Верхняя треугольная часть U'); U=triu(PA)-eye(n).*PA+eye(n) %______________________________________________________ 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 По формуле вычисляем k-ю строку матрицы U. Выбираем среди элементов k-й строки главный элемент. По формуле 1 вычисляем k-й столбец матрицы L. Формула 1.
Выбираем главный элемент:
Стратегия выбора главного элемента по столбцу. В качестве главного на k -м шаге метода Гаусса выбирается максимальный по модулю элемент первого столбца активной подматрицы. Затем этот элемент меняется местами с диагональным элементом, что соответствует перестановке строк матрицы A ( k –1). Строки матрицы A ( k –1) остаются на своих местах, переставляются только элементы вектора перестановок p, в котором хранятся номера строк исходной матрицы, т. е. элементы вектора перестановок. Все обращения к элементам нижней треугольной части L и верхней треугольной части матрицы A осуществляются через вектор перестановок. По формуле 2 вычисляем k – ю строку матрицы U. Формула 2 Программный код функции факторизации ( - разложение по схеме Краута с выбором главного элемента по столбцу):
function [A, p, znak, op_numb] = factorization(A, n)
p = zeros(n,1); max = 0; imax = 0; buf = 0; znak=1; op_numb = 0; % Заполняем массивы перестановок начальными значениями перед началом факторизации for i=1:n p(i)=i; end
for k=1:n for i=k:n s=0; for r=1:(k-1) s=s+A(p(i),r)*A(p(r),k); end A(p(i),k)=A(p(i),k)-s; end
%---------Выбор главного элемента--------------
max=abs(A(p(k),k)); imax=k; for i=(k+1):n
if(abs(A(p(i),k))>max) max=abs(A(p(i),k)); imax=i; end end % Обмен if(imax ~= k) znak=-znak; buf=p(k); p(k)=p(imax); p(imax)=buf; end %----------------------------------------------
for j=(k+1):n s=0; for r=1:(k-1) s=s + A(p(k),r)*A(p(r),j); end A(p(k),j)=(A(p(k),j)-s)/A(p(k),k); end end
Вспомогательная функция для учёта перестановок строк в разложенной матрице:
function [ PA ] = printPA(A, n, p) % Вспомогательная функция for i=1:n for j=1:n PA(i,j)=A(p(i),j); end end
end
Рис 2. Вывод результатов факторизации
Задание №3 “Функция решения системы линейных алгебраических уравнений” Поиск решения СЛАУ осуществляется по разложению. Решение находим в два этапа: сначала находится решение для системы с нижней треугольной матрицей, затем находится решение с верхней треугольной матрицей. Математическое обоснование:
Программный код функции решения СЛАУ:
function [x, op_numb] = SLAE(A, b, n, p)
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),j)*y(j); op_numb = op_numb + 1; end y(i)=(b(p(i))-sum)/A(p(i),i); op_numb = op_numb + 1; end
for i=n:-1:1 sum=0; for j=i+1:n sum =sum+A(p(i),j)*x(j); op_numb = op_numb + 1; end x(i)=y(i)-sum; end
end
Рис 3. Вывод результатов решения СЛАУ
Задание №4 “Функция вычисления определителя матрицы” Определитель находится как произведение диагональных элементов в разложенной матрице, поскольку (определитель треугольной матрицы с единичной диагональю равен 1, а определитель треугольной матрицы равен произведению диагональных элементов). Так как используется первая стратегия выбора главного элемента, то при работе алгоритма факторизации происходит обмен строк. При этом каждый раз при обмене знак определителя меняется на противоположный. Для учёта знака в программе используется переменная-флаг znak, в котором сохраняется значение -1 или 1, в зависимости от того, чётное или нечётное количество перестановок было сделано на данный момент.
Программный код функции вычисления определителя:
function [det] = determinant(A, n, p, znak) det=1; for i=1:n det=det*A(p(i),i); end det= det*znak; end
Рис 4. Вывод результатов вычисления определителя
Задание №5 “Функция обращения матрицы через решение системы AX=E” Решение осуществляется исходя из равенства AX=E, где X – обратная матрица, E – единичная матрица. Разобьём матрицы X и E на столбцы и запишем систему СЛАУ: , где xk – k -ый столбец матрицы X, ek – k -ый столбец матрицы E. Учитывая, что известно -разложение матрицы A, каждая из СЛАУ решается алгоритмом, описанном в задании № 3. Столбцы обратной матрицы находятся по очереди и затем собираются вместе.
Программный код функции обращения матрицы через решение системы AX=E:
function [V, count_ops_1] = inversion_1(A, n, p) 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); 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)
count_ops_2 = 0; % Первый этап - подготовка (обращение элементарных матриц)
for i=1:n for j=(i+1):n A(p(i),j)=-A(p(i),j); end end
for j=1:n A(p(j),j)=1/A(p(j),j); count_ops_2=count_ops_2+1; for i=(j+1):n A(p(i),j)=-A(p(i),j)*A(p(j),j); count_ops_2=count_ops_2+1; end end
disp('Результат выполнения первого этапа'); PA=printPA(A,n,p)
% Считаем матрицу U^{-1} for k=n:(-1):2 for i=1:(k-2) for j=k:n A(p(i),j)= A(p(i),j)+A(p(i),k-1)*A(p(k-1),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),j)=A(p(i),j) + A(p(i),k+1)*A(p(k+1),j); count_ops_2=count_ops_2+1; end end for j=1:k A(p(k+1),j)=A(p(k+1),j) * A(p(k+1),k+1); count_ops_2=count_ops_2+1; end end
disp('Результат выполнения второго этапа'); PA=printPA(A,n,p)
% Перемножаем в одной матрице 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),j); count_ops_2=count_ops_2+1; end end if(i>=j) sum =A(p(i),j); for k=i+1:n sum = sum +A(p(i),k)*A(p(k),j); count_ops_2=count_ops_2+1; end end A(p(i),j)=sum; end end
disp('Результат выполнения третьего этапа'); % Учитываем перестановки срок A_ob=printPA1(A,n,p);
end
Вспомогательная функция для учёта обратной перестановки столбцов:
function [ PA1 ] = printPA1(A,n,p) % Учитываем перестановки строк при обращении
% Вектор обратных перестановок строк 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, znak, op_numb] = factorization(A, n); [x, op_numb_2] = SLAE(A, b, n, p); 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);
% Вычисляем точность найденного решения 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.
|