Politechnika Śląska Rok akademicki 2007/2008 semestr zimowy

Laboratorium Metod Numerycznych

ROZWIĄZYWANIE UKŁADÓW RÓWNAŃ LINIOWCYH

TEMAT LABORATORIUM : METODA CHOLESKIEGO ELIMINACJĄ GAUSSA

kier. Informatyka.

Grupa I

Sekcja 2

Łukasz Głomb

Wstęp.

Celem ćwiczenia było rozwiązanie układu typu A X = B metodą Choleskiego (eliminacja Gaussa).

Stosując tą metodę należy przekształcić macierz A na macierze L i U, by zachodził następujący związek:

A=LU

gdzie L = U = .

Algorytm rozwiązania tego układu jest następujący:

Do równania AX = B za A podstawiamy LU i otrzymujemy LUX = B,
Wcześniej otrzymane równanie przekształcamy na następujące: LY=B i UX=Y,
Z pierwszego równania wyznaczamy wektor Y a z drugiego wektor X oraz rozwiązujemy następujące układy równań:

Macierze L i U otrzymujemy przy kolejnych etapach eliminacji Gaussa macierzy A i obliczamy je według następującego wzorów:

Kod źródłowy programu:

#include <stdio.h>

#include <math.h>

#include <stdlib.h>

#include <string.h>

double **A = NULL; //Macierze A i wektor B

double *B = NULL; // A - macierz wspolczynnikow, B - wektor wyrazow wolnych

double **L = NULL; //Macierze L i U powstale po przeksztalceniu

double **U = NULL; //macierzy A

double **A2 = NULL; //Pomocnicza macierz (uzywany jej przy elminacji Gaussa)

int N; // Ilosc rownan

double *X = NULL; //Wektory X i Y stanowia

double *Y = NULL; //rozwiazanie ukladu

// Funkcja a:

// - wczytuje dane z pliku zrodlowego

// - oznaczenia:

// N - ilosc rownan

void a(N)

{

FILE *plik;

double tmp;

int i,

plik = fopen("macierz.txt","r");

// Poczatek rezerowania pamieci na macierz A i wektor B

A = (double**) malloc ((N+1)*sizeof(double*));

for (i=0; i<=N; i++)

{

A[i] = (double*)malloc((N+1)*sizeof(double));

}

B = (double*) malloc ((N+1)*sizeof(double));

// koniec rezerwowania pamieci

for (i=1; i<=N; i++) //zczytywanie watrości z pliku

for (j=1; j<= N+1; j++)

{

fscanf(plik, "%lf", &tmp);

if (j<=N)

A[i][j] = tmp;

else

B[i] = tmp;

}

fclose(plik);

return;

}

// Funkcja b tworzy macierze L i U metodą Gaussa oraz

// sprawdza warunek utworzenia tych macierzy

int b(N)

