plik

SPRAWOZDANIE #2

TEMAT:

Rozwiązywanie układu równań liniowych: eliminacja Gaussa, eliminacja
Gaussa-Jordana, eliminacja rozkładem LU. Obliczanie wyznacznika
macierzy. Wyznaczanie macierzy odwrotnej.

ROK/SEM.

2/3

IMIĘ I NAZWISKO:

Michał Kręzel

ODDANO:

2009.12.16

OCENA:

1. Cel ćwiczenia.

Seria ćwiczeń miała za zadanie opanować najbardziej popularne metody rozwiązywania układów

równań liniowych z jednym rozwiązaniem, a także przyswojenie rozkładu LU, co pozwalało na łatwe

obliczenie wyznacznika. Przy pomocy poznanych metod mogliśmy także łatwo wyznaczyć macierz

odwrotną.

2. Przy kodowaniu kolejnych funkcji zdefiniowano następujące struktury i funkcje pomocnicze:

2.1.

Struktury opisujące macierz współczynników i jednokolumnową macierz wyników równań.

Zmienna typu Twsp określa macierz współczynników w kolejnych równaniach dla kolejnych

zmiennych. Zmienna typu Twynik określa stałe po prawej stronie równań.

struct

Twsp

{

int

n, m;

//wiersze, kolumny

double

**t;

};

struct

Twynik

{

int

double

t[100];

};

2.2.

Funkcje odpowiedzialne za wczytywanie zmiennych typu Twsp i Twynik.

Twsp load_wsp(

double

tab[n_MAX][m_MAX],

int

n = 0,

int

m = 0)

