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