Rozwiązywanie układów równań
Będzie nam potrzebny cały kod z implementacją macierzy z poprzednich zajęć.
Jeżeli nie było Cię na poprzednich zajęciach, nie udało Ci się dokończyć wszystkich zadań
lub masz wątpliwości, czy twoje rozwiązania są poprawne, możesz pobrać
stąd naszą implementację klasy Matrix .
Następnie obok pliku matrices.py gauss.py
from matrices import Matrix , read_matrix_file , write_matrix_file
import sys
def reduce_to_echelon ( m ):
... # TODO: funkcja schodkująca macierze
def compute_solution ( echelon_m ):
... # TODO: funkcja wyliczająca rozwiązanie z wyschodkowanej macierzy
if __name__ == "__main__" and len ( sys . argv ) >= 2 :
m = read_matrix_file ( sys . argv [ 1 ])
reduce_to_echelon ( m )
sol = compute_solution ( m )
print ( f "Rozwiązanie: { sol } " )
if len ( sys . argv ) >= 3 :
write_matrix_file ( m , sys . argv [ 2 ])
Celem na dziś będzie napisanie programu rozwiązującego układ równań:
\[\begin{split}
\begin{cases}
m_{1\,1} x_1 + m_{1\,2} x_2 + \dots + m_{1\,n} x_n = m_{1\,n+1} \\
m_{2\,1} x_1 + m_{2\,2} x_2 + \dots + m_{2\,n} x_n = m_{2\,n+1} \\
\dots \\
m_{n\,1} x_1 + m_{n\,2} x_2 + \dots + m_{n\,n} x_n = m_{n\,n+1} \\
\end{cases}
\end{split}\]
Czyli układu opisanego macierzą:
\[\begin{split}
M = \begin{cases}
m_{1\,1} & m_{1\,2} & \cdots & m_{1\,n} & m_{1\,n+1} \\
m_{2\,1} & m_{2\,2} & \cdots & m_{2\,n} & m_{2\,n+1} \\
\vdots & \vdots & \ddots & \vdots & \vdots \\
m_{n\,1} & m_{n\,2} & \cdots & m_{n\,n} & m_{n\,n+1} \\
\end{cases}
\end{split}\]
Napiszemy do tego dwie funkcje:
Wtedy przy pomocy powyższego programu będziemy mogli rozwiązywać układy równań:
python gauss.py plik_z_macierza.txt
Ważne
Jeżeli w którymkolwiek momencie coś nie będzie Ci działać, możesz użyć funkcji
write_matrix_file(m) m print
Schodkowanie macierzy
Zajmiemy się teraz implementacją metody Gaussa.
Zadanie 1: redukcja kolumny
Zadanie 1
Napisz funkcję reduce_column_down(m, k) m Matrix.add_row k k+1
Wskazówka
Dla każdego wiersza o numerze \(i\colon n \ge i > k\) dodawaj \(k\) -ty wiersz przemnożony przez
\(-\frac{m_{i\,k}}{m_{k\,k}}\) , gdzie \(n\) to liczba wierszy w \(M\) .
Przetestuj swoją funkcję na jakiejś przykładowej macierzy.
Zadanie 2: bardzo naiwny Gauss
Zadanie 2
Korzystając z funkcji reduce_column_down reduce_to_echelon(m)
Zadanie 2.5
Co może pójść nie tak? Wymyśl macierz, dla której to podejście nie zadziała i zobacz, jaki wynik
otrzymasz.
Wskazówka
Kiedy ten program wykona dzielenie przez zero?
Zadanie 3: pół-pełny Gauss
Żeby uniknąć dzielenia przez zero (o ile to w ogóle możliwe — może się okazać, że rozwiązanie
nie jest jednoznaczne i w którymś momencie dostaniemy same zera).
Przykładowo mając w drugim kroku (kolumna o indeksie 1
\[\begin{split}\begin{pmatrix}
1 & 20 & 4 & 8 & 7 & 10 \\
0 & 1 & 5 & 2 & 7 & -3 \\
0 & -6 & 3 & 2 & 1 & 0 \\
0 & 5 & 2 & 1 & 1 & -7
\end{pmatrix}\end{split}\]
Zanim wywoła reduce_column_down(m, 1) m.swap_rows(1, 2)
Zadanie 3
Zmodyfikuj poprzednie rozwiązanie tak, żeby zamiast brać po prostu kolejny wiersz, najpierw
zamieniało wiersze miejscami w taki sposób, żeby do wyzerowania użyć wiersza mającego
w danej kolumnie największą co do wartości bezwzględnej wartość.
Wskazówka
Dla wywołania reduce_column_down(m, k)
Najpierw w pętli znajdź takie \(i\) , że \(n \ge i \ge k\) oraz:
$\( \lvert m_{i\,k}\rvert = max\{\lvert m_{k\,k}\rvert, \lvert m_{k+1\,k}\rvert, \dots, \lvert m_{n\,k}\rvert\} \) $
Następnie zamień wiersze \(i\) -ty z \(k\) -tym.
Wykonaj resztę operacji tak jak w poprzedniej wersji.
Obliczanie rozwiązania z wyschodkowanej macierzy
Zadanie 4: obliczanie pojedynczej zmiennej
Zaczniemy od funkcji, która mając obliczone wartości zmiennych \(x_{k+1}, \dots, x_n\) oraz
pojedyncze równanie w już wyschodkowanej macierzy obliczy z niego kolejną zmienną.
Ta sytuacja jest przedstawiona na rysunku:
Podświetlony wiersz odpowiada równaniu:
\[ 0x_1 + 0x_2 + \dots + m_{k\,k} x_k + m_{k\,k+1} x_{k+1} + m_{k\,k+2} x_{k+2} + \dots + m_{k\,n} x_n = m_{k\,n+1} \]
Na tym etapie rozwiązywania układu równań wartości zmiennych \(x_{k+1}, \dots, x_n\) są już wyliczone.
Zatem kolejną zmienną \(x_k\) wyliczamy jako:
\[ x_k = \frac{1}{m_{k\,k}}\left(m_{k\,n+1} - m_{k\,k+1} x_{k+1} - m_{k\,k+2} x_{k+2} - \dots - m_{k\,n} x_n\right) \]
W programie współczynniki będziemy brać z macierzy — \(m_{k\,j} = \texttt{m.get(k - 1, j - 1)}\) .
Wartości zmiennych będziemy pobierać i zapisywać do pojedynczej listy \(x_j = \texttt{sol[j - 1]}\) .
Później możemy kilka razy wywołać tę funkcję, aby wyliczyć całe rozwiązanie układu z wyschodkowanej
macierzy:
m = Matrix ( 3 , 4 )
m . set ( 0 , 0 , 1 ); m . set ( 0 , 1 , 5 ); m . set ( 0 , 2 , 4 ); m . set ( 0 , 3 , 9 )
m . set ( 1 , 1 , 1 ); m . set ( 1 , 2 , - 3 ); m . set ( 1 , 3 , 2 )
m . set ( 2 , 2 , 1 ); m . set ( 2 , 3 , - 1 )
sol = [ None ] * 3 # Na początek nie znamy żadnej z trzech zmiennych
sol [ 2 ] = extend_solution ( m , 2 , sol )
sol [ 1 ] = extend_solution ( m , 1 , sol )
sol [ 0 ] = extend_solution ( m , 0 , sol )
if sol == [ 18 , - 1 , - 1 ]:
print ( "extend_solution: OK" )
else :
print ( "extend_solution: ERROR" )
Zadanie 4
Napisz funkcję extend_solution(m, k, x)
Zadanie 5: obliczanie wszystkich zmiennych
Zadanie 5
Teraz korzystając z funkcji extend_solution compute_solution
Wskazówka
Niech \(n = \texttt{m.num_rows()}\) .
Zacznij od zainicjowania listy na wynik jako sol = [None] * n
Zaczynając od \(n\) -tego wiersza, kończąc na pierwszym przypisuj kolejne wartości rozwiązania
korzystając z funkcji extend_solution
Poniżej masz kilka macierzy schodkowych, których możesz użyć do przetestowania:
Skopiuj i zapisz wartości macierzy do pliku tekstowego.
Nie kopiuj pierwszej linijki z rozwiązaniem .
Wykonaj polecenie w terminalu: python gauss.py <nazwa_pliku>
Porównaj otrzymane rozwiązanie.
# Rozwiązanie: [4, 2]
1 3 10
0 1 2
# Rozwiązanie: [0, 3, 2]
1 2 1 8
0 1 4 11
0 0 1 2
# Rozwiązanie: [2, -1, 3, 1, 3]
1 1 2 0 1 10
0 1 3 1 2 15
0 0 1 2 1 8
0 0 0 1 3 10
0 0 0 0 1 3
# Rozwiązanie: [-11868, 4832, -1446, 249, -55, 18, -1, -1]
1 2 - 1 3 0 1 4 - 2 5
0 1 3 - 2 1 4 0 1 12
0 0 1 5 - 3 2 1 3 - 4
0 0 0 1 4 - 1 2 2 7
0 0 0 0 1 3 0 - 1 0
0 0 0 0 0 1 5 4 9
0 0 0 0 0 0 1 - 3 2
0 0 0 0 0 0 0 1 - 1
Zadania domowe
Zadanie g.err [1 punkt]
Podaj przykład macierzy układu, który: a) nie ma rozwiązań; b) ma nieskończenie wiele rozwiązań.
Uruchom na tych macierzach zaimplementowaną na ćwiczeniach funkcję do schodkowania.
Jaki jest efekt?
Zadanie g.solve [1 punkt]
Napisz funkcję solve(m) m
Zadanie g.better_gauss [2 punkty]
Popraw implementację funkcji reduce_to_echelon raise ValueError("nieobsługiwany układ") \(10^{-6}\) traktuj jako zera.
Uwaga: możesz dostać punkty tylko za jedno z zadań g.better_gauss g.better_gauss2
Zadanie g.coord [2 punkty]
Napisz funkcję compute_coordinates(vector, base)
Argument vector
Argument base
Wynik to lista współrzędnych, gdzie \(i\) -ty element odpowiada \(i\) -temu wektorowi z listy base
Załóż, że podany argument jest poprawną bazą przestrzeni liniowej.
Na przykład compute_coordinates([3, 6], [[1, 0], [0, 2]]) [3, 3]
Zadanie g.better_gauss2 [3 punkty]
Rozwiąż zadanie g.better_gauss \(10^{-6}\)
wyznacz tę granicę na podstawie oryginalnej macierzy:
Przed rozpoczęciem schodkowania znajdź największy co do wartości bezwzględnej element macierzy.
Traktuj jako zera te wartości, które są mniejsze od \(10^{-6}m\) , gdzie \(m\) to wartość bezwzględna
znalezionego wcześniej największego elementu.
Uwaga: możesz dostać punkty tylko za jedno z zadań g.better_gauss g.better_gauss2
Zadanie g.raycast [3 punkty]
Uzupełnij implementację poniższej funkcji, która oblicza przecięcie płaszczyzny z prostą:
class Point3D :
# Punkt w przestrzeni trójwymiarowej.
def __init__ ( self , x , y , z ):
self . x = x
self . y = y
self . z = z
class Plane3D :
# Płaszczyzna w przestrzeni trójwymiarowej dana za pomocą trzech punktów na niej leżących.
def __init__ ( self , p1 , p2 , p3 ):
self . p1 = p1
self . p2 = p2
self . p3 = p3
class Line3D :
# Prosta w przestrzeni trójwymiarowej przechodząca przez podane punkty.
def raycast ( plane , line ):
... # TODO...
return intersection # zwraca Point3D, w którym prosta przecina linię.
Załóż, że prosta nie jest równoległa do tej płasczyzny.
Argumenty — płaszczyzna i prosta — są podawane to funkcji raycast
Uwaga: częścią rozwiązania jest również matematyczne uzasadnienie działania tej funkcji.
Możesz je dołączyć jako zdjęcie.
