Sprawozdanie I

Metody numeryczne Rok akademicki

2008/2009

Termin

Czwartek 30

30

11 – 14

Rodzaj studiów

Kierunek

Prowadzący

Grupa

Sekcja

Dzienne

INF

dr inż. Ewa Starzewska-Karwan

1

2

Sprawozdanie z ćwiczenia numer 1 Data wykonania ćwiczenia: 2009-05-07

Temat ćwiczenia:

Układy równań liniowych

Skład sekcji: Gabryś Adam Łopuszański Łukasz

1. Kod programu #include #include #include #include void // bool void void



wczytaj_dane(int &n); zwraca true, jezeli udalo sie stworzyc macierz trojkatna tworz_macierz_trojkatna(const int &n); wyznacz_niewiadome(const int &n); raportuj(const int &n, bool &error);

using namespace std; // // // // // // //

rozmiar macierzy A i E oraz wektorow B i X const int rozmiar = 30; macierz zawierajaca wczytane wspolczynniki double macierzA_oryginal[rozmiar][rozmiar]; wektor zawierajacy wczytane wyrazy wolne double wektorB_oryginal[rozmiar]; wektor zawierajacy wyliczone wartosci niewiadomych double wektorX[rozmiar]; macierz zawierajaca wspolczynniki modyfikowane w programie double macierzA[rozmiar][rozmiar]; wektor zawierajacy wyrazy wolne modyfikowane w programie double wektorB[rozmiar]; wektor zawierajacy dokladne rozwiazania (jesli poda je uzytkownik) double wektorX_dokladne[rozmiar];

//

macierz odwrotna: A^(-1) (jest wyliczana jesli uzytkownik zazyczy sobie wyznaczenia wspolczynnikow uwarunkowania)

//

double macierzA_odwrotna[rozmiar][rozmiar]; macierz jednostkowa double macierzE[rozmiar][rozmiar];

int main(void) { // n - rozmiar macierzy wspolczynnikow oraz wektorow: wyrazow wolnych i z niewiadomymi int n; // zmienna przyjmuje wartosc false jesli nie wystapi blad obliczen (proba podzielenia przez 0) bool error = true; wczytaj_dane(n); if (tworz_macierz_trojkatna(n)) { wyznacz_niewiadome(n); error = false; } raportuj(n, error); return 0; } // funkcja wczytujaca dane z pliku do: // - macierzA i macierzA_oryginal // - wektorB i wektorB_oryginal // pobiera takze rozmiar macierzy A oraz wektorow B i X (przechowuje je zmienna n) // w przypadku podania nazwy nieistniejacego pliku wypisuje stosowny komunikat i probuje pobrac nazwe jeszcze raz

// funkcja zaklada, ze uzytkownik pwskazal plik z poprawnymi danymi void wczytaj_dane(int &n) { cout << " Podaj rozmiar macierzy kwadratowej A, wektora B i wektora X: "; cin >> n; string filename; cout << "\n Podaj nazwe pliku z danymi: "; cin >> filename; ifstream input(filename.c_str(), ios_base::in); while(!input.is_open()) { cout << "\n Plik o podanej nazwie nie istnieje!!!\n Podaj nazwe istniejacego pliku z danymi: ";

cin >> filename; input.clear(); input.open(filename.c_str(), ios_base::in); } for (int i = 0; i < n; ++i) { for (int j = 0; j < n; ++j) { input >> macierzA_oryginal[i][j]; macierzA[i][j] = macierzA_oryginal[i][j]; } input >> wektorB_oryginal[i]; wektorB[i] = wektorB_oryginal[i]; } input.close(); }

