POLITECHNIKA LUBELSKA WYDZIAŁ ELEKTRYCZNY

KATEDRA INFORMATYKI

METODY NUMERYCZNE

Temat :Obliczanie wartości residuów, wskaźnika uwarunkowania oraz macierzy odwrotnej przy wykorzystaniu odpowiednich metod numerycznych

Opracował :

Chmarzyński Robert

Gr. E D 7.5

Lublin styczeń 1996

Temat projektu.

W układzie równań Ax = b

dokładnym rozwiązaniem jest wektor x^T = [1,-1]. Dodatkowo dane są dwa rozwiązania przybliżone: x^T₁ = [0,999 ; -1,001] , x^T₂ = [0,341 ; -0,087]

obliczyć wartości residuów r(x₁) , r(x₂),

wyznaczyć wskaźnik uwarunkowania cond (A) korzystając z normy maksimum oraz z macierzy

wyznaczyć macierz odwrotną A^-1 korzystając z odpowiednich metod numerycznych.

Punkt 1

Obliczanie wartości residuów.

Rozwiązując układ równań liniowych (ręcznie lub za pomocą komputera) popełnia się błędy zaokrągleń. Komputer dla każdego układu wykonuje nieraz wiele milionów działań a każde z nich wprowadza błąd zaokrąglenia. Aby oszacować błąd posługujemy się residuami.

Jeżeli x' jest obliczonym rozwiązaniem układu Ax = b to określamy dla niego residuum jako wektor:

r = b - Ax'

Obliczamy residuum dla rozwiązania przybliżonego x₁:

r₁ = b - Ax₁'

Obliczamy residuum dla rozwiązania przybliżonego x₂:

r₂ = b - Ax₂'

Punkt 2

Wyznaczenie wskaźnika uwarunkowania.

Wskaźnik euklidesowy uwarunkowania k₂(A) dla macierzy prostokątnej A o rozmiarach m x n i kolumnach niezależnych liniowo określa się wzorem:

Jeżeli m = n (tzn. macierz jest kwadratowa) to wzór przyjmuje postać:

gdzie ||A||₂ oznaczono normę macierzy która wyraża się wzorem:

Obliczam zgodnie z powyższym wzorem normę macierzy A:

Oraz normę macierzy odwróconej A^-1:

Obliczam wskaźnik uwarunkowania cond(A) :

Punkt 3

Wyznaczenie macierzy odwrotnej A^-1 za pomocą metod numerycznych.

Macierz odwrotną wyznaczyłem korzystając z zależności:

Do obliczania macierzy odwrotnej wykorzystałem program napisany na laboratorium metod numerycznych (zamieszczony na końcu sprawozdania).

a)Odwracanie macierzy metodą dekompozycji LU (opcja 4).

Rozwiązanie:

a[1,1]= 6.5900013652*10⁵

a[1,2]= -5.6300011663*10⁵

a[2,1]= -9.1300018914*10⁵

a[2,2]= 7.8000016159*10⁵

b)Odwracanie macierzy według wzoru definicyjnego z wykorzystaniem eliminacji Gaussa do liczenia wyznaczników (opcja 3).

Rozwiązanie dokładnie takie jak w punkcie (a).

Do obliczenia dokładnego rozwiązania równania Ax = b wykorzystałem procedury zamieszczone w opracowaniu „Metody numeryczne” Barona. Warszawa 1995.

Rozwiązanie:

x[1] = 1.0000005460

x[2] = -1. 0000007564

Otrzymane wyniki są jednakowe dla metod:

a)eliminacji Gaussa

b)Jacobiego

c)Gaussa - SeidlaTEKST ŹRÓDŁOWY PROGRAMU DO OBLICZANIA MACIERZY ODWROTNEJ METODĄ ELIMINACJI GAUSSA

I DEKOMPOZYCJI LU

program det_mac;

{$M 64000,00000,600000}

uses crt;

type tab=array[1..15,1..15] of real;

indeks=array[1..15] of byte;

var n,l_p,blad,opcja:byte;

A,B:tab;

wpis:boolean;

zero,zm:real;

{***"procedura stop" zatrzymuje ekran do naciini©cia klawisza ***}

procedure stop;

var z:char;

begin

writeln;

write('przycisnij dowolny klawisz...');

