Eliminacja Gaussa bez wyboru elementu głównego:
function [L, U] = bezwyboru(A)
b=size(A);
b=b(1);
U=A;
L=eye(size(A));
for a=1:(b-1)
for i=a:(b-1)
L(i+1,a)=U(i+1,a)/U(a,a);
U(i+1,:)=U(i+1,:)-U(i+1,a)/U(a,a)*U(a,:);
end
end
Eliminacja Gaussa z częściowym wyborem elementu głównego:
function [L, U] = zwyborem(A)
b=size(A);
b=b(1);
U=A;
L=eye(b);
for a=1:(b-1)
C=0;
for j=1:(b-a+1)
C(j)=abs(U(j-1+a,a));
end
for j=a:b
if max(C)==abs(U(j,a))
index=j;
end
end
if max(C)~=0
D=U(a,:);
U(a,:)=U(index,:);
U(index,:)=D;
for k=-a+1:-1
D=L(a,-k);
L(a,-k)=L(index,-k);
L(index,-k)=D;
end
for i=a:(b-1)
L(i+1,a)=U(i+1,a)/U(a,a);
U(i+1,:)=U(i+1,:)-U(i+1,a)/U(a,a)*U(a,:);
end
end
end
Pomiar czasu
function [iloczyn,wyznacznik,slad,LU]=czas(N)
A=tstcnd(N)
B=tstcnd(N)
tic
for i=1:100
end
toc
p=toc
tic
for i=1:100
C=A*B
end
toc
a=toc
iloczyn=a/100-p;
tic
for i=1:100
det(A)
end
b=toc
wyznacznik=b/100-p;
tic
for i=1:100
trace(A)
end
toc
c=toc
slad=c/100-p;
tic
for i=1:100
lu(A)
end
toc
d=toc
LU=d/100-p;
disp(iloczyn)
disp(wyznacznik)
disp(slad)
disp(LU)
| Wielkość macierzy | Wyznacznik | Ślad | Iloczyn | LU |
|---|---|---|---|---|
| 2 | 8.7980e-005 | 1.0470e-004 | 2.3030e-005 | 1.0910e-004 |
| 4 | 2.5573e-005 | 1.9276e-005 | 8.3811e-005 | 2.4059e-004 |
| 8 | 1.6075e-004 | 8.4109e-005 | 0.0020 | 2.4363e-004 |
| 16 | 4.8006e-005 | 2.1959e-005 | 7.8761e-004 | 7.6916e-004 |
| 32 | 5.3571e-005 | 2.1010e-005 | 0.0039 | 0.0039 |
Macierz dobrze uwarunkowana
A=tstmat(5)
A =
77 -20 -44 46 -198
-3 -106 -126 -187 -346
-5 201 -166 -128 -215
-118 -33 375 -26 49
151 -87 246 117 58
Cond(A)= 8.9687
x=[-3 ;8 ;-4; -2; -4]
b=[485; 1423 ;3403; -1554; -2599]
function []= zaburzenie1()
for i=1:30
%A-macierz ,b-wektor wynikowy
A=[77 -20 -44 46 -198;-3 -106 -126 -187 -346;-5 201 -166 -128 -215; -118 -33 375 -26 49;151 -87 246 117 58];
x=[-3 ;8 ;-4; -2; -4];
b=[485; 1423 ;3403; -1554; -2599];
%zmiana w A
macierzA=inv(A)*b;
epsilon(1,i)=0.001*i;
A(1,1)=(A(1,1)*(1+epsilon(1,i)));
macierzA2=inv(A)*b;
zab(1,i)=norm(macierzA-macierzA2);
%zmiana elementu w b
b(1,1)=(b(1,1)*(1+epsilon(1,i)));
macierzA3=inv(A)*b
zab2(2,i)=norm(macierzA-macierzA3);
end
%wykresy
subplot(2,1,1);plot(epsilon(1,:),zab(1,:),'.r')
xlabel('zaburzenie');
ylabel('norma macierzy');
title('Zaburzenie macierzy (dobrze uwarunkowana)');
subplot(2,1,2);plot(epsilon(1,:),zab2(2,:),'.k')
xlabel('zaburzenie');
ylabel('norma macierzy');
title('Zaburzenie wektora b');
end
:
Macierz źle uwarunkowana
A=tstcnd(5)
A =
-625964 -21959 -119877 -455246 792434
611825 20965 115874 442731 -771741
18887 474 3127 12892 -22852
309459 10879 59323 225163 -391882
617222 21466 117718 448054 -780323
Cond(A)= 8.7748e+006
x =[-3;8;-4;-2;-4]
b=[-77516; 70251;247; 3856;74374]
m-plik
function []= metody3b()
for i=1:30
%A-macierz ,b-wektor wynikowy
A=[-625964 -21959 -119877 -455246 792434; 611825 20965 115874 442731 -771741; 18887 474 3127 12892 -22852; 309459 10879 59323 225163 -391882; 617222 21466 117718 448054 -780323];
x =[-3;8;-4;-2;-4];
b=[-77516; 70251;247; 3856;74374];
%zmiana w A
macierzAA1=inv(A)*b;
epsilon(1,i)=0.001*i;
A(1,1)=(A(1,1)*(1+epsilon(1,i)));
macierzAA2=inv(A)*b;
zab(1,i)=norm(macierzAA1-macierzAA2);
%zmiana elementu w b
b(1,1)=(b(1,1)*(1+epsilon(1,i)));
macierzAA3=inv(A)*b
zab2(2,i)=norm(macierzAA1-macierzAA3);
%wykresy
subplot(2,1,1);plot(epsilon(1,:),zab(1,:),'.r')
xlabel('zaburzenie');
ylabel('norma macierzy');
title('Zaburzenie macierzy (źle uwarunkowana)');
subplot(2,1,2);plot(epsilon(1,:),zab2(2,:),'.m')
xlabel('zaburzenie');
ylabel('norma macierzy');
title('Zaburzenie wektora b');
end