// funkcja tworzy macierz trojkatna i odpowiednio modyfikuje wektor B bool tworz_macierz_trojkatna(const int &n) { for(int k = 0; k < n-1; ++k) { for (int i = k+1; i < n; ++i) { wektorB[i] = wektorB[i] - (macierzA[i][k]*wektorB[k]) / macierzA[k][k]; for (int j = k+1; j < n ; ++j) { if (macierzA[k][k] == 0) return false; macierzA[i][j] = macierzA[i][j] - (macierzA[i][k]*macierzA[k][j]) / macierzA[k][k];

} macierzA[i][k] = 0; } } return true; } // funkcja wyznacza niewiadome z macierzy trojkatnej i zmodyfikowanego wektora B void wyznacz_niewiadome(const int &n) { double suma; // wyznaczenie pierwszej niewiadomej wektorX[n-1] = wektorB[n-1] / macierzA[n-1][n-1]; for (int i = n-2; i >= 0; --i) { suma = 0.0; for (int j = i+1 ; j < n; ++j) { suma += macierzA[i][j] * wektorX[j]; } wektorX[i] = (wektorB[i] - suma) / macierzA[i][i]; } } // funkcja zapisujaca raport z przeprowadzonych operacji do pliku // dodatko umozliwia generowanie raportu zgodnego z upodobaniami uzytkownika: // - mozna dodac wartosci wyliczonego bledu bezwzglednego (jezeli zostana wprowadzone wartosci dokladne) // - mozna dodac wspolczynniki uwarunkowania (w tym celu zostaje wyznaczona macierz A^(-1)) void raportuj(const int &n, bool &error) { string filename; cout << "\n Podaj nazwe pliku do ktorego zostanie zapisany raport: "; cin >> filename; ofstream output(filename.c_str(), ios_base::out); if (error) { output << "Blad - obliczenia przerwano!!!\nPowod: niedozwolona operacja arytmetyczna (dzielenie przez 0)\n";

return; } //

jezeli uzytkownik zna dokladne rozwiazania, to w raporcie zostana dodane wartosci bledu bezwzglednego char bledy; cout << " Czy znasz dokladne rozwiazania rownania (t - tak, n - nie)?: "; cin >> bledy; if (bledy == 't') { cout << " Podaj dokladne rozwiazania:\n"; for (int i = 0; i < n; ++i) { cout << " x[" << i+1 << "] = "; cin >> wektorX_dokladne[i]; } } char kryterium; cout << " Czy chcesz sprawdzic kryterium zlego uwarunkowania (t - tak, n - nie)?: "; cin >> kryterium;

//

ustawienie notacji naukowej oraz precyzji mantysy na 10 pozycji output.setf(ios_base::scientific); output.precision(10); output << "Macierz wspolczynnikow A i wektor wyrazow wolnych B:\n"; for (int i = 0; i < n; ++i) { for (int j = 0; j < n; ++j) { output << macierzA_oryginal[i][j] << " "; } output << " \t" << wektorB_oryginal[i] << endl; }

output << "\n\nMacierz trojkatna i wektor zmodyfikowanych wyrazow wolnych:\n"; for (int i = 0; i < n; ++i) { for (int j = 0; j < n; ++j) { output << macierzA[i][j] << " "; } \t" << wektorB[i] << endl; output << " } //

jezeli zostaly podane wartosci dokladne, to zostana dolaczone wartosci bledow bezwzglednych output << "\n\nRozwiazanie ukladu"; if (bledy == 't') { double temp; output << " i blad bezwzgledny:\n"; for (int i = 0; i < n; ++i) { temp = wektorX[i] - wektorX_dokladne[i]; output << "x[" << i+1 << "] = " << wektorX[i] << "\tblad: " << (temp >= 0 ? temp : -temp) << endl;

} } else { output << ":\n"; for (int i = 0; i < n; ++i) output << "x[" << i+1 << "] = " << wektorX[i] << endl; } // //

