Сравнение методов решения систем линейных уравнений
1. Представление матрицы
1.
LU-разложение - представление матрицы A в виде LU, где L - нижняя треугольная
матрица с диагональными элементами, равными единице, а U - верхняя треугольная
матрица. LU-разложение еще называют LU-факторизацией. Алгоритм LU-разложения
лежит в основе широко распространённого метода решения систем линейных
алгебраических уравнений (СЛАУ) - метода Гаусса. Матрица L является
нижнетреугольной, с единичной диагональю, поэтому её определитель равен 1.
Матрица U - верхнетреугольная матрица, значит её определитель равен
произведению элементов, расположенных на главной диагонали.
2.
QR-разложение матрицы - представление матрицы в виде произведения унитарной
(или ортогональной матрицы) и верхнетреугольной матрицы. где Q - унитарная
матрица размера, а R - верхнетреугольная матрица размера. В частном случае,
когда матрица A состоит из действительных чисел, Q является ортогональной
матрицей (то есть QTQ = I, где I - единичная матрица). По аналогии, можно
определить варианты этого разложения: QL-, RQ-, и LQ-разложения, где L -
нижнетреугольная матрица.
3.
Разложе́ние Холе́цкого - представление симметричной положительно-определённой матрицы A
в виде A = LLT, где L - нижняя треугольная матрица со строго положительными
элементами на диагонали. Иногда разложение записывается в эквивалентной форме:
A = UTU, где U = LT - верхняя треугольная матрица. Разложение Холецкого всегда
существует и единственно для любой симметричной положительно-определённой
матрицы. Существует также обобщение этого разложения на случай
комплекснозначных матриц. Если A - положительно-определённая эрмитова матрица,
то существует разложение, где L - нижняя треугольная матрица с положительными
действительными элементами на диагонали, а - эрмитово-сопряжённая к ней
матрица.
Листинг программы
clear all;
n=100;
A1=rand (n,
n);=rand(n);=A1*X1;=rand(n);=rand(n);=A2*X2;i=1:n (i, i)=rand(1); =orth
(rand(n));=S31*J31*(S31');=rand (n, 1);=A3*X3;i=1:10,
[S1,
J1]=eig(A1);(1,1)=J1 (1,1)*10^(i);=S1*J1*S1^(-1);=T1*X1;
[L, U]=lu(T1);=inv(U)*inv(L)*B11;(i)=norm
(X11-X1)/norm(X1);
[S2,
J2]=eig(A2);(1,1)=J2 (1,1)*10^(i);=S2*J2*S2^(-1);=T2*X2;
[Q,
R]=qr(T2);=inv(R)*inv(Q)*B21;(i)=norm (X21-X2)/norm(X2);=S31;=J31;(1,1)=J3
(1,1)*10^i;=S3*J3*S3';=T3*X3;=chol(T3);=inv(R)*inv(R')*B31;(i)=norm (X31-X3)/norm(X3);(pogr1,'green');on;on;(pogr2,'red');
plot (pogr3,'blue');
Результаты и анализ
Матрица 10*10
Матрица 20*20
Матрица 100*100
При решении систем линейных
алгебраических уравнений, метод Холецкого более точный и дает меньшую погрешность.
Причиной этого может служить тот факт, что для применения данного разложения
необходимо использовать матрицу специального вида (симметричная, положительно
определенная). В остальных же двух методах, QR разложение обычно
давало результат лучше, чем LU разложение. Также характерной чертой для всех видов разложения
является тот факт, что при увеличении размерности матрицы, погрешность
возрастает.
3. Приближенные методы
решения алгебраических систем
Метод Зейделя
Чтобы пояснить суть метода,
перепишем задачу в виде:
Здесь в j-м уравнении мы перенесли в
правую часть все члены, содержащие xi, для i > j. Эта запись может быть
представлена:
где в принятых обозначениях D
означает матрицу, у которой на главной диагонали стоят соответствующие элементы
матрицы A, а все остальные нули; тогда как матрицы U и L содержат верхнюю и
нижнюю треугольные части A, на главной диагонали которых нули.
Итерационный процесс в методе
Гаусса-Зейделя строится по формуле
после выбора соответствующего
начального приближения. 
Метод Гаусса-Зейделя можно рассматривать
как модификацию метода Якоби. Основная идея модификации состоит в том, что
новые значения 
используются
здесь сразу же по мере получения, в то время как в методе Якоби они не
используются до следующей итерации:
где 
Таким образом i-тая компонента (k +
1) - го приближения вычисляется по формуле:
Условие сходимости
Приведём достаточное условие
сходимости метода. Теорема:
Пусть 
, где 
-
матрица, обратная к 
. Тогда
при любом выборе начального приближения 
:
5.
скорость сходимости метода равна скорости сходимости
геометрической прогрессии со знаменателем 
;
6.
верна оценка погрешности: 
.
Условие окончания
Условие окончания итерационного
процесса Зейделя при достижении точности в упрощённой форме имеет вид:
Более точное условие окончания
итерационного процесса имеет вид
и требует больше вычислений. Хорошо
подходит для разреженных матриц.
Метод простых итераций
Введя в рассмотрение
матрицы, эту систему можно коротко записать в виде матричного уравнения 
Предполагая, что
диагональные коэффициенты 
разрешим первое
уравнение системы (1) относительно x1,
второе - относительно x2,
и т.д. Тогда получим эквивалентную
систему
где
Любое (k+1) - е приближение вычисляют по формуле
То этот предел решением
второй системы. Запишем формулы приближений в развернутом виде:
Листинг программы
clear all
clc
%паарметры
dx=10^(-5);=10;=0;=0;=diag
(diag(rand (n, n)));=orth (rand(n, n));
%Задание числа
обусловленностей матирцы и цикл обработки
for
k=0:4=i+1;=J1;(1,1)=J (1,1)*10^k;=S*J*(S');
%Зейдель=eye(n);=rand (n, 1);_old=zeros (n,
1);=A*X;=A;=A;Index_1=1:nIndex_2=1: Index_1(Index_1,
Index_2)=0;Index_2=(Index_1+1):n(Index_1,
Index_2)=0;=inv(K);=-InvK*M;=InvK*z;_new=C;norm (X-x_new,
inf)>dx_new=B*x_old+C;_old=x_new;=step+1;
%ограничение на
количество итерацийstep>=10^8
stop;(i)=step;
XX(i)=k;=0;
%Простые Итерации
A=A/(norm(A)*1.1);=E-A;_old=zeros
(n, 1);=norm(B);q>=0.999999;=A*X;
x_new=z;
%высчитывание корней с
заданной точностью
while
norm (X-x_new, inf)>dx_new=B*x_old+z;_old=x_new;=step+1;
%ограничение на
количество итерацийstep>=10^8
stop;(i)=step;(i)=k;=0;(XX,
YY, 'red');on;(XX1, YY1,'blue');off;;
Метод
Прогонки
clear
all;=10;i=1: (n-1)(i, i+1)=rand (1,1);(i+1, i)=rand (1,1);(i, i)=A (i, i+1)+A
(i+1, i)+rand (1,1);(n, n)=A (n-1, n)+rand (1,1);=rand (n, 1);=A*X;(1)=(-1)*A
(2,1)/A (1,1);(1)=B(1)/A (1,1);i=2: (n-1)(i)=(-1)*A (i, i+1)/(delta (i-1)*A (i,
i-1)+A (i, i));(i)=(B(i) - A (i, i-1)*alpha (i-1))/(delta (i-1)*A (i, i-1)+A
(i, i));(n)=0;(n)=(B(i) - A (i+1, i)*alpha (i-1))/(delta (i-1)*A (i+1, i)+A (i,
i));(n, 1)=alpha(n);i=1:n-1(n-i, 1)=delta (n-i)*X1 (n-i+1,1)+alpha(i);=norm
(X1-X)/norm(X);
матрица листинг холецкий зейдель