Układy równań liniowych sprawozdanie poprawione

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




  1. 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:




  1. 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,

j;


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,

j,

k;


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;

}



  1. 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ć:


a11 a12 ... a1n b1

a21 a22 a2n b2

. . .

. . .

. . .

an1 an2 ann bn




  1. 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




  1. 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

1.0000000000e+000

1.0000000000e+000

1.0000000000e+000

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

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

1.0000000000e+000

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 X4 i X5 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ą:

, , , ,







  1. Wnioski


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


12



Wyszukiwarka

Podobne podstrony:
Zestaw 12 Macierz odwrotna, układy równań liniowych
lab8 1 uklady rownan liniowych
Układy równań liniowych
2011 lab 02, Uklady rownan liniowych
Układy równań liniowych
układy równań liniowych 2
Układy równań liniowych z parametrem
Matematyka I (Ćw) Lista 05 Układy m równań liniowych z n niewiadomymi
Układy równań liniowych, Matematyka dla ekonomistów
Uklady rownan liniowych
02. Układy równań liniowych
2011 lab 02 Uklady rownan liniowychid 27450
02 Układy równań liniowychid 3448
Zestaw uklady rownan liniowych
Układy równań liniowych z trzema niewiadomymi
Układy równań liniowych
matematyka, Układy równań liniowych, Układy równań liniowych o dwóch niewiadomych
6-MACIERZE, WYZNACZNIKI, UKŁADY RÓWNAŃ LINIOWYCH, MACIERZE I WYZNACZNIKI
W2 RZAD MACIERZY UKLADY ROWNAN LINIOWYCH, UEP lata 2014-2019, Ekonometria

więcej podobnych podstron