Η μέθοδος απαλοιφής του Gauss

Η μέθοδος απαλοιφής του Gauss είναι η πιο γνωστή τεχνική για την επίλυση γραμμικών συστημάτων. Αποτελεί αναπόσπαστο κομμάτι κάθε μαθήματος γραμμικής άλγεβρας σε οποιαδήποτε σχολή. Μια στοιχειώδης μορφή του αλγορίθμου εμφανίζεται και στο βιβλίο της Άλγεβρας Β΄ Λυκείου. Επομένως, οι περισσότεροι πρωτοετείς φοιτητές (είτε σε ΑΕΙ, είτε σε ΤΕΙ είτε στο ΕΑΠ) πιθανότατα κάπου θα την έχουν ακούσει. Όσον αφορά το ΕΑΠ διδάσκεται στην ενότητα ΠΛΗ12. Όσοι δεν τον γνωρίζετε μπορείτε να τον διαβάσετε στις σημειώσεις μας (Γραμμικά Συστήματα). Σε αυτό το άρθρο θα δώσουμε τον βασικό αλγόριθμο της μεθόδου καθώς επίσης και κώδικα σε python και Matlab για όσους ενδιαφέρονται να το τρέξουν. Ο κώδικας που δίνουμε μπορεί να σας βοηθήσει να ελέγξετε τα αποτελέσματά που βρίσκετε σε τυχόν ασκήσεις που λύνετε στο χαρτί αφού δίνει τη μετατροπή σε κλιμακωτό πίνακα βήμα-βήμα.

Α. Κλιμακωτοί πίνακες (row-Echelon forms)

Ένας πίνακας θα λέγεται κλιμακωτός αν

Όλες οι μη μηδενικές γραμμές του πίνακα βρίσκονται πάνω από τις μηδενικές γραμμές του πίνακα.
Το πρώτο μη μηδενικό στοιχείο κάθε γραμμής βρίσκεται δεξιότερα του αντίστοιχου πρώτου μη μηδενικού στοιχείου της προηγούμενης γραμμής
Το πρώτο μη μηδενικό στοιχείο κάθε γραμμής είναι ίσο με 1.

Β. Ο Αλγόριθμος

Ο αλγόριθμος απαλοιφής του Gauss (Gauss elimination) μπορεί να μετατρέψει κάθε πίνακα σε κλιμακωτό πίνακα. Αν λοιπόν πάρουμε τον επαυξημένο πίνακα του συστήματος και τον φέρουμε σε κλιμακωτή μορφή, τότε μπορούμε εύκολα να λύσουμε το σύστημα.

-----------------------------------------------------------------------------------------

Δίνεται ένας πίνακας Α με διαστάσεις ΜxΝ.

Θέτω i0=1

Για κάθε j=1,2,...N επανάλαβε

Βρες το πρώτο μη μηδενικό στοιχείο της στήλης j (*)

Αν τέτοιο στοιχείο δεν υπάρχει τότε

προχώρα στην επόμενη στήλη

Αλλιώς

Έστω ότι το πρώτο μη μηδενικό στοιχείο της στήλης j

είναι το Α[t,j]

Αντιμετάθεσε τις γραμμές i0 και t.

Διαίρεσε όλα τα στοιχεία της γραμμής i0 με το Α[i0,j]

Για i=i0+1,i0+2,...,M επανάλαβε

Από τα στοιχεία της γραμμής i (R_i) αφαίρεσε τα

στοιχεία της γραμμής i0 (R0) πολλαπλασιασμένα με

το A[i,j]. Δηλαδή:

R_i <-- R_i - A[i,j]*R0

τέλος_επανάληψης

i0 <-- i0 + 1

τέλος_αν

Αν i0>M τότε

τερματισμός.

τέλος_επανάληψης

------------------------------------------------------------------------------------------

