W tym poście, zostaną omówione jedne z najczęściej spotykanych metod numerycznych dzięki którym można rozwiązać układ n równań liniowych z n niewiadomymi. Zostaną omówione metody:
- Eliminacja Gaussa
- Eliminacja Gaussa - Jordana
- LU
- Metoda Jacobiego
- Gauss - Seidel
1. Eliminacja Gaussa
Na początek zaczniemy od krótkiego przykładu, który pokaże w jaki sposób obliczyć dany układ tą metodą.
Mamy następujący układ równań:
Zapiszemy powyższe równanie w postaci macierzy:
Mając tak zapisaną macierz przystępujemy do eliminacji. Naszym celem jest uzyskanie macierzy:
Gdzie x jest wartością którą obliczymy. Aby doprowadzić macierz do takiej postaci najpierw eliminujemy x1 z równania drugiego. W tym celu do równania drugiego musimy dodać pierwsze pomnożone przez -2/3:
Następnie z równania trzeciego eliminujemy także x1 poprzez pomnożenie równania pierwszego przez -1/3 i dodanie do równania trzeciego:
Teraz czas na eliminację z równania trzeciego x2. Dokonujemy tego poprzez pomnożenie równania drugiego przez -1/14:
W tym momencie posiadamy tzw. Macierz trójkątną górną. Postępowanie prowadzimy od końca najpierw obliczając x3 następnie wprowadzając do równania drugiego otrzymaną wartość i obliczenie x2 itd.
Ostatecznie wynik jaki otrzymamy będzie następujący:
Przejdźmy teraz do obliczeń od strony programistycznej.
W pierwszym kroku z równania drugiego eliminowaliśmy x1. Aby tego dokonać obliczyliśmy przez co należy pomnożyć równanie pierwsze, aby po odjęciu zlikwidowało się x1 w równaniu drugim. Szukaną wartością było -a_21/a11, otrzymaną wartość mnożyliśmy przez wszystkie wyrazy równania pierwszego i dodawaliśmy bezpośrednio do drugiego:
W ogólności dla drugiego równania mamy:
W kolejnych krokach eliminujemy z równania trzeciego x1 a następnie z równania trzeciego x2. W ogólności mamy więc n – 1 kroków. Ogólnie nasz wzór na obliczenie macierzy trójkątnej górnej wygląda następująco:
Postępowanie odwrotne:
Z ostatniego równania obliczamy x3:
Następnie otrzymaną wartość wstawiamy do równania drugiego i obliczamy x2. W ogólności wzór na obliczenie xn jest następujący:
Przykładowa implementacja podanego algorytmu w języku C#:
public double[] GaussElimination(double[,] A, double[] b, int n)
{
double[] x = new double[n];
double[,] tmpA = new double[n, n + 1];
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
tmpA[i, j] = A[i, j];
}
tmpA[i, n] = b[i];
}
double tmp = 0;
for (int k = 0; k < n - 1; k++)
{
for (int i = k + 1; i < n; i++)
{
tmp = tmpA[i, k] / tmpA[k, k];
for (int j = k; j < n + 1; j++)
{
tmpA[i, j] -= tmp * tmpA[k, j];
}
}
}
for (int k = n - 1; k >= 0; k--)
{
tmp = 0;
for (int j = k + 1; j < n; j++)
{
tmp += tmpA[k, j] * x[j];
}
x[k] = (tmpA[k, n] - tmp) / tmpA[k, k];
}
return x;
}
2. Eliminacja Gaussa - Jordana
Tak samo jak w poprzednim przypadku - rozpoczniemy od przykładu rozwiązanego tą metodą:
Zapisujemy układ jako macierz:
Krok 1: Z równań drugiego i trzeciego usuwamy x1, w tym celu równanie pierwsze mnożymy najpierw przez -2 a następnie przez 3:
Krok 2: Z równań pierwszego i trzeciego usuwamy x2, w tym celu równanie drugie mnożymy najpierw przez 2/5, a następnie przez 4/5:
Krok 3: Równania drugie i trzecie upraszczamy poprzez podzielenie przez 5 oraz 3:
Krok 3: Z równań pierwszego i drugiego eliminujemy x3, w tym celu równanie trzecie mnożymy najpierw przez -1 a następnie przez 1:
Rozwiązaniem układu są więc wartości:
Przejdźmy teraz do części programistycznej powyższej metody. Jeżeli ktoś wcześniej zapoznał się z metodą Gaussa, to nie wiele mamy do zmian. Więc jeżeli jeszcze nie zapoznałeś się z metodą eliminacji Gaussa – możesz to zrobić teraz. Jakie zmiany musimy nanieść w metodzie Gaussa aby uzyskać metodę Gaussa-Jordana? Właściwie tylko kosmetyczne:
- Każdy wiersz normalizujemy poprzez dzielenie go przez element główny
- Zmienne są eliminowane ze wszystkich równań a nie tylko następnych
- Po n krokach otrzymujemy macierz jednostkową – nie ma więc postępowania odwrotnego
Kod C#:
public double[] GaussJordanElimination(double[,] A, double[] b, int n)
{
double[] x = new double[n];
double[,] tmpA = new double[n, n + 1];
double tmp = 0;
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
tmpA[i, j] = A[i, j];
}
tmpA[i, n] = b[i];
}
for (int k = 0; k < n; k++)
{
tmp = tmpA[k, k];
for (int i = 0; i < n + 1; i++)
{
tmpA[k, i] = tmpA[k, i] / tmp;
}
for (int i = 0; i < n; i++)
{
if (i != k)
{
tmp = tmpA[i, k] / tmpA[k, k];
for (int j = k; j < n + 1; j++)
{
tmpA[i, j] -= tmp * tmpA[k, j];
}
}
}
}
for (int i = 0; i < n; i++)
{
x[i] = tmpA[i, n];
}
return x;
}
3. LU
W metodzie tej macierz A przedstawiamy jako iloczyn macierzy LU. Macierz L to dolna macierz trójkątna a U górna macierz trójkątna.
Do rozkładu na powyższe macierze L i U można użyć metody Doolittle’a. W metodzie równość A=LU traktuje się jako układ n2 równań z n2 niewiadomymi. Elementy macierzy L i U są wyznaczane naprzemiennie tzn. kolumnę macierzy L a później wiersz macierzy U itd. Wzory ogólne na rozkład:
Przykład:
Zapisujemy w postaci macierzy:
W pierwszym kroku obliczamy pierwszy wiersz macierzy U:
Teraz obliczymy kolumnę macierzy L:
Teraz znów wiersz macierzy U:
Kolejno obliczymy teraz ostatnią kolumnę macierzy L a następnie ostatni wiersz macierzy U:
Ostateczna postać macierzy L i U po rozkładzie:
Rozwiązaniem układu odbywa się kolejno przez wykonanie działań:
Tak więc pierwsze to po prostu rozwiązanie macierzy trójkątnej górnej, a kolejne równanie to rozwiązanie macierzy trójkątnej górnej (jak w metodzie Gaussa).
Przykładowy kod w C# tworzący macierz L, U oraz rozwiązujący układ równań:
public void LUDecompozition(double[,] A, double[,] L, double[,] U, int n)
{
for (int i = 0; i < L.GetLength(0); i++)
{
L[i, i] = 1;
}
double tmp;
for (int i = 0; i < n; i++)
{
for (int j = i; j < n; j++)
{
tmp = 0;
for (int k = 0; k < i; k++)
{
tmp += L[i, k] * U[k, j];
}
U[i, j] = A[i, j] - tmp;
}
for (int j = i + 1; j < n; j++)
{
tmp = 0;
for (int k = 0; k < i; k++)
{
tmp += L[j, k] * U[k, i];
}
L[j, i] = (A[j, i] - tmp) / U[i, i];
}
}
}
public double[] SolvWithLUMethod(double[,] L, double[,] U, double[] y, int n)
{
double[] z = new double[n];
double[] x = new double[n];
z = SolvLowerMatrix(L, y, n);
x = SolvUpperMatrix(U, z, n);
return x;
}
public double[] SolvLowerMatrix(double[,] L, double[] b, int n)
{
double[] x = new double[n];
x[0] = b[0] / L[0, 0];
double tmp;
for (int i = 1; i < n; i++)
{
tmp = 0;
for (int j = 0; j < i; j++)
{
tmp += L[i, j] * x[j];
}
x[i] = b[i] - tmp;
}
return x;
}
public double[] SolvUpperMatrix(double[,] U, double[] b, int n)
{
double[] x = new double[n];
double tmp;
for (int k = n - 1; k >= 0; k--)
{
tmp = 0;
for (int j = k + 1; j < n; j++)
{
tmp += U[k, j] * x[j];
}
x[k] = (b[k] - tmp) / U[k, k];
}
return x;
}
4. Metoda Jacobiego
Metoda ta jest zaliczona do metod iteracyjnych (Gauss, Gauss-Jordan, LU to metody skończone). Metody iteracyjne nadają się do rozwiązywania bardzo dużych układów równań – np. powyżej miliona równań.
W metodzie tej tworzy się ciąg przybliżeń według wzoru:
Przykładowa implementacja w C#:
public double[] JacobiMethod(double[,] A, double[] b, int n, double eps)
{
double tmp;
double[] x = new double[n];
double[] x1 = new double[n];
double sumx, sumx1;
int count = 0;
for (int i = 0; i < n; i++)
{
x1[i] = b[i] / A[i, i];
}
do
{
for (int i = 0; i < n; i++)
{
x[i] = x1[i];
}
for (int i = 0; i < n; i++)
{
tmp = 0;
for (int j = 0; j < n; j++)
{
if (j != i)
{
tmp += A[i, j] * x[j];
}
}
x1[i] = (1.0 / A[i, i]) * (b[i] - tmp);
}
sumx = 0;
sumx1 = 0;
for (int i = 0; i < n; i++)
{
sumx += x[i];
sumx1 += x1[i];
}
++count;
} while (Math.Abs(sumx - sumx1) > eps);
return x1;
}
5. Metoda Gaussa - Seidela
Metoda ta, jak i Jacobiego należy do metod iteracyjnych. W porównaniu do metody Jacobiego, metoda Gaussa – Seidela każda modyfikacja rozwiązania próbnego korzysta ze wszystkich dostępnych przybliżeń. Dzięki temu można uzyskać praktycznie dwukrotne zmniejszenie ilości powtórzeń w celu otrzymania zadanej dokładności w porównaniu z metodą Jacobiego.
Kolejne przybliżenia w tej metodzie oblicza się korzystając ze wzoru:
Przykładowy kod w C# :
public double[] GaussaSeidela(double[,] A, double[] b, int n, double eps)
{
double tmp1;
double tmp2;
double[] x = new double[n];
double[] x1 = new double[n];
double sumx, sumx1;
int count = 0;
for (int i = 0; i < n; i++)
{
x1[i] = b[i] / A[i, i];
}
do
{
for (int i = 0; i < n; i++)
{
x[i] = x1[i];
}
for (int i = 0; i < n; i++)
{
tmp1 = 0;
tmp2 = 0;
for (int j = 0; j < i; j++)
{
tmp1 += A[i, j] * x1[j];
}
for (int j = i + 1; j < n; j++)
{
tmp2 += A[i, j] * x[j];
}
x1[i] = (1.0 / A[i, i]) * (b[i] - tmp1 - tmp2);
}
sumx = 0;
sumx1 = 0;
for (int i = 0; i < n; i++)
{
sumx += x[i];
sumx1 += x1[i];
}
++count;
} while (Math.Abs(sumx - sumx1) > eps);
return x1;
}
Dla przykładowego układu równań:
Korzystając z metody Jacobiego i wymaganej dokładności eps = 0.0000001 wymagane są 24 iteracje. W przypadku metody Gaussa – Seidela tylko 13.
Prosty programik implementujący powyższe metody można pobrać stąd, a kod programu stąd. Należy wziąć pod uwagę aspekt, że nie uwzględniliśmy sytuacji gdy na głównej przekątnej macierzy znajduje się 0. Wtedy wykonanie kodu w takiej postaci jak jest spowoduje błąd podczas dzielenia. W takich sytuacjach można przykładowo przestawić wiersze i/lub kolumny. Ale to może kiedy indziej opiszę.
Prosty programik implementujący powyższe metody można pobrać stąd, a kod programu stąd. Należy wziąć pod uwagę aspekt, że nie uwzględniliśmy sytuacji gdy na głównej przekątnej macierzy znajduje się 0. Wtedy wykonanie kodu w takiej postaci jak jest spowoduje błąd podczas dzielenia. W takich sytuacjach można przykładowo przestawić wiersze i/lub kolumny. Ale to może kiedy indziej opiszę.
Brak komentarzy:
Prześlij komentarz