obliczanie kryterium zlego uwarunkowania (na zyczenie uzytkownika) if (kryterium == 't') { wypelnianie macierzy jednostkowej for (int i = 0; i < n; ++i) { for (int j = 0; j < n; ++j) { if (i == j) macierzE[i][j] = 1; else macierzE[i][j] = 0; } } for (int i = 0; i < n; ++i) {

//

przypisanie macierzA oryginalnch wartosci (potrzebne do wyznaczania macierzy odwrotnej)

for(int j = 0; j < n; ++j) for(int k = 0; k < n; ++k) macierzA[j][k] = macierzA_oryginal[j][k]; //

wypelnianie wektora B

//

k-ta kolumna macierzy E to wektor wyrazow wolnych, z niego liczymy k-ta kolumne macierzy odwroconej

for (int j = 0; j < n; ++j) wektorB[j] = macierzE[j][i]; rozwiazanie ukladu rownan tworz_macierz_trojkatna(n); wyznacz_niewiadome(n); for (int j = 0; j < n; ++j) macierzA_odwrotna[j][i] = wektorX[j];

//

} //

wypisanie do pliku macierzy odwrotnej output << "\n\nMacierz odwotna A^-1:\n"; for (int i = 0; i < n; ++i) { for (int j = 0; j < n; ++j) { output << macierzA_odwrotna[i][j] << " } output << endl; }

";

//

wyznaczanie norm macierzy A i A^(-1) wykorzystywanych do obliczenia wspolczynnika uwarunkowania

//

wybralismy normy I double normaA_I = 0; double normaA_II = 0; double normaA_III = 0; double normaA_odwrotna_I = 0; double normaA_odwrotna_II = 0; double normaA_odwrotna_III = 0; double temp = 0; double temp_odwrotna = 0;

//

obliczanie norm I for (int i = 0; i < n; ++i) { for (int j = 0; j < n; ++j) { temp += (macierzA_oryginal[i][j] >= 0 ? macierzA_oryginal[i][j] : -macierzA_oryginal[i][j]); temp_odwrotna += (macierzA_odwrotna[i][j] >= 0 ? macierzA_odwrotna[i][j] : -macierzA_odwrotna[i][j]);

} if (temp > normaA_I) normaA_I = temp; if (temp_odwrotna > normaA_odwrotna_I) normaA_odwrotna_I = temp_odwrotna; temp = 0; temp_odwrotna = 0; //

} obliczanie norm II for (int j = 0; j < n; ++j) { for (int i = 0; i < n; ++i) { temp += (macierzA_oryginal[i][j] >= 0 ? macierzA_oryginal[i][j] : -macierzA_oryginal[i][j]); temp_odwrotna += (macierzA_odwrotna[i][j] >= 0 ? macierzA_odwrotna[i][j] : -macierzA_odwrotna[i][j]);

} if (temp > normaA_II) normaA_II = temp; if (temp_odwrotna > normaA_odwrotna_II) normaA_odwrotna_II = temp_odwrotna; temp = 0; temp_odwrotna = 0; //

} obliczanie norm III for (int i = 0; i < n; ++i) for (int j = 0; j < n; ++j) { normaA_III += pow(macierzA_oryginal[i][j], 2); normaA_odwrotna_III += pow(macierzA_odwrotna[i][j], 2); } normaA_III = sqrt(normaA_III); normaA_odwrotna_III = sqrt(normaA_odwrotna_III); double double double output output output output output output output output

wynik_I = normaA_I * normaA_odwrotna_I; wynik_II = normaA_II * normaA_odwrotna_II; wynik_III = normaA_III * normaA_odwrotna_III; << "\n\nNormy macierzy A:\n"; << "Norma I: " << normaA_I << "\nNorma II: " << normaA_II; << "\nNorma III: " << normaA_III << endl; << "\n\nNormy macierzy odwrotnej A^-1:\n"; << "Norma I: " << normaA_odwrotna_I << "\nNorma II: " << normaA_odwrotna_II; << "\nNorma III: " << normaA_odwrotna_III << endl; << "\n\nWspółczynniki uwarunkowania:\n"; << "I: " << wynik_I << "\nII: " << wynik_II << "\nIII: " << wynik_III;

} output.close(); }

*pliki projektu zamieściliśmy pod adresem http://www.adamgabrys.ovh.org/mn/url.zip

2. Wprowadzanie danych oraz wygenerowanie raportów Użytkownik podczas korzystania z programu jest pytany o kilka ważnych informacji niezbędnych do przeprowadzenia obliczeń: • Wielkość macierzy współczynników, wektora wyrazów wolnych i wektora niewiadomych (gdzie Anxn, Xnx1 i Bnx1 - n to wielkość podana przez użytkownika). • Nazwę pliku, z którego zostaną pobrane dane. • Nazwę pliku, do którego zostanie zapisany raport z przeprowadzonych obliczeń. • Informację, czy znane są wartości dokładne niewiadomych (jeśli tak, to użytkownik zostanie poproszony o ich wprowadzenie). • Informację, czy mają zostać wyliczone współczynniki uwarunkowania. Po podaniu tych wszystkich informacji na dysku zostanie utworzony plik z raportem. Poniżej przedstawiamy przykładowy screen z pełnym przebiegiem wprowadzania danych:

Plik z danymi wejściowymi ma następującą postać (zestaw 1):

Liczby są oddzielone od siebie pojedynczymi spacjami. Ostatnia kolumna stanowi kolumnę wyrazów wolnych (czyli wektor B).

Otrzymany raport (dla zestawu 1) wygląda następująco:

W przypadku zażądania sprawdzenia kryterium złego uwarunkowania macierzy, raport jest wzbogacony o dodatkowe dane:

Program może także zwrócić raport zawierający informację o przerwaniu obliczeń. Jeżeli taka sytuacja wystąpi aplikacja nie wyświetla zapytań o dotyczących wartości dokładnych niewiadomych i wyznaczania współczynnika uwarunkowania:

Raport dla zestawu 4:

Raport dla zestawu 5:

*pliki z danymi wejściowymi i raporty zamieściliśmy pod adresem http://www.adamgabrys.ovh.org/mn/raporty.zip

3. Sprawdzenie kryterium złego uwarunkowania dla zestawu 2. W celu sprawdzenia kryterium złego uwarunkowania należy wyznaczyć macierz odwrotną do macierzy współczynników, oraz normy obu macierzy. Macierz będzie źle uwarunkowana, jeżeli przynajmniej jedna z niżej podanych nierówności będzie prawdziwa: ||A||I * ||A-1||I >> 1 ||A||II * ||A-1||II >> 1 ||A||III * ||A-1||III >> 1 Macierz A-1 oraz normy i ich iloczyny wyznacza program na żądanie użytkownika (chęć wyznaczenia kryterium złego uwarunkowania):

Obliczenia przeprowadzimy także ręcznie w celu zaprezentowania metody. Najpierw należy wyznaczyć macierz odwrotną (np. za pomocą programu MS Excel):

Następnie wyznaczamy normy dla obu macierzy z następujących wzorów:

Dla macierzy A:

suma[i=1, j] = |10,0001| + |2| + |3| + |1| + |4| = 20,0001 suma[i=2, j] = |13| + |20| + |2| + |2| + |3| = 40 suma[i=3, j] = |5| + |11| + |40| + |9| + |15| = 80 suma[i=4, j] = |3| + |2| + |5| + |10| + |0| = 20 suma[i=5, j] = |10| + |2| + |3| + |1| + |4| = 20 ||A||I = max suma[const i, j] = suma[i=3, j] = 80 suma[i, j=1] = |10,0001| + |13| + |5| + |3| + |10| = 41,0001 suma[i, j=2] = |2| + |20| + |11| + |2| + |2| = 37 suma[i, j=3] = |3| + |2| + |40| + |5| + |3| = 53 suma[i, j=4] = |1| + |2| + |9| + |10| +|1| = 23 suma[i, j=5] = |4| + |3| + |15| + |0| + |4| = 26 ||A||II = max suma[i, const j] = suma[i, j=3] = 53

10,0001ଶ + 2ଶ + 3ଶ + 1ଶ + 4ଶ + 13ଶ + 20ଶ + 2ଶ + 2ଶ + 3ଶ + 5ଶ + 11ଶ + ||‫||ܣ‬ூூூ = ට = 55,0999274 +40ଶ + 9ଶ + 15ଶ + 3ଶ + 2ଶ + 5ଶ + 10ଶ + 0ଶ + 10ଶ + 2ଶ + 3ଶ + 1ଶ + 4ଶ Dla macierzy A-1:

suma[i=1, j] = |10000| + |1,94833E-13| + |-4,65984E-14| + |-5,54039E-14| + |-10000| = 20000 suma[i=2, j] = |-2173,095996| + |0,054549215| + |0,001827008| + |-0,007777836| + |2173,048233| = = 4346,208 suma[i=3, j] = |13246,85494| + |-0,005115624| + |0,038106175| + |-0,019366289| + |-13246,994| = = 26493,91 suma[i=4, j] = |-9188,808269| + |-0,008352038| + |-0,019418489| + |0,111238712| + |9188,887352| = = 18377,83 suma[i=5, j] = |-31551,39114| + |-0,021349898| + |-0,024638513| + |-0,009396043| + |31551,74954| = = 63103,2 ||A-1||I = max suma[const i, j] = suma[i=5,j] = 63103,2

suma[i, j=1] = |10000| + |-2173,095996| + |13246,85494| + |-9188,808269| + |-31551,39114| = 66160,15034 suma[i, j=2] = |1,94833E-13| + |0,054549251| + |-0,005115624| + |-0,008352038| + |-0,021349898| = 0,089366811 suma[i, j=3] = |-4,65984E-14| + |0,001827008| + |0,038106175| + |-0,019418189| + |-0,024638513| = 0,083990186 suma[i, j=4] = |-5,54039E-14| + |-0,007777836| + |-0,019366289| + |0,111238712| + |-0,009396043| = 0,14777888 suma[i, j=5] = |-10000| + |2173,048233| + |-13246,994| + |9188,887352| + |31551,74954| = 66160,67913 ||A-1||II = max suma[i, const j] = suma[i, j=5] = 66160,67913

||‫ିܣ‬ଵ ||ூூூ

ለ 10000ଶ + ሺ1,94833 ∗ 10ିଵଷ ሻଶ + ሺ−4,65984 ∗ 10ିଵସ ሻଶ + ሺ−5,54039 ∗ 10ିଵସ ሻଶ + ሺ−10000ሻଶ + ള ള ଶ + 0,054549251ଶ + 0,001827008ଶ + ሺ−0,007777836ሻଶ + 2173,048233ଶ + ള +ሺ−2173,095996ሻ ള ള ଶ + ሺ−0,005115624ሻଶ + 0,038106175ଶ + ሺ−0,019366289ሻଶ + +13246,85494 ള =ള = 52156,29905 ള ଶ +ሺ−13246,994ሻ + ሺ−9188,808269ሻଶ + ሺ−0,008352038ሻଶ + ሺ−0,019418489ሻଶ + ള ള ള +0,111238712ଶ + 9188,887352ଶ + ሺ−31551,39114ሻଶ + ሺ−0,021349898ሻଶ + ള +ሺ−0,024638513ሻଶ + ሺ−0,009396043ሻଶ + 31551,74954ଶ ۣ

Po wyznaczeniu norm obliczamy wartości współczynników: ||A||I * ||A-1||I = 5048255,685 ||A||II * ||A-1||II = 3506515,994 ||A||III * ||A-1||III = 2873808,291

Obliczone wartości są zbliżone do zapisanych w raporcie (różnice na mniej znaczących pozycjach spowodowane operacjami zaokrąglania). Wszystkie trzy współczynniki są znacznie dużo dużo większe od jedynki, co oznacza, że macierz jest źle uwarunkowana i małe zmiany w wartościach współczynników powodują duże błędy.

4. Wnioski Po przeanalizowaniu wygenerowanych raportów doszliśmy do następujących wniosków: •

•

•

•

Maszyna cyfrowa nie jest idealnym narzędziem do przeprowadzania obliczeń matematycznych. Powodem tego są ograniczone możliwości technologiczne, które można zauważyć w raporcie wygenerowanym dla zestawu 1 (dwa błędy bezwzględne są różne od 0, co jest spowodowane skończoną dokładnością liczb typu double – błędy zaokrąglania) Metoda Gaussa nie jest idealnym sposobem rozwiązywania układów równań liniowych. Nie radzi sobie z pewnymi kombinacjami rozłożenia współczynników, powodują próbę wykonania dzielenia przez 0 (zestaw 3) lub bardzo małe liczby (pojawiają się duże błędy zaokrąglenia). Aby tego uniknąć stosuje się metodę wyboru elementu podstawowego, która w znacznym stopniu rozwiązuje ten problem, ale znowu z kolei zwiększa zapotrzebowanie na pamięć. Wartości współczynników mają ogromne znaczenie na dokładność otrzymanego rozwiązania. Widać to dokładnie porównując raporty dla zestawu 4 i 5. Zmiany na pozycjach części dziesięciotysięcznych spowodowały, że błędy bezwzględne rzędu jedności, czy części dziesiętnych, zmalały do wartości rzędu 10-16. Zestawy 2, 4 i 5 mają bardzo wysokie wartości współczynników uwarunkowania. Oznacza to, że są to macierze źle uwarunkowane, czyli małe zmiany wartości współczynników powodują duże błędy w otrzymanych wynikach. Widać to dokładnie po przeanalizowaniu wartości błędów bezwzględnych: dla zestawu 4 największy błąd wyniósł ok. 91, a dla zestawu 2 ok. 2,8 * 10-11. Należy zwrócić uwagę, że różnica między tymi zestawami polega jedynie na zmianie wartości jednego współczynnika w wektorze B z 20 na 20,0001.

Dodatkowo zbadaliśmy zestaw 1 pod kątem uwarunkowania macierzy:

Jak widać, wartości współczynników są małe, czyli macierz nie jest źle uwarunkowana.