Παρατηρείστε ότι σε κάθε επαναληπτικό βήμα, ο αλγόριθμος διαιρεί όλα τα στοιχεία μιας γραμμής με το πρώτο μη μηδενικό στοιχείο της στήλης που βρίσκει (βήμα *). Αν αυτό το στοιχείο είναι πολύ κοντά στο μηδέν, τότε ο αλγόριθμος μπορεί να οδηγήσει σε σοβαρά υπολογιστικά σφάλματα. Για να αντιμετωπιστούν τέτοια θέματα συνήθως εφαρμόζεται η μέθοδος της μερικής περιστροφής (partial pivoting) κατά την οποία αντικαθισούμε το βήμα (*) με:

Βρες το μεγαλύτερο (κατά απόλυτη τιμή) μη μηδενικό στοιχείο της στήλης j (**)

Γ. Ο Αλγόριθμος σε γλώσσα python

Παρακάτω δίνουμε έναν απλό κώδικα που υλοποιεί τον αλγόριθμο απαλοιφής του Gauss σε γλώσσα python. Στον συγκεκριμενο κώδικα ένας πίνακας υλοποιείται ως μια λίστα από λίστες. Π.χ. ο πίνακας

υλοποιείται ως

A=[[1,2,3],[4,5,6]]

Επιλέξαμε αυτή τη μορφή εισόδου, ώστε να μπορείτε να τρέξετε τον αλγόριθμο χωρίς την είσοδο άλλων βιβλιοθηκών (όπως της NumPy).

Στην αρχή του κώδικα μπορείτε να αλλάξετε τον πίνακα Α (δίνονται παραδείγματα) και να ελέγξετε τις μεταβλητές what_to_print (0,1,2 - επιλέγετε πόσο αναλυτικά θα δοθεί η λύση) και move_largest_to_pivot (1,0 - επιλέγετε αν θα εφαρμοστεί το partial pivoting ή όχι). Δεν χρειάζονται ειδικές γνώσεις python για να κάνετε τις αλλαγές, απλώς κατεβάστε τον διερμηνευτή της γλώσσας python και τρέξτε τον κώδικα που δίνουμε. Να σημειώσουμε ότι ο κώδικας που δίνουμε σε python εκτελεί τους μετασχηματισμούς μεταξύ των γραμμών, σε όλο το εύρος μιας γραμμής. Δηλαδή ακόμα και στις θέσεις που οι γραμμές έχουν μηδενικά στοιχεία. Θα μπορούσαμε να αλλάξουμε τον κώδικα ώστε να αποφεύγει τις άχρηστες πράξεις, αλλά τότε θα καταλήγαμε σε ένα πιο δυσανάγνωστο κώδικα.

#This function prints the elements of the list #as they were elements of a matrix def PrintMatrix( A ): for line in A: for element in line: print('{:2.2f}'.format(element), end='\t') else: print("\n",end='') return; def LinComb(li1, li2, a): li1[:] = [x-a*y for x,y in zip(li1,li2)] def Rescale(li,a): li[:] = [x/a for x in li] #Gauss Elimination Method #Change the variables what_to_print, move_largest_to_pivot to your liking #Give the matrix as a list of lists A = [[1, -1, 1, -2, -2], [2,-1,0,1,8],[1,1,-1,0,0],[3,2,2,-1,3]] #A = [[2,-1,1,1,6],[1,1,1,1,2],[-1,-1,1,-1,4]] #A = [[1,1,1,2],[1,-1,-1,0],[-2,1,1,0]] #A = [[1,1,1],[1,2,2],[1,-1,1]] #What do you want to print? #0 - print only the final matrix #1 - print the final stage of each step #2 - print everything what_to_print = 2 #move the entry with the largest absolute value to the pivot position #0 no #1 yes move_largest_to_pivot = 0 #This is the algorithm print("Original Matrix:") PrintMatrix(A) print("**********************************\n") B = A M = len(A) N = len(A[0]) i0=0; for j in range(0,N): if move_largest_to_pivot==1: #Find maximum element of column j maxim = abs(B[i0][j]) thesis_max = 0 for i,x in enumerate(B[i0:]): if abs(x[j])>maxim : maxim = abs(x[j]) thesis_max = i #swap row thesis with row j thesis_max += i0 temp = B[thesis_max] B[thesis_max]=B[i0] B[i0] = temp if what_to_print>=2: print("-----Move to Pivot Step {}-----".format(j+1)) PrintMatrix(B) else: #Find first non-zero element of column j for i,x in enumerate(B[i0:]): if x[j]!=0 : thesis = i break; else: thesis = 0 thesis += i0 #swap row thesis with row j temp = B[thesis] B[thesis]=B[i0] B[i0] = temp if what_to_print>=2: print("-----Move to Pivot Step {}-----".format(j+1)) PrintMatrix(B) if B[i0][j]!=0: Rescale(B[i0],B[i0][j]) if what_to_print>=2: print("--Normalize row {}--".format(i0+1)) PrintMatrix(B) pivot = B[i0] for i,line in enumerate(B[i0+1:]): if line[j]!=0:

LinComb(line,pivot,line[j])

if what_to_print>=2:

print("--Step {}.{}--".format(j+1,i+1))

PrintMatrix(B) if what_to_print==1: print("--Step {}--".format(j+1)) PrintMatrix(B) i0+=1 if i0==M: break print("**********************************\n") print("Transformed Matrix:") PrintMatrix(B)

Αν τρέξετε τον κώδικα χωρίς καμμιά αλλαγή θα δώσει το παρακάτω output

Original Matrix: 1.00 -1.00 1.00 -2.00 -2.00 2.00 -1.00 0.00 1.00 8.00 1.00 1.00 -1.00 0.00 0.00 3.00 2.00 2.00 -1.00 3.00 ********************************** -----Move to Pivot Step 1----- 1.00 -1.00 1.00 -2.00 -2.00 2.00 -1.00 0.00 1.00 8.00 1.00 1.00 -1.00 0.00 0.00 3.00 2.00 2.00 -1.00 3.00 --Normalize row 1-- 1.00 -1.00 1.00 -2.00 -2.00 2.00 -1.00 0.00 1.00 8.00 1.00 1.00 -1.00 0.00 0.00 3.00 2.00 2.00 -1.00 3.00 --Step 1.1-- 1.00 -1.00 1.00 -2.00 -2.00 0.00 1.00 -2.00 5.00 12.00 1.00 1.00 -1.00 0.00 0.00 3.00 2.00 2.00 -1.00 3.00 --Step 1.2-- 1.00 -1.00 1.00 -2.00 -2.00 0.00 1.00 -2.00 5.00 12.00 0.00 2.00 -2.00 2.00 2.00 3.00 2.00 2.00 -1.00 3.00 --Step 1.3-- 1.00 -1.00 1.00 -2.00 -2.00 0.00 1.00 -2.00 5.00 12.00 0.00 2.00 -2.00 2.00 2.00 0.00 5.00 -1.00 5.00 9.00 -----Move to Pivot Step 2----- 1.00 -1.00 1.00 -2.00 -2.00 0.00 1.00 -2.00 5.00 12.00 0.00 2.00 -2.00 2.00 2.00 0.00 5.00 -1.00 5.00 9.00 --Normalize row 2-- 1.00 -1.00 1.00 -2.00 -2.00 0.00 1.00 -2.00 5.00 12.00 0.00 2.00 -2.00 2.00 2.00 0.00 5.00 -1.00 5.00 9.00 --Step 2.1-- 1.00 -1.00 1.00 -2.00 -2.00 0.00 1.00 -2.00 5.00 12.00 0.00 0.00 2.00 -8.00 -22.00 0.00 5.00 -1.00 5.00 9.00 --Step 2.2-- 1.00 -1.00 1.00 -2.00 -2.00 0.00 1.00 -2.00 5.00 12.00 0.00 0.00 2.00 -8.00 -22.00 0.00 0.00 9.00 -20.00 -51.00 -----Move to Pivot Step 3----- 1.00 -1.00 1.00 -2.00 -2.00 0.00 1.00 -2.00 5.00 12.00 0.00 0.00 2.00 -8.00 -22.00 0.00 0.00 9.00 -20.00 -51.00 --Normalize row 3-- 1.00 -1.00 1.00 -2.00 -2.00 0.00 1.00 -2.00 5.00 12.00 0.00 0.00 1.00 -4.00 -11.00 0.00 0.00 9.00 -20.00 -51.00 --Step 3.1-- 1.00 -1.00 1.00 -2.00 -2.00 0.00 1.00 -2.00 5.00 12.00 0.00 0.00 1.00 -4.00 -11.00 0.00 0.00 0.00 16.00 48.00 -----Move to Pivot Step 4----- 1.00 -1.00 1.00 -2.00 -2.00 0.00 1.00 -2.00 5.00 12.00 0.00 0.00 1.00 -4.00 -11.00 0.00 0.00 0.00 16.00 48.00 --Normalize row 4-- 1.00 -1.00 1.00 -2.00 -2.00 0.00 1.00 -2.00 5.00 12.00 0.00 0.00 1.00 -4.00 -11.00 0.00 0.00 0.00 1.00 3.00 ********************************** Transformed Matrix: 1.00 -1.00 1.00 -2.00 -2.00 0.00 1.00 -2.00 5.00 12.00 0.00 0.00 1.00 -4.00 -11.00 0.00 0.00 0.00 1.00 3.00 Με λίγα λόγια, καλούμε τον αλγόριθμο να μετατρέψει τον πίνακα