{

Twsp tmp;

int

i, j;

(n != 0 && m != 0)

{

tmp.n = n;
tmp.m = m;
tmp.t =

new

double

*[n];

for

(i=0; i<m; i++) tmp.t[i] =

new

double

[m];

for

(i=0; i<n; i++)

for

(j=0; j<m; j++)

tmp.t[i][j] = tab[i][j];

}

1/9

else

{

cout <<

"rows: "

; cin >> tmp.n;

cout <<

"cols: "

; cin >> tmp.m;

tmp.t =

new

double

*[tmp.n];

for

(i=0; i<m; i++) tmp.t[i] =

new

double

[tmp.m];

for

(i=0; i<tmp.n; i++)

for

(j=0; j<tmp.m; j++)

{

cout <<

"give_["

<< i <<

"]["

<< j <<

"]: "

;

cin >> tmp.t[i][j];

}

return

tmp;

}

Twynik load_wynik(

double

tab[],

int

n=0)

{

Twynik tmp;

int

(n != 0)

{

tmp.n = n;

for

(i=0; i<n; i++)

tmp.t[i] = tab[i];

}

else

{

cout <<

"rows: "

; cin >> tmp.n;

for

(i=0; i<n; i++)

{

cout <<

"give_["

<< i <<

"]: "

;

cin >> tmp.t[i];

}

return

tmp;

}

2.3.

Funkcje odpowiedzialne za wypisywanie zmiennych Twsp i Twynik na ekran.

void

Print_Twsp(Twsp x,

char

* sep=

"\t"

bool

margins=

true

)

{

(margins) cout << endl;

for

(

int

i=0; i<x.n; i++)

{

for

(

int

j=0; j<(x.m - 1); j++)

{

cout << setw(6) << setprecision(3) << x.t[i][j] << sep;

}
cout << setw(6) << setprecision(3) << x.t[i][x.m - 1] << endl;

}

(margins) cout << endl;

return

;

}

2/9

void

Print_Twynik(Twynik x,

char

* sep=

"\t"

bool

margins=

true

)

{

(margins) cout << endl;

for

(

int

i=0; i<(x.n - 1); i++)

cout << setw(6) << setprecision(3) << x.t[i] << sep;

cout << setw(6) << setprecision(3) << x.t[x.n - 1] << endl;

(margins) cout << endl;

return

;

}

2.4.

Ponieważ bazę macierzy współczynników tworzy tablica dynamiczna, przekazywanie

zmiennych Twsp nawet przez wartość wiąże się ze zmianą elementów tej tablicy. W pewnych

przypadkach będziemy musieli użyć funkcji kopiowania zmiennych Twsp.

void

Copy_Twsp(Twsp A_original, Twsp &A)

{

for

(

int

i=0; i<A.n; i++)

for

(

int

j=0; j<A.m; j++)

A.t[i][j] = A_original.t[i][j];

}

3. Weźmy przykładowy układ równań:

{



Zatem:

[

1 2 3
3 1 2
2 3 1

]

[

1
2
3

]

4. Eliminacja Gaussa.

Metoda ta polega na doprowadzeniu macierzy współczynników układu równań (Twsp A) do

macierzy trójkątnej górnej (eliminacji zmiennych), po czym w 'postępowaniu odwrotnym' następuje obliczenie

kolejnych zmiennych za pomocą macierzy stałych równań (Twynik B). W algorytmie zastosowano zamianę

równań w przypadku, gdy współczynnik przy kolejnej zmiennej znaczącej jest równy zero.

4.1.

Kod:

Twynik Eliminacja_Gaussa(Twsp A, Twynik B)
{

int

licznik,i,j;

int

k = 0;

Twynik X; X.n = A.m;

//Eliminacja zmiennych

for

(licznik=1; licznik<A.n; licznik++)

{

3/9

(A.t[k][k] == 0)

{

double

max = A.t[k][k];

int

wiersz = k;

double

temp;

for

(

int

p=k; p<A.n; p++)

(fabs(A.t[p][k]) > fabs(max))

{max = A.t[p][k]; wiersz = p;}

//zamiana równań

for

(

int

p=0; p<A.m; p++)

{

temp = A.t[k][p];
A.t[k][p] = A.t[wiersz][p];
A.t[wiersz][p] = temp;

}
temp = B.t[k]; B.t[k] = B.t[wiersz]; B.t[wiersz] = temp;

}

for

(i=licznik; i<A.n; i++)

{

double

temp = A.t[i][k];

for

(j=k; j<A.m; j++)

A.t[i][j] = A.t[i][j] - temp*A.t[k][j]/A.t[k][k];

B.t[i] = B.t[i] - temp*B.t[k]/A.t[k][k];

}

//Print_Twsp(A);

k++;

}

//Postępowanie odwrotne

double

suma;

for

(i=(A.n - 1); i>=0; i--)

{

suma = 0;

for

(j=(i+1); j<A.m; j++)

suma += A.t[i][j] * X.t[j];

X.t[i] = (B.t[i] - suma)/A.t[i][i];

}

return

}

4.2.

Kroki algorytmu.

[

1 2 3
3 1 2
2 3 1

]

⇒

[

0 −5 −7
0 −1 −5

]

⇒

[

0 −5

−

3.6

]

X =

[

0.667 0.667 −0.333

]

5. Eliminacja Gaussa-Jordana.

Metoda ta bardzo przypomina eliminację Gaussa, ale macierz współczynników układu równań

(Twsp A) jest doprowadzana do macierzy, której wszystkie elementy poza przekątną są równe zero. Po

zakończeniu algorytmu, postępowanie odwrotne jest już zbyteczne. Bazę algorytmu stanowi

4/9

'wyjedynkowanie' pierwszego współczynnika kolejnego równania i postępowanie Jordana dla wszystkich

równań poza uprzednio 'wyjedynkowanym'. W algorytmie zastosowano zamianę równań w przypadku, gdy

współczynnik przy kolejnej zmiennej znaczącej jest równy zero.

5.1.

Kod:

Twynik Eliminacja_Gaussa_Jordana(Twsp A, Twynik B)
{

int

licznik,i,j;

int

k = 0;

Twynik X; X.n = A.m;

//Eliminacja zmiennych

for

(licznik=0; licznik<A.n; licznik++)

{

(A.t[k][k] == 0)

{

double

max = A.t[k][k];

int

wiersz = k;

double

temp;

for

(

int

p=k; p<A.n; p++)

(fabs(A.t[p][k]) > fabs(max))

{max = A.t[p][k]; wiersz = p;}

//zamiana równań

for

(

int

p=0; p<A.m; p++)

{

temp = A.t[k][p];
A.t[k][p] = A.t[wiersz][p];
A.t[wiersz][p] = temp;

}
temp = B.t[k]; B.t[k] = B.t[wiersz]; B.t[wiersz] = temp;

}

//wyjedynkowanie pierwszego współczynnika

double

temp = A.t[k][k];

for

(j=k; j<A.m; j++)

A.t[k][j] = A.t[k][j] / temp;

B.t[k] = B.t[k] / temp;

//Postępowanie Jordana

for

(i=0; i<A.n; i++)

(i != k)

{

double

temp = A.t[i][k];

for

(j=k; j<A.m; j++)

A.t[i][j] = A.t[i][j] - temp*A.t[k][j];

B.t[i] = B.t[i] - temp*B.t[k];

}

k++;

}

//Postępowanie odwrotne

for

(i=0; i<A.n; i++)

X.t[i] = B.t[i];

return

}

5/9

5.2.

Kroki algorytmu.

[

1 2 3
3 1 2
2 3 1

]

⇒

[

0 −5 −7
0 −1 −5

]

⇒

[

1 0

0.2

0 1

1.4

0 0 −3.6

]

⇒

[

1 0 0
0 1 0
0 0 1

]

X =

[

0.667 0.667 −0.333

]

6. Eliminacja rozkładem LU.

Metoda ta polega na znalezieniu macierzy trójkątnej dolnej L (której elementy na przekątnej są

równe 1) i macierzy trójkątnej górnej U, których iloczyn daje macierz współczynników układu równań

(Twsp A). Po ich uzyskaniu korzystamy z zależności:

{

L⋅y=B

U⋅x= y

gdzie y, x są typu (Twynik). Ponieważ macierze L i U są trójkątne, to do uzyskania y, x korzystamy z

'postępowania odwrotnego'.

6.1.

Kod:

Twynik Eliminacja_LU(Twsp A, Twynik B)
{

int

i,j,k;

double

suma;

Twynik X; X.n = A.m;
Twynik Y; Y.n = A.m;

double

tab[n_MAX][m_MAX] = {0};

Twsp L = load_wsp(tab, A.n, A.m);
Twsp U = load_wsp(tab, A.n, A.m);

//L: jedynki na przekątnej

k = 0;

for

(i=0; i<L.n; i++) {L.t[i][k] = 1; k++;}

//Rozkład LU

for

(i=0; i<L.n; i++)

{

for

(j=i; j<U.m; j++)

//kolejne kolumny

{

suma = 0;

for

(k=0; k<=j; k++)

suma += (L.t[i][k]*U.t[k][j]);

U.t[i][j] = (A.t[i][j] - suma);

}

for

(j=i+1; j<L.n; j++)

//kolejne wiersze

{

suma = 0;

for

(k=0; k<=j; k++)

suma += (L.t[j][k]*U.t[k][i]);

L.t[j][i] = (A.t[j][i] - suma)/U.t[i][i];

}

6/9

//LY = B. Obliczanie Y

for

(i=0; i<L.n; i++)

{

suma = 0;

for

(j=0; j<i; j++)

suma += L.t[i][j] * Y.t[j];

Y.t[i] = (B.t[i] - suma)/L.t[i][i];

}

//UX = Y. Obliczanie X

for

(i=(U.n - 1); i>=0; i--)

{

suma = 0;

for

(j=(i+1); j<U.m; j++)

suma += U.t[i][j] * X.t[j];

X.t[i] = (Y.t[i] - suma)/U.t[i][i];

}

return

}

6.2.

Kroki algorytmu.

[

1 0 0
0 1 0
0 0 1

]

⇒

[

1 0 0
3 1 0
2 0 1

]

⇒

[

2 0.2 1

]

U =

[

0 0 0
0 0 0
0 0 0

]

⇒

[

1 2 3
0 0 0
0 0 0

]

⇒

[

0 −5 −7
0

]

⇒

[

0 −5

−

3.6

]

X =

[

0.667 0.667 −0.333

]

7. Obliczanie wyznacznika macierzy.

Korzystając z rozkładu macierzy współczynników układu równań (Twsp A) na iloczyn macierzy

L i U, możemy zauważyć:

det  A=det  L⋅U =∣L∣⋅∣U∣

det  L=1 ⇒ det A=det U =

∏

n−1

gdzie i numerujemy od zera, a n to liczba wierszy macierzy U.

7.1.

Kod:

double

det(Twsp A)

{

int

i,j,k;

double

suma, temp=1;

double

tab[n_MAX][m_MAX] = {0};

Twsp L = load_wsp(tab, A.n, A.m);

7/9

Twsp U = load_wsp(tab, A.n, A.m);

//L: jedynki na przekątnej

k = 0;

for

(i=0; i<L.n; i++) {L.t[i][k] = 1; k++;}

//Rozkład LU

for

(i=0; i<L.n; i++)

{

for

(j=i; j<U.m; j++)

//kolejne kolumny

{

suma = 0;

for

(k=0; k<=j; k++)

suma += (L.t[i][k]*U.t[k][j]);

U.t[i][j] = (A.t[i][j] - suma);

}

for

(j=i+1; j<L.n; j++)

//kolejne wiersze

{

suma = 0;

for

(k=0; k<=j; k++)

suma += (L.t[j][k]*U.t[k][i]);

L.t[j][i] = (A.t[j][i] - suma)/U.t[i][i];

}

//det(A) = iloczyn el. przekątnej U

for

(i=0; i<U.n; i++) temp *= U.t[i][i];

return

temp;

}

7.2.

Kroki algorytmu.

[

1 0 0
0 1 0
0 0 1

]

⇒

[

1 0 0
3 1 0
2 0 1

]

⇒

[

2 0.2 1

]

U =

[

0 0 0
0 0 0
0 0 0

]

⇒

[

1 2 3
0 0 0
0 0 0

]

⇒

[

0 −5 −7
0

]

⇒

[

0 −5

−

3.6

]

det  A=

∣

1 2 3
3 1 2
2 3 1

∣

8. Wyznaczanie macierzy odwrotnej.

Z definicji iloczyn macierzy

i macierzy do niej odwrotnej

−

daje macierz jednostkową

, zatem traktując kolejną kolumnę macierzy odwrotnej jako wektor zmiennych, a kolejeną kolumnę

macierzy jednostkowej jako wektor stałych równań układu, możemy w kolejnych iteracjach obliczyć te

zmienne, czyli skłądowe macierzy odwrotnej. Do tego wykorzystujemy jedną z metod rozwiązywania układu

równań liniowych. Przy każdej iteracji musimy odświeżyć macierz A za pomocą metody Copy_Twsp,

ponieważ tablica w strukturze Twsp jest wskaźnikiem.

8/9

8.1.

Kod:

Twsp Macierz_Odwrotna(Twsp A)
{

int

i,j;

Twynik tmp;
Twynik x; x.n = A.n;

for

(i=0; i<x.n; i++) x.t[i] = 0;

double

tab[n_MAX][m_MAX] = {0};

Twsp A_original = load_wsp(tab, A.n, A.m);
Copy_Twsp(A, A_original);
Twsp A_ = load_wsp(tab, A.n, A.m);

//A*A_=I, to A*A_[][i]=I[][i]

for

(j=0; j<A.m; j++)

{

//kolumna macierzy jednostkowej

(j == 0) x.t[j] = 1;

else

{x.t[j-1] = 0; x.t[j] = 1;}

tmp = Eliminacja_Gaussa_Jordana(A, x);

for

(i=0; i<A_.n; i++) A_.t[i][j] = tmp.t[i];

Copy_Twsp(A_original, A);

//bo struktura Twsp, chociaż przekazywana

//przez wartość, to zawiera wskaźnik

}

return

A_;

}

8.2.

Kroki algorytmu.

−

[

0 0 0
0 0 0
0 0 0

]

⇒

[

−

0.278 0 0

0.0556 0 0

0.389 0 0

]

⇒

[

−

0.278

0.389 0

0.0556 −0.278 0

0.389

0.0556 0

]

⇒

[

−

0.278

0.389

0.0556

0.0556 −0.278

0.389

0.0556 −0.278

]

9/9