{

int i,

L = (double**) malloc ((N+1)*sizeof(double*));

for (i=0; i<=N; i++)

{

L[i] = (double*)malloc((N+1)*sizeof(double));

//rezerwoanie pamięci na macierze

} // L i U

U = (double**) malloc ((N+1)*sizeof(double*));

for (i=0; i<=N; i++)

{

U[i] = (double*)malloc((N+1)*sizeof(double));

}

A2 = (double**) malloc ((N+1)*sizeof(double*)); //rezerwoanie pamięci

for (i=0; i<=N; i++) //na macierz pomocną

{

A2[i] = (double*)malloc((N+1)*sizeof(double));

}

printf("Zarezerwowano pamięć\n");

for (i=1; i<=N; i++)

for (j=1; j<=N; j++)

{

L[i][j] = 0; //zerowanie macierzy L i U

U[i][j] = 0;

A2[i][j] = A[i][j]; //kopiowanie macierzy A do macierzy

} //pomocniczej A2

for (i=1; i<=N; i++) //wstawienie jedynek po przekatnej

{ //macierzy L

L[i][i] = 1;

}

for (i=1; i<=N; ++i)

{

for (j=1; j<=N; ++j)

printf("A2[%d][%d] = %f ",i, j, A2[i][j]);

printf("\n");

}

printf("\n\n\n");

for (k=1; k<N; k++)

{

if ((A2[k][k]) == 0)

{

printf("Na przekotnej wystapilo zero\n");

printf("Nie został utworzony rapot");

exit(1); //Zakonczenie pracy programu

}

for (i=k+1; i<=N; ++i)

{

for (j=k+1; j<=N; ++j)

{

A2[i][j]= (A2[i][j]) - (A2[k][j])*((A2[i][k])/A2[k][k]);

}

L[i][k] = A2[i][k]/A2[k][k];

A2[i][k] = 0;

}

for (k = 1; k<=N; k++)

for (i = k; i<=N; i++)

{

U[k][i]= A2[k][i];

}

for (i=1; i<=N; ++i)

{

for (j=1; j<=N; ++j)

printf("L[%d][%d] = %f ",i, j, L[i][j]);

printf("\n");

}

printf("\n\n");

for (i=1; i<=N; ++i)

{

for (j=1; j<=N; ++j)

printf("U[%d][%d] = %f ",i, j, U[i][j]);

printf("\n");

}

return 0;

}

// Funkcja c oblicza wektory X i Y

void c()

{

double suma1,

suma2;

int i,k;

Y = (double*) malloc ((N+1)*sizeof(double));

X = (double*) malloc ((N+1)*sizeof(double));

Y[1]=B[1];

for (i=2; i<=N; i++)

{

suma1=0;

for (k=1; k<=i-1; k++)

{

suma1 += L[i][k]* Y[k] ;

}

Y[i] = B[i]-suma1;

}

X[N] = Y[N] / U[N][N];

for (i=N-1; i>=1; i--)

{

suma2=0;

for (k=i+1; k <= N; k++)

{

suma2 += U[i][k]*X[k];

}

X[i] = (Y[i] - suma2) / U[i][i];

}

//Funkcja d tworzy raport do pliku, wektory X i Y w postaci wykładniczej

//mantysa ma 10 cyfr znaczacych

void d()

{

FILE *plik;

int i,j;

char zapis_plik[20];

printf("Podaj nazwe pliku z rozszerzeniem do zapisu");

scanf("%s",&zapis_plik);

plik = fopen (zapis_plik, "w");

fprintf(plik, "Macierz A:\n");

for (i=1; i<=N; i++)

{

for (j=1; j<=N; j++)

{

fprintf(plik,"%f ", A[i][j]);

}

fprintf(plik,"\n");

}

fprintf(plik, "Wektor B:\n");

for (i=1; i<=N; i++)

{

fprintf(plik,"%f \n", B[i]);

}

fprintf(plik, "Macierz L:\n");

for (i=1; i<=N; i++)

{

for (j=1; j<=N; j++)

{

fprintf(plik,"%f ", L[i][j]);

}

fprintf(plik,"\n");

}

fprintf(plik, "Macierz U:\n");

for (i=1; i<=N; i++)

{

for (j=1; j<=N; j++)

{

fprintf(plik,"%f ", U[i][j]);

}

fprintf(plik,"\n");

}

fprintf(plik, "Wektor X:\n");

for (i=1; i<=N; i++)

{

fprintf(plik,"%.10f\n", X[i]);

}

fprintf(plik, "Wektor Y:\n");

for (i=1; i<=N; i++)

{

fprintf(plik,"%.10e\n", Y[i]);

}

fclose(plik);

}

int main()

{

printf("Podaj liczbe ukladow rownan:\n");

scanf("%d",&N);

a(N);

b(N);

c(N);

d(N);

return 0;

}

Dane wejściowe.

Plik z danymi wejściowymi powinien znajdować się w katalogu razem z programem. Jego nazwa powinna być: macierz.txt . Użytkownik podaję na początku ilość, rozmiar macierzy A oraz wektora B. Dane w pliku powinny być oddzielone spacjami oraz mieć postać:

a₁₁ a₁₂ ... a_1n b₁

a₂₁ a₂₂ … a_2n b₂

. . .

a_n1 a_n2 … a_nn b_n

Dane wyjściowe.

Jeżeli udało się utworzyć macierze L i U, a co za tym idzie wyznaczyć wektory X i Y to jest tworzony raport do pliku wskazanego przez użytkownika programu (ścieżka do pliku nie może przekroczyć 20 znaków). Pozostałe dane są zapisywane w typie double.

Raport ma postać:

Macierz A:

//Zawartość macierzy A

Wektor B:

//Zawartość wektora B

Macierz L:

//Zawartość macierzy L

Macierz U:

//Zawartość macierzy U

Wektor X:

//Zawartość wektora X

Wektor Y:

//Zawartość wektora Y

Raporty.

Dla wszystkich zestawów (oprócz 4.) rozwiązaniem dokładnym był wektor X=[1 1 1 1 1]^T, błędy bezwzględne zostały obliczone według wzoru:

, gdzie X – rozwiązanie dokładne, x- rozwiązanie jakie otrzymaliśmy

Zestaw 1:

Macierz A:

10.000000 1.000000 1.000000 1.000000 2.000000

3.000000 20.000000 4.000000 2.000000 1.000000

5.000000 1.000000 40.000000 9.000000 5.000000

3.000000 0.000000 1.000000 10.000000 1.000000

6.000000 2.000000 1.000000 1.000000 20.000000

Wektor B:

15.000000

30.000000

60.000000

15.000000

30.000000

Macierz L:

1.000000 0.000000 0.000000 0.000000 0.000000

0.300000 1.000000 0.000000 0.000000 0.000000

0.500000 0.025381 1.000000 0.000000 0.000000

0.300000 -0.015228 0.019194 1.000000 0.000000

0.600000 0.071066 0.003478 0.026117 1.000000

Macierz U:

10.000000 1.000000 1.000000 1.000000 2.000000

0.000000 19.700000 3.700000 1.700000 0.400000

0.000000 0.000000 39.406091 8.456853 3.989848

0.000000 0.000000 0.000000 9.563571 0.329512

0.000000 0.000000 0.000000 0.000000 18.749091

Wektor X:

1.0000000000e+000

Wektor Y:

1.5000000000e+001

2.5500000000e+001

5.1852791878e+001

9.8930825712e+000

1.8749090811e+001

Jak możemy zaobserwować z powyższego rozwiązania błąd bezwzględny dla każdego X jest równy 0. Możliwe, że jest błąd jest bardzo mały, lecz z powodu na dokładność (mantysa ma 10 cyfr znaczących) program nie zanotował żadnych odchyleń

Zestaw 2:

Macierz A:

10.000100 2.000000 3.000000 1.000000 4.000000

13.000000 20.000000 2.000000 2.000000 3.000000

5.000000 11.000000 40.000000 9.000000 15.000000

3.000000 2.000000 5.000000 10.000000 0.000000

10.000000 2.000000 3.000000 1.000000 4.000000

Wektor B:

20.000100

40.000000

80.000000

20.000000

Macierz L:

1.000000 0.000000 0.000000 0.000000 0.000000

1.299987 1.000000 0.000000 0.000000 0.000000

0.499995 0.574712 1.000000 0.000000 0.000000

0.299997 0.080460 0.107418 1.000000 0.000000

0.999990 0.000001 0.000001 0.000000 1.000000

Macierz U:

10.000100 2.000000 3.000000 1.000000 4.000000

0.000000 17.400026 -1.899961 0.700013 -2.199948

0.000000 0.000000 39.591946 8.097699 14.264357

0.000000 0.000000 0.000000 8.773843 -2.555226

0.000000 0.000000 0.000000 0.000000 0.000032

Wektor X:

1.0000000000e+000

9.9999999999e-001

9.9999999996e-001

Wektor Y:

2.0000100000e+001

1.4000129999e+001

6.1954002253e+001

6.2186168227e+000

3.1693963550e-005

W tym przypadku błąd bezwzględny dla X₄ i X₅ jest bardzo mały i jest rzędu . Dla pozostałych przypadków również nie jesteśmy wstanie zaobserwować różnicy po między oryginalnym rozwiązaniem.

Dla tego zestawu musieliśmy jeszcze sprawdzić kryterium złego uwarunkowania macierzy ( to macierz jest źle uwarunkowana). W naszym przypadku ten współczynnik wynosi około , więc ten zestaw jest źle uwarunkowany. Po mimo tego wyniki wyszły bardzo dokładnie.

Zestaw 3:

Nie został wygenerowany raport do pliku, ponieważ na przekątnej wystąpiło 0.

Zestaw 4:

Macierz A:

10.000100 2.000000 3.000000 1.000000 4.000000

13.000000 20.000000 2.000000 2.000000 3.000000

5.000000 11.000000 40.000000 9.000000 15.000000

3.000000 2.000000 5.000000 10.000000 0.000000

10.000000 2.000000 3.000000 1.000000 4.000000

Wektor B:

20.000100

40.000000

80.000000

20.000000

20.000100

Macierz L:

1.000000 0.000000 0.000000 0.000000 0.000000

1.299987 1.000000 0.000000 0.000000 0.000000

0.499995 0.574712 1.000000 0.000000 0.000000

0.299997 0.080460 0.107418 1.000000 0.000000

0.999990 0.000001 0.000001 0.000000 1.000000

Macierz U:

10.000100 2.000000 3.000000 1.000000 4.000000

0.000000 17.400026 -1.899961 0.700013 -2.199948

0.000000 0.000000 39.591946 8.097699 14.264357

0.000000 0.000000 0.000000 8.773843 -2.555226

0.000000 0.000000 0.000000 0.000000 0.000032

Wektor X:

1.0768522747e-011

1.2173048233e+000

-3.2469939968e-001

1.9188887352e+000

4.1551749543e+000

Wektor Y:

2.0000100000e+001

1.4000129999e+001

6.1954002253e+001

6.2186168227e+000

1.3169396355e-004

W tym zestawie dokładnie rozwiązanie nie ma postaci X=[1 1 1 1 1]^T, ponieważ został zaburzony ostatni element wektora B. Nie jesteśmy w stanie policzyć dokładnego błędu bezwzględnego, ponieważ nie znamy dokładnych rozwiązań.

Zestaw 5:

Macierz A:

10.000100 2.000000 3.000000 1.000000 4.000000

13.000000 20.000000 2.000000 2.000000 3.000000

5.000000 11.000000 40.000000 9.000000 15.000000

3.000000 2.000000 5.000000 10.000000 0.000000

10.000000 2.000100 3.000100 1.000000 4.000000

Wektor B:

20.000100

40.000000

80.000000

20.000000

20.000200

Macierz L:

1.000000 0.000000 0.000000 0.000000 0.000000

1.299987 1.000000 0.000000 0.000000 0.000000

0.499995 0.574712 1.000000 0.000000 0.000000

0.299997 0.080460 0.107418 1.000000 0.000000

0.999990 0.000007 0.000004 -0.000003 1.000000

Macierz U:

10.000100 2.000000 3.000000 1.000000 4.000000

0.000000 17.400026 -1.899961 0.700013 -2.199948

0.000000 0.000000 39.591946 8.097699 14.264357

0.000000 0.000000 0.000000 8.773843 -2.555226

0.000000 0.000000 0.000000 0.000000 -0.000003

Wektor X:

9.9999999988e-001

1.0000000000e+000

9.9999999985e-001

1.0000000001e+000

1.0000000004e+000

Wektor Y:

2.0000100000e+001

1.4000129999e+001

6.1954002253e+001

6.2186168227e+000

-3.4037597914e-006

Błędy bezwzględne dla tego układu wynoszą:

, , , ,

Wnioski

Na podstawnie wyników nasuwają się następujące wnioski:

metoda Choleskiego (macierze L i U otrzymywane z eliminacji Gaussa) jest metodą bardzo dokładną, błędy bezwzględne są bardzo małe, prawie nie wpływające na otrzymany wynik
złe uwarunkowanie macierzy nie wpływa znacząco na poprawność otrzymywanych wyników (przykładem zestaw 2.)
jeżeli nie zostania zastosowana metoda wyboru elementu podstawowego istnieje możliwość nierozwiązana układu po mimo istnienia rozwiązania macierzy (przykładem jest zestaw 3.). By uniknąć tego problemu należałoby zaopatrzyć program w algorytm wyboru elementu podstawowego
nawet małe zaburzenie dowolnego elementu może wpłynąć znacząco na otrzymane wyniki (przykładem zestaw 4.)