σε κλιμακωτό πίνακα. Ο κώδικας εκτυπώνει όλα τα βήματα και στο τέλος μας δίνει την κλιμακωτή μορφή του πίνακα:

Δ. Ο Αλγόριθμος σε Matlab

Παρακάτω δίνουμε τον ίδιο αλγόριθμος σε υλοποίηση Matlab. Οι παράμετροι είναι παρόμοιοι. Ο κώδικας μπορεί να εκτυπώσει όλα τα βήματα τις μεθόδου. Ο κώδικας σε Matlab είναι αρκετά πιο αποτλεσματικός γιατί δεν κάνει άχρηστες πράξεις (όπως στην περίπτωση της python).

%Gauss Elimination Method %Change the variables what_to_print, move_largest_to_pivot to your liking %Give the matrix as a list of lists A = [1, -1, 1, -2, -2; 2,-1,0,1,8; 1,1,-1,0,0 ; 3,2,2,-1,3]; %A = [2,-1,1,1,6; 1,1,1,1,2; -1,-1,1,-1,4]; %A = [1,1,1,2; 1,-1,-1,0; -2,1,1,0]; %What do you want to print? %0 - print only the final matrix %1 - print the final stage of each step %2 - print everything what_to_print = 1; %move the entry with the largest absolute value to the pivot position %0 no %1 yes move_largest_to_pivot = 0; %This is the algorithm A disp('**********************************'); B = A; [M,N] = size(A); i0=1; for j=1:N if move_largest_to_pivot==1 %Find maximum element of column j max = abs(B(i0,j)); thesis_max = i0; for i=i0+1:M if abs(B(i,j))>max max = abs(B(i,j)); thesis_max = i; end end %swap row thesis with row j temp = B(thesis_max,:); B(thesis_max,:)=B(i0,:); B(i0,:) = temp; if what_to_print>=2 disp(['-----Move to Pivot Step ', num2str(j), '-----']); B end else %Find maximum element of column j thesis = i0; for i=i0:M if abs(B(i,j))~=0 thesis = i; break; end end %swap row thesis with row j temp = B(thesis,:); B(thesis,:)=B(i0,:); B(i0,:) = temp; if what_to_print>=2 disp(['-----Move to Pivot Step ', num2str(j), '-----']); B end end if B(i0,j)~=0 B(i0,:) = B(i0,:)/B(i0,j); if what_to_print>=2 disp(['--Normalize row', num2str(i0), '--']); B end for i=i0+1:M

if B(i,j)~=0 B(i,j+1:N) = B(i,j+1:N) - B(i,j)*B(i0,j+1:N); if what_to_print>=2 print(['--Step ', num2str(j), '.', num2str(i), '--']); B end

end end if what_to_print==1 disp(['--Step', num2str(j), '--']); B end i0=i0+1; end if i0==M+1 break; end end disp('**********************************'); disp('Transformed Matrix:'); B