while keypressed do z:=readkey;

z:=readkey;

writeln

end;

{***"procedura czytaj" czyta elementy macierzy z klawiatury ***}

procedure czytaj(n:byte;var A:tab);

var i,j:byte;

begin

writeln('Wpisz elementy macierzy...');

writeln('format: a[wiersz,kolumna]...');

for i:=1 to n do

for j:=1 to n do

begin

write('a[',i,',',j,']=');

read(A[i,j])

end

end;

{***"procedura il_mac" mnoľy przez siebie macierze A i B ***}

procedure il_mac(A,B:tab;n:byte;var C:tab);

var i,j,k:byte;

il,s:real;

begin

for i:=1 to n do

for j:=1 to n do

begin

s:=0;

for k:=1 to n do

s:=s+A[i,k]*B[k,j];

C[i,j]:=s

end;

end;

{***"procedura wypisz" drukuje na ekranie elementy macierzy ***}

procedure wypisz(A:tab;n:byte);

var i,j:byte;

begin

writeln('Elementy sa drukowane w formacie : a[wiersz,kolumna]...');

writeln('Po kazdym wierszu zatrzymanie drukowania');

for i:=1 to n do

begin

for j:=1 to n do writeln('a[',i,',',j,']=',A[i,j]);

stop

end

end;

{***"procedura trans_mac" transponuje macierz wejsciowĄ ***}

procedure trans_mac(A:tab;n:byte;var B:tab);

var i,j:byte;

begin

for i:=1 to n do

for j:=1 to n do

B[i,j]:=A[j,i]

end;

{***"procedure min" eliminuje z macierzy A wiersz s i kolumn© h oraz ***}

{*** przeksztaaca macierz A tak,aby policzy† minor M[s,h] macierzy A ***}

procedure min(s,h,n:byte;A:tab;var B:tab);

var i,j,k,r:byte;

begin

k:=1;

for i:=1 to n do

begin

r:=1;

if k=s then k:=k+1;

for j:=1 to n do

begin

if r=h then r:=r+1;

B[i,j]:=A[k,r];

r:=r+1

end;

k:=k+1

end;

end;

{***"funkcja detA_def" liczy wyznacznik z macierzy A weddug definicji ***}

function detA_def(A:tab;n:byte):real;

var g:integer;

j:byte;

b:tab;

detB,w:real;

begin

if n=2 then detA_def:=A[1,1]*A[2,2]-A[1,2]*A[2,1]

else if n>2 then

begin

g:=1;

detB:=0;

for j:=1 to n do

begin

min(1,j,n-1,A,B);

w:=a[1,j]*detA_def(B,n-1)*g;

detB:=w+detB;

g:=g*(-1)

end;

detA_def:=detB

end;

end;

{**"procedura wybor_el_g" dokonuje cz©©ciowego wyboru elementu gg˘wnego ***}

{** podczas dekompozycji LU,zapami©tuje liczb© przestawieä(zmienna l_p) ***}

{** oraz usytuowanie wierszy po przestawieniu(zmienna przst) co jest wy-***}

{** korzystywane p˘«niej przy obliczaniu wyznacznika i odwracaniu macierzy*}

procedure wybor_el_g(n,k:byte;zero:real;var L,U:tab;var l_p:byte;

var przst:indeks;var osobliwa:boolean);

var A,B:array[1..15] of real;

h,e,wiersz,p:byte;

max:real;

begin

osobliwa:=false;

max:=abs(zero);

wiersz:=k;

for h:=k+1 to n do { *************** }

if abs(U[h,k])>max then begin { wyb˘r elementu }

max:=abs(U[h,k]); { gg˘wnego }

wiersz:=h

end; { *************** }

if wiersz<>k then

begin

inc(l_p);

for e:=1 to n do { *************** }

begin

A[e]:=U[k,e];

U[k,e]:=U[wiersz,e]; { przestawienie }

U[wiersz,e]:=A[e]; { wierszy }

B[e]:=L[k,e];

L[k,e]:=L[wiersz,e];

L[wiersz,e]:=B[e];

end; { *************** }

p:=przst[wiersz]; { zapami©tanie }

przst[wiersz]:=przst[k]; { nowych ustawieä }

przst[k]:=p;

end else begin { *************** }

