Politechnika Śląska Rok akademicki 2007/2008 semestr zimowy
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:
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,
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;
}
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
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
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ą:
, , , ,
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.)