osobliwa:=true;exit

end;

end;

{*** "procedura wroc" przywraca szyk kolumn w macierzy odwr˘conej ***}

{*** przestawiony podczas dekompozycji LU macierzy wejjciowej przy***}

{*** wyborze elementu gg˘wnego,zapamietany przez zmiennĄ "przst" ***}

procedure wroc(A:tab;n:byte;przst:indeks;var B:tab);

var i,j:byte;

begin

move(A,B,sizeof(A));

for i:=1 to n do

for j:=1 to n do

if i<>przst[i] then B[j,przst[i]]:=A[j,i];

end;

{**"procedura dek_LU_z_wyb" dokonuje dekompozycji LU na macierzy wejjciowej*}

{** wewnatrz wywooywana jest procedura wyboru elementu gg˘wnego ************}

procedure dek_LU_z_wyb(A:tab;n:byte;zero:real;var L,U:tab;var l_p:byte;

var przst:indeks;var osobliwa:boolean);

var i,j,k,h,wiersz,e:byte;

c,suma,max:real;

d:array[1..15] of real;

begin

l_p:=0;

for i:=1 to n do przst[i]:=i;

osobliwa:=false;

move(A,U,sizeof(A));

for i:=1 to n-1 do

begin

if abs(U[i,i])<abs(zero) then wybor_el_g(n,i,zero,L,U,l_p,przst,osobliwa);

if osobliwa then exit;

for j:=i+1 to n do { ****************** }

begin { dekompozycja }

c:=U[j,i]/U[i,i]; { w/g zmodyfikowanej }

L[j,i]:=c; { eliminacji Gaussa }

for k:=1 to n do U[j,k]:=U[j,k]-c*U[i,k]; { ****************** }

end;

end;

if U[n,n]=0 then begin

osobliwa:=true;exit;

end;

for i:=1 to n do

for j:=1 to n do

if i=j then L[i,j]:=1 else

if j>i then L[i,j]:=0

end;

{***"procedura odwroc_mac_L" odwraca macierz tr˘jkĄtnĄ dolnĄ ***}

procedure odwroc_mac_L(L:tab;n:byte;var R:tab);

var i,j,d,k:byte;

licznik,g:real;

begin

fillchar(R,sizeof(R),0);

for d:=1 to n do

begin

k:=d;

for i:=1 to n-d+1 do

begin

if i=k then R[k,i]:=1 { indeksy i,j,k,n zgodne }

else if k>i then { z indeksami we wzorze }

begin { podanym w konspekcie }

licznik:=0;

for j:=i+1 to k do

licznik:=licznik+R[k,j]*L[j,i];

R[k,i]:=(-1)*licznik/L[i,i]

end;

inc(k);

end;

end;

end;

{***"procedura odwroc_mac_U" odwraca macierz tr˘jkĄtnĄ g˘rna ***}

procedure odwroc_mac_U(U:tab;n:byte;var R:tab);

var i,j,k,d:byte;

licznik,g:real;

begin

fillchar(R,sizeof(R),0);

for d:=1 to n do

begin

k:=d;

for i:=1 to n-d+1 do { indeksy i,j,k,n zgodne }

begin { z indeksami we wzorze }

if i=k then R[i,k]:=1/U[i,i] { podanym w konspekcie }

else if i<k then

begin

licznik:=0;

for j:=i+1 to k do

licznik:=licznik+R[j,k]*U[i,j];

R[i,k]:=(-1)*licznik/U[i,i]

end;

inc(k);

end;

end;

end;

{***"procedura odwroc_mac_LU" odwraca macierz wykorzystujĄc jej rozkkad ***}

{*** czynniki tr˘jkĄtne,informuje o osobliwooci macierzy:blad=1 ***********}

procedure odwroc_mac_LU(A:tab;n:byte;zero:real;var B:tab;var blad:byte);

var l_p:byte;

t_p:indeks;

osobliwa:boolean;

L,U:tab;

begin

dek_LU_z_wyb(A,n,zero,L,U,l_p,t_p,osobliwa);

if osobliwa then begin

osobliwa:=true;

exit

end;

odwroc_mac_L(L,n,L);

odwroc_mac_U(U,n,U);

il_mac(U,L,n,B);

wroc(B,n,t_p,B);

blad:=0;

end;

{***"funkcja detA_Gaus" liczy wyznacznik macierzy na podstawie macierzy***}

{*** tr˘jkĄtnej g˘rnej otrzymanej z dekompozycji LU w/g eliminacji Gaussa*}

function detA_Gaus(A:tab;n:byte;zero:real):real;

var l_p,i:byte;

t_p:indeks;

k:shortint;

s:real;

L:tab;

osobliwa:boolean;

begin

dek_LU_z_wyb(A,n,zero,L,A,l_p,t_p,osobliwa);

if osobliwa then begin

detA_Gaus:=0;

exit

end;

s:=1;

for i:=1 to n do

s:=s*A[i,i];

if odd(l_p) then k:=-1 else k:=1;

detA_Gaus:=k*s;

end;

{***"procedura odwroc_mac_def" odwraca macierz weddug wzoru definicyjnego***}

{*** ale do liczenia wyznacznik˘w uľywa funkcji detA_Gaus(w/g el. Gaussa****}

procedure odwroc_mac_def(A:tab;n:byte;zero:real;var C:tab;var blad:byte);

var h:real;

k,g:shortint;

i,j:byte;

B:tab;

begin

blad:=0;

h:=detA_Gaus(A,n,zero);

if h=0 then begin

blad:=1;

exit

end;

h:=1/h;

trans_mac(A,n,A);

k:=1;

for i:=1 to n do

begin

g:=k;

for j:=1 to n do

begin

min(i,j,n-1,A,B);

c[i,j]:=g*detA_Gaus(B,n-1,zero);

g:=g*(-1);

end;

k:=k*(-1);

end;

for i:=1 to n do

for j:=1 to n do

C[i,j]:=h*C[i,j]

end;

begin

clrscr;

writeln('>> Przykladowy program zastosowania metod numerycznych <<');

writeln('>> do obliczania wyznacznikow i odwracania macierzy... <<');

gotoxy(wherex,wherey+3);

write('Podaj rozmiar macierzy [3..15] n=');read(n);

write('Podaj ulamkowa liczbe reprezentujaca zero =');read(zero);

czytaj(n,A);

repeat

clrscr;gotoxy(wherex,wherey+2);

writeln('>> Mozliwe opcje do wyboru:');

writeln('>> 1.Liczenie wyznacznika w/g definicji.');

writeln('>> 2.Liczenie wyznacznika z eliminacji Gaussa.');

writeln('>> 3.Odwrocenie macierzy w/g wzoru definicyjnego z wykorzystaniem');

writeln('>> eliminacji Gaussa do liczenia wyznacznikow.');

writeln('>> 4.Odwrocenie macierzy z dekompozycji LU.');

writeln('>> 5.Wpisanie nowej macierzy.');

writeln('>> 6.Koniec programu.');

gotoxy(wherex,wherey+1);

writeln('>> Uwaga!!! Wybranie opcji 1 dla macierzy o duzym rozmiarze (n>8)');

writeln('>> powoduje bardzo dlugi czas oczekiwania na wynik.');

writeln('>> Uwaga!!! Zbyt duze "zero" powoduje w opcjach 2,3,4 ze detA=0');

gotoxy(wherex,wherey+1);

write('>> Twoj wybor:');read(opcja);

gotoxy(wherex,wherey+2);

case opcja of

1:begin

write('czekaj...');

zm:=detA_def(a,n);writeln;

writeln('Wyznacznik detA=',zm);stop;

end;

2:begin

writeln('Wyznacznik z eliminacji Gaussa detA=',detA_Gaus(a,n,zero));

stop;

end;

3,4:begin

write('czekaj...');

if opcja=3 then

odwroc_mac_def(a,n,zero,b,blad)

else odwroc_mac_LU(a,n,zero,b,blad);

writeln;

if blad=1 then begin

writeln('Macierz osobliwa...nie odwroce');stop;

end else begin

writeln('Macierz odwrocona...');

wypisz(b,n);

end;

end;

5:begin

write('Podaj rozmiar macierzy [3..15] n=');read(n);

write('Podaj ulamkowa liczbe reprezentujaca zero =');read(zero);

czytaj(n,A);

end;

6:begin

clrscr;exit

end;

end

until false

end.

PROCEDURY WYKORZYSTANE W OPRACOWANIU

ELIMINACJA GUSSA

PROCEDURE RRAL(VAR A :MAC;

VAR B,X :WEK;

N :BYTE;

EPS :FLOAT;

VAR BLAD :BYTE);

{ Rozwiazywanie liniowego ukladu rownan A*X=B }

{ metoda eliminacji Gaussa }

{ A - macierz kwadratowa rzedu N }

{ B - wektor wyrazow wolnych o N - elementach }

{ X - wektor rozwiazania o N - rozwiazaniach }

{ EPS - dokladnosc rozroznienia zerowego wiersza }

{ (np. EPS=1E-36) }

{ BLAD - nr bledu; 0 - brak bledu }

PROCEDURE RRAL;

VAR i,j,k :BYTE;

T,S :FLOAT;

BEGIN

BLAD:=0;

{ Skalowanie macierzy A i wektora wyrazow wolnych B }

SKALROW(A,B,N);

{ Konstrukcja ciagu macierzy A(i) wzory (1.20),(1.21) i (1.23)

oraz ciagu wektorow B(i) wzory (1.22) }

FOR i:=1 TO N DO

BEGIN

{ Wybor elementu glownego wg wzoru (1.24) }

T:=ABS(A[i,i]); k:=i;

FOR j:=i+1 TO N DO

IF ABS(A[j,i])>T THEN

BEGIN

T:=ABS(A[j,i]); k:=j { Zamiana wiersza k-tego z i-tym }

END;

IF T<EPS THEN

{ Warunek przerwania poszukiwania elementu glownego (1.27)}

{ Warunek nie istnienia rozwiazania rownania }

BEGIN

BLAD:=1;

EXIT

END;

IF i=k

THEN T:=A[i,i]

ELSE BEGIN

{ Zamiana elementu k-tego z i-tym wektora B }

T:=B[k]; B[k]:=B[i]; B[i]:=T;

{ Zamiana wiersza k-tego z i-tym macierzy A }

FOR j:=N DOWNTO i DO

BEGIN

T:=A[k,j]; A[k,j]:=A[i,j]; A[i,j]:=T

END

END;

T:=1/T; A[i,i]:=T;

FOR j:=i+1 TO N DO

BEGIN

S:=-A[j,i]*T; { wzor (1.23) }

B[j]:=B[j]+B[i]*S; { wzor (1.22) }

FOR k:=i+1 TO N DO

A[j,k]:=A[j,k]+A[i,k]*S; { wzor (1.21) }

END

END;

{ Rozwiazanie ukladu trojkatnego metoda postepowania

wstecz wzor (1.26) }

FOR i:=N DOWNTO 1 DO

BEGIN

T:=B[i];

FOR j:=i+1 TO N DO

T:=T-A[i,j]*x[j];

X[i]:=T*A[i,i]

END

END { RRAL };

JACOBIEGO

PROCEDURE RRALIJR(VAR A :MAC;

VAR b,x :WEK;

N :BYTE;

EPS,ALFA :FLOAT;

maxit :WORD;

VAR BLAD :BYTE);

{ Rozwiazywanie ukladu rownan liniowych A*x=b metoda }

{ iteracji prostej Jacobiego ALFA=1 wzory (1.59), (1.60) }

{ Richardsona ALFA<>1 wzor (1.62) }

{ A - macierz kwadratowa rzedu N }

{ b - wektor wyrazow wolnych o N elementach }

{ x - wektor rozwiazan - wymaga sie pewnych wartosci }

{ poczatkowych }

{ EPS - dokladnosc bezwzgledna iteracji (np. EPS=1E-8) }

{ ALFA - parametr regularyzacji Richardsona }

{ maxit - zadana maksymalna liczba iteracji }

{ BLAD - nr bledu; 0 - brak bledu }

PROCEDURE RRALIJR;

VAR i,j :BYTE;

k :WORD;

R,S,T:FLOAT;

y :WEK;

CH :CHAR;

BEGIN

SKALROW(A,b,N); { Skalowanie rownania A*x=b }

{ Przeksztalcenie macierzy A do postaci z mozliwie }

{ maksymalnymi co do modulu elementami na glownej }

{ przekatnej metoda przestawiania wierszy }

FOR i:=1 TO N DO

BEGIN

T:=ABS(A[i,i]); k:=i;

FOR j:=i+1 TO N DO

IF ABS(A[j,i])>T THEN

BEGIN

T:=ABS(A[j,i]);

k:=j

END;

IF i<>k THEN

BEGIN

T:=b[k]; b[k]:=b[i]; b[i]:=T;

FOR j:=1 TO N DO

BEGIN

T:=A[k,j]; A[k,j]:=A[i,j]; A[i,j]:=T

END

END

END;

R:=MWWMA(A,N,1E-4,ALFA,maxit,BLAD); k:=0;

writeln('najw jacobi wlasna',R);

{ Najwieksza co do modulu wartosc wlasna macierzy A1-ALFA*A }

IF R>=1

THEN BEGIN

BLAD:=6;

EXIT

END

ELSE BLAD:=0;

REPEAT

{ Iteracja prosta Jacobiego }

y:=x; R:=0.0; k:=k+1;

write('=',k);

FOR i:=1 TO N DO

BEGIN

S:=b[i];

FOR j:=1 TO N DO

S:=S-A[i,j]*y[j];

x[i]:=y[i]+ALFA*S; S:=ABS(S);

IF S>R THEN

R:=S

END;

UNTIL (R<EPS) OR (k>maxit);

IF k>maxit THEN

BLAD:=6;

END { RRALIJR };

GAUSSA - SEIDLA

PROCEDURE RRALIGS(VAR A :MAC;

VAR b,x :WEK;

N :BYTE;

EPS,OMEGA :FLOAT;

maxit :WORD;

VAR BLAD :BYTE);

{ Rozwiazywanie ukladu rownan liniowych A*x=b metoda }

{ iteracji Gaussa-Seidela OMEGA=1 wzory (1.64), (1.65) }

{ nadrelaksacji 0<OMEGA<>1 wzory (1.66), (1.67)

{ A - macierz kwadratowa rzedu N }

{ b - wektor wyrazow wolnych o N elementach }

{ x - wektor rozwiazan - wymaga sie pewnych wartosci

{ poczatkowych }

{ EPS - dokladnosc bezwzgledna iteracji (np. EPS=1E-8) }

{ OMEGA - parametr relaksacji }

{ maxit - zadana maksymalna liczba iteracji }

{ BLAD - nr bledu; 0 - brak bledu }

PROCEDURE RRALIGS;

VAR i,j :BYTE;

k :WORD;

R,S,T:FLOAT;

y :WEK;

CH :CHAR;

BEGIN

SKALROW(A,b,N); { Skalowanie rownania A*x=b }

{ Przeksztalcenie macierzy A do postaci z mozliwie }

{ maksymalnymi co do modulu elementami na glownej }

{ przekatnej metoda przestawiania wierszy }

FOR i:=1 TO N DO

BEGIN

T:=ABS(A[i,i]); k:=i;

FOR j:=i+1 TO N DO

IF ABS(A[j,i])>T THEN

BEGIN

T:=ABS(A[j,i]); k:=j

END;

IF i<>k THEN

BEGIN

T:=b[k]; b[k]:=b[i]; b[i]:=T;

FOR j:=1 TO N DO

BEGIN

T:=A[k,j]; A[k,j]:=A[i,j]; A[i,j]:=T

END

END

END;

R:=MWWMA(A,N,1E-4,OMEGA,maxit,BLAD); k:=0;

writeln('najw wartosc seidel wlasna = ',r);

{ Najwieksza co do modulu wartosc wlasna macierzy A1-OMEGA*A}

IF R>=1

THEN BEGIN

BLAD:=7;

EXIT

END

ELSE BLAD:=0;

REPEAT

y:=x; R:=0.0; k:=k+1;

write('=',k);

FOR i:=1 TO N DO

BEGIN

S:=b[i];

FOR j:=1 TO i-1 DO

S:=S-A[i,j]*x[j];

FOR j:=i TO N DO

S:=S-A[i,j]*y[j];

x[i]:=y[i]+OMEGA*S; S:=ABS(S);

IF S>R THEN

R:=S

END;

UNTIL (R<EPS) OR (k>maxit);

IF k>maxit THEN

BLAD:=8

END { RRALIGS };

---