Γραμμικά συστήματα

4. Γραμμικά συστήματα#

Ένα σύστημα \(m\) γραμμικών εξισώσεων με \(n\) αγνώστους εκφράζεται στην γενική του μορφή ως:

\[\begin{split} \begin{aligned} &a_{1,1} x_1 & + & a_{1,2} x_2 && {\ldots}& +&a_{1,n} x_n &=& b_1,\\ &a_{2,1} x_1 & + & a_{2,2} x_2 &&{\ldots}& +&a_{2,n} x_n &=& b_2, \\ &{\ldots} & &{\ldots} &&{\ldots}&&{\ldots}&&{\ldots} \\ &a_{m,1} x_1 & + & a_{m,2}x_2 &&{\ldots}& +& a_{m,n} x_n &=& b_{m}. \end{aligned} \end{split}\]

ή ισοδύναμα:

\[Α x =b \]

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

\(m=n\), αλλά κάποιες εξισώσεις είναι γραμμικά εξαρτημένες,
\(m<n\), δηλαδή οι εξισώσεις δεν επαρκούν,
\(m>n\), δηλαδή έχουμε πλεονάζουσες εξισώσεις.

Παρακάτω θα ασχοληθούμε με προβλήματα όπου ο πίνακας \(Α\) είναι τετραγωνικός \((m=n)\). Η κατάστρωση του προβλήματος συντομεύεται με την χρήση πινάκων ως:

\[\begin{split} \begin{bmatrix} a_{1,1} & a_{1,2} & {\ldots}& a_{1,n} \\ a_{2,1} & a_{2,2} & {\ldots}& a_{2,n} \\ {\ldots} &{\ldots} & {\ldots}&{\ldots} \\ a_{n,1} & a_{n,2} & {\ldots}& a_{n,n} \\ \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ {\ldots}\\ x_n\\ \end{bmatrix} = \begin{bmatrix} b_1\\ b_2\\ {\ldots}\\ b_n\\ \end{bmatrix} \end{split}\]

4.1. Απαλοιφή Gauss#

Η απαλοιφή Gauss (Gauss elimination) στοχεύει στην μετατροπή του πίνακα \(Α\) σε άνω τριγωνικό, δηλαδή όλες οι τιμές κάτω από την διαγώνιο να έχουν μηδενική τιμή.

Για την πρώτη στήλη \((j=1)\) αυτό επιτυγχάνεται αφαιρώντας την πρώτη εξίσωση από όλες τις υπόλοιπες, αφότου πολλαπλασιασθεί με τον κατάλληλο συντελεστή \(\frac{a_{j,1}}{a_{1,1}}\). Αυτό επαναλαμβάνεται με όλες τις στήλες \((j=1,n)\).

Tip

Η γραμμή ή εξίσωση, που χρησιμοποιούμε για να μηδενίσουμε τις τιμές κάτω από την διαγώνιο, ονομάζεται οδηγός (pivot).

Στο τέλος το παραπάνω σύστημα μετασχηματίζεται σε:

\[\begin{split} \begin{aligned} &a_{1,1} x_1 & + & a_{1,2} x_2 && {\ldots}& +&a_{1,n} x_n &=& b_1,\\ & & + & {a_{2,2}}' x_2 &&{\ldots}& +& {a_{2,n}}' x_n &=& {b_2}' \\ &&&&&{\ldots} &&{\ldots}&& {\ldots}\\ & & & & & & & {a_{n,n}}' x_n &=& {b_{n}}' \end{aligned} \end{split}\]

με αλλαγές σε όλους τους συντελεστές εκτός από την πρώτη γραμμή.

Exercise 4.1

Επιλύστε το παρακάτω γραμμικό σύστημα:

\( \begin{aligned} & 4x_1 & +&2x_2 & +& x_3 &=& 6 \\ & x_1 & -& x_2 & +&2x_3 &=& 5 \\ & 2x_1 & +&2x_2 & -& x_3 &=& 0 \\ \end{aligned} \)

Λύση 1

Το αρχικό σύστημα γράφεται με την μορφή του επαυξημένου πίνακα \(\left[Α|b\right]\)

\[\begin{split} \begin{bmatrix} 4 & 2 & 1 & | & 6\\ 1 & -1 & 2 & | & 5\\ 2 & 2 & -1 & | & 0\\ \end{bmatrix} \end{split}\]

Τα στοιχεία της πρώτης στήλη κάτω από την διαγώνιο μηδενίζονται:

\[ \text{Γραμμή}2 = \text{Γραμμή}2-\frac{1}{4}\text{Γραμμή}1 \]

\[ \text{Γραμμή}3 = \text{Γραμμή}3-\frac{2}{4}\text{Γραμμή}1 \]

\[\begin{split} \begin{bmatrix} 4 & 2 & 1 & | & 6\\ 0 & -1.5 & 1.75 & | & 3.5\\ 0 & 1 & -1.5 & | & 3\\ \end{bmatrix} \end{split}\]

Τα στοιχεία της δεύτερης στήλης κάτω από την διαγώνιο μηδενίζονται:

\[ \text{Γραμμή3} = \text{Γραμμή3}-\left(\frac{1}{-1.5}\right)\text{Γραμμή2} \]

\[\begin{split} \begin{bmatrix} 4 & 2 & 1 & | & 6\\ 0 & -1.5 & 1.75 & | & 3.5\\ 0 & 0 & -0.33 & | & -0.66\\ \end{bmatrix} \end{split}\]

Ξεκινώντας από τελευταία σχέση, βρίσκουμε τον άγνωστο και τον αντικαθιστούμε στην προηγούμενες. Έτσι προκύπτει διαδοχικά:

\[x_3=\frac{-0.66}{-0.33}=2\]

\[ x_2=\frac{3.5-2*1.75}{-1.5}=0\]

\[x_1=\frac{6-2*0-1*2}{4}=1\]

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

Η μέθοδος που προκύπτει με την εφαρμογή αντιμετάθεσης γραμμών μόνο ονομάζεται απαλοιφή Gauss με μερική καθοδήγηση. Η επίλυση

4.2. Καταμέτρηση πράξεων#

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

import numpy as np

def solveGauss(A:np.array,b:np.array)->np.array:
    #   Original arrays are destroyed, no checks for array sizes and zero values

    # Forward elimination
    n=len(b)
    x=np.zeros(n,dtype=np.float64)

    for k in range(0,n-1):                  #   k=column to eliminate (last column is excluded)
        for i in range(k+1,n):              #   i=line to subtract from 
            factor=A[i,k]/A[k,k]                #   1 multiplication
            for j in range(k+1,n):          #   j=column in line to subtract from
                A[i,j]=A[i,j]-factor*A[k,j]     #   1 multiplication + 1 addition
            b[i]=b[i]-factor*b[k]               #   1 multiplication + 1 addition

    # Reverse substitution
    x[n-1]=b[n-1]/A[n-1,n-1]
    for i in range(n-2,-1,-1):           #   i=line
        sum=b[i]
        for j in range(i+1,n):        #   j=column
            sum=sum-A[i,j]*x[j]
        x[i]=sum/A[i,i]

    return x

A = np.array([[4, 2, 1], [1, -1, 2], [2, 2, -1]],dtype=np.float64)
b = np.array([6, 5, 0],dtype=np.float64)

print(solveGauss(A,b))

[ 1. -0.  2.]

Για την καταμέτρηση είναι χρήσιμα τα παρακάτω αθροίσματα:

\[\begin{split} \begin{align} \sum_{i=1}^{n} 1 &=& 1 + 1 + 1 + \dots + 1 &=& n\\ \sum_{i=1}^{n} i &=& 1 + 2 + 3 + \dots + n &=& \frac{n(n+1)}{2}&=&\frac{n^2}{2}+\mathcal{O}(n)\\ \sum_{i=1}^{n} i^2 &=& 1 + 2^2 + 3^2 + \dots + n^2 &=& \frac{n(n+1)(2n+1)}{6}&=&\underbrace{\frac{n^3}{3}}_\text{Αποτέλεσμα 3ης τάξης}+\underbrace{\mathcal{O}(n^2)}_\text{Σφάλμα 2ης τάξης} \end{align} \end{split}\]

Tip

Ο συμβολισμός \(\mathcal{O}(\dots)\) δείχνει τάξη μεγέθους. Ένα μέγεθος με \(\mathcal{O}(n^m)\) έχει τάξη μεγέθους ανάλογη με το \(n\) στην δύναμη \(m\). Οι όροι με δύναμη μικρότερης της \(m\) θεωρούνται αμελητέοι.

Το σύνολο των προσθέσεων και αφαιρέσεων που απαιτεί η φάση της απαλοιφής της μεθόδου Gauss δίνεται από τα αθροίσματα που αντιστοιχούν στους 3 βρόχους επανάληψης:

\[\begin{split} \begin{align} \sum_{k=1}^{n} \sum_{i=k+1}^{n} \sum_{j=k+1}^{n} 1 &=& \sum_{k=1}^{n} \sum_{i=k+1}^{n} \left(\sum_{j=1}^{n} 1-\sum_{j=1}^{k} 1\right)=\sum_{k=1}^{n} \sum_{i=k+1}^{n} (n-k)\\ &=&\sum_{k=1}^{n} \left(\sum_{j=1}^{n} (n-k)-\sum_{j=1}^{k} (n-k)\right)= \sum_{k=1}^{n} \left(n (n-k)-k (n-k)\right)= \sum_{k=1}^{n} \left(n-k\right)^2\\ &=&\sum_{k=1}^{n}\left(n^2 -2nk + k^2\right)=n^2n-2n[\frac{n^2}{2}+\mathcal{O}(n)]+\frac{n^3}{3}+\mathcal{O}(n^2)\\ &=&\frac{n^3}{3}+\mathcal{O}(n^2) \end{align} \end{split}\]

Η καταμέτρηση για τους πολλαπλασιασμούς και τις διαιρέσεις δίνει το ίδιο αποτέλεσμα επομένως το τελικό σύνολο των πράξεων είναι \(\frac{3}{2} n^3+\mathcal{O}(n^2)\). Με την εξαίρεση πολύ μικρών συστημάτων, ο όρος τάξης \(n^2\) είναι πολύ μικρότερος του όρου \(n^3\) και για αυτό συνήθως αγνοείται. Επίσης αυτό που μας ενδιαφέρει είναι κυρίως ο εκθέτης του \(n\) και λιγότερο ο συντελεστής (εδώ \(3/2\)). Επομένως το πλήθος πράξεων της απαλοιφής Gauss είναι της τάξεως \(\mathcal{Ο}(n^3)\). Η προς τα πίσω αντικατάσταση απαιτεί πλήθος πράξεων τάξης \(\mathcal{O}(n^2)\), επομένως αγνοείται όταν αναφέρεται H τάξη μεγέθους της συνολικής διαδικασίας.

4.3. Απαλοιφή Gauss-Jordan#

Στόχος της απαλοιφής Gauss-Jordan είναι ο μετασχηματισμός του πίνακα \(Α\) σε μοναδιαίο και του συστήματος στην μορφή:

\[\begin{split} \begin{bmatrix} 1 & 0 & {\ldots}&0 \\ 0 & 1 & {\ldots}& 0 \\ {\ldots} &{\ldots} & {\ldots}&{\ldots} \\ 0 & 0 & {\ldots}& 1 \\ \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ {\ldots}\\ x_n\\ \end{bmatrix} = \begin{bmatrix} {b_1}'\\ {b_2}'\\ {\ldots}\\ {b_n}'\\ \end{bmatrix} \end{split}\]

Η διαδικασία διαφοροποιείται σε δύο σημεία σε σχέση με την απαλοιφή Gauss:

Οι γραμμές κανονικοποιούνται με βάση το στοιχείο της διαγωνίου, δηλαδή η γραμμή πολλαπλασιάζεται έτσι ώστε να προκύψει μονάδα στην διαγώνιο.
Μετά την διαδικασία απαλοιφής των στοιχείων κάτω από την διαγώνιο, ακολουθεί η διαδικασία απαλοιφής των στοιχείων πάνω από την διαγώνιο.

Η μέθοδος Gauss-Jordan μπορεί να χρησιμοποιηθεί για την εύρεση του αντίστροφου ενός πίνακα. Η ίδια διαδικασία εφαρμόζεται στον επαυξημένο πίνακα \( \left[A|I\right]\) και καταλήγει στον πίνακα \( \left[Ι|Α^{-1}\right]\).

Για την επίλυση της Ασκησης 4.1, ακολουθούμε τα παρακάτω βήματα:

Λύση 2

Αρχικά εφαρμόζεται η κανονικοποίηση στην πρώτη γραμμή:

\[ \text{Γραμμή}1 =\frac{1}{4} \text{Γραμμή}1 \]

και χρησιμοποιείται για την απαλοιφή όλων των υποκείμενων στοιχείων στην πρώτη στήλη:

Στην συνέχεια εφαρμόζεται στην δεύτερη γραμμή:

\[ \text{Γραμμή}2 =\frac{1}{-1.5} \text{Γραμμή}2 \]

και χρησιμοποιείται για την απαλοιφή του υποκείμενου στοιχείου.

Τέλος κανονικοποιείται η τρίτη γραμμή:

\[ \text{Γραμμή}3 =\frac{1}{2.667} \text{Γραμμή}3 \]

\[\begin{split} \begin{bmatrix} 1 & 0.5 & 0.25 & | & 1.5\\ 0 & 1 & 1.167 & | & 2.333\\ 0 & 0 & 1 & | & 2\\ \end{bmatrix} \end{split}\]

Τα στοιχεία πάνω από την διαγώνιο στην τελευταία στήλη μηδενίζονται όπως στην απαλοιφή Gauss:

\[ \text{Γραμμή}2 =\text{Γραμμή}2 - \frac{1.167}{1.0} \text{Γραμμή}3 \]

\[ \text{Γραμμή}1 =\text{Γραμμή}1 - \frac{0.25}{1.0} \text{Γραμμή}3 \]

\[\begin{split} \begin{bmatrix} 1 & 0.5 & 0 & | & 1\\ 0 & 1 & 0 & | & 0\\ 0 & 0 & 1 & | & 2\\ \end{bmatrix} \end{split}\]

Μηδενίζονται τα στοιχεία στην δεύτερη στήλη:

\[ \text{Γραμμή}1 =\text{Γραμμή}1 - \frac{0.5}{1} \text{Γραμμή}2 \]

\[\begin{split} \begin{bmatrix} 1 & 0 & 0 & | & 1\\ 0 & 1 & 0 & | & 0\\ 0 & 0 & 1 & | & 2\\ \end{bmatrix} \end{split}\]

4.4. Παραγοντοποίηση LU#

Πολύ συχνά απαντώνται προβλήματα που απαιτούν επίλυση του ίδιου γραμμικού προβλήματος πολλές φορές με διαφοροποίηση μόνο στον πίνακα \(b\).

Σε αυτή την περίπτωση χρησιμοποιείται η παραγοντοποίηση LU, η οποία εξοικονομεί έναν σημαντικό αριθμό πράξεων. Το πρώτο βήμα είναι ο μετασχηματισμός του \(A\) ως

\[Α=LU\]

όπου

\(L\) ένας άνω τριγωνικός πίνακας
\(U\) ενας κάτω τριγωνικός πίνακας

Στην συνέχεια, αρκεί να αντικατασταθεί \(Ux=m\) στην εξίσωση \(LUx=b\), οπότε προκύπτει το σύστημα

\(Lm=b\)

Το διάνυσμα \(m\) μπορεί να λυθεί εύκολα με απαλοιφή καθώς ο πίνακας \(L\) είναι τριγωνικός.

Με γνωστό το \(m\), λύνεται και το σύστημα \(Ux=m\) πάλι με απλή απαλοιφή στον τριγωνικό πίνακα \(U\).

Το μόνο που λείπει πλέον είναι να βρεθούν οι πίνακες \(L\) και \(U\). Ένας τρόπος είναι η απαλοιφή Gauss. Στην λύση 1 της Ασκησης 4.1 που παρουσιάστηκε ήδη, ο πίνακας \(U\) προκύπτει κατά την απαλοιφή των στοιχείων κάτω από την διαγώνιο:

\[\begin{split} U= \begin{bmatrix} 4 & 2 & 1 \\ 0 & -1.5 & 1.75 \\ 0 & 0 & -0.33 \\ \end{bmatrix} \end{split}\]

και ο πίνακας \(L\) είναι οι συντελεστές που χρησιμοποιήθηκαν στην αφαίρεση των γραμμών:

\[\begin{split} L= \begin{bmatrix} 1 & 0 & 0 \\ 1/4 & 1 & 0 \\ 1/2 & -2/3 & 1 \\ \end{bmatrix} \end{split}\]

Η συνάρτηση της Python numpy.linalg.solve χρησιμοποιεί εσωτερικά μια υλοποίηση της παραγοντοποίησης LU που προέρχεται από την κλασσική βιβλιοθήκη LAPACK. Στην μονάδα κώδικα scipy διατίθεται η scipy.linalg.lu η οποία επιστρέφει τους τριγωνικούς πίνακες ξεχωριστά για περαιτέρω χρήση.

Exercise 4.2

Επιλύστε το γραμμικό σύστημα της Ασκησης 4.1 με την χρήση της μονάδας κώδικα numpy και εφαρμόστε την LU στον πίνακα \(A\) με την χρήση της μονάδας κώδικα scipy.

from scipy.linalg import lu
import numpy as np

A = np.array([[4, 2, 1], [1, -1, 2], [2, 2, -1]])
b = np.array([6, 5, 0])

x = np.linalg.solve(A, b)
print(x)

P, L, U = lu(A)     # Permutation matrix, Lower triangular, Upper triangular
print(f"\n{P}\n\n{L}\n\n{U}\n\n{L@U}")

[ 1. -0.  2.]

[[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]

[[ 1.          0.          0.        ]
 [ 0.25        1.          0.        ]
 [ 0.5        -0.66666667  1.        ]]

[[ 4.          2.          1.        ]
 [ 0.         -1.5         1.75      ]
 [ 0.          0.         -0.33333333]]

[[ 4.  2.  1.]
 [ 1. -1.  2.]
 [ 2.  2. -1.]]

Όταν χρειάζεται να επιλύσουμε πολλές φορές ένα σύστημα με σταθερό πίνακα \(A\) και μεταβλητό \(b\), ενδείκνυται η χρήση της παραγοντοποίησης LU. Πριν την πρώτη επίλυση υπολογίζονται οι τριγωνικοί πίνακες \(L\) και \(U\) με πλήθος πράξεων τάξης \(\mathcal{O}(n^3)\). Οι ίδιοι πίνακες χρησιμοποιούνται στην συνέχεια για όσες επιλύσεις απαιτούνται, καθώς παραμένουν ίδιοι. Η επίλυση με τους τριγωνικούς πίνακες είναι μια προς τα πίσω αντικατάσταση με πλήθος πράξεων τάξης \(\mathcal{O}(n^2)\). Επομένως όταν χρησιμοποιείται η παραγοντοποίηση LU o κύριος όγκος των πράξεων είναι τάξης \(\mathcal{O}(n^2)\) με σημαντική εξοικονόμηση πράξεων σε σχέση με την Gauss.

4.5. Επαναληπτικές Μέθοδοι#

Οι παραπάνω μέθοδοι ονομάζονται άμεσες καθώς απαιτούν ένα πεπερασμένο πλήθος αριθμητικών πράξεων για την επίλυση του γραμμικού συστήματος. Μία δεύτερη κατηγορία μεθόδων είναι οι επαναληπτικές ή εμμεσες. Η λογική αυτών των μεθόδων είναι η χρήση μια αρχικής εκτίμησης για την λύση και η σταδιακή βελτίωσή της μέχρι να περιοριστούν τα σφάλματα σε ικανοποιητικό βαθμό.

4.5.1. Μέθοδος Jacobi#

Έστω το σύστημα \(Α x =b\) και μια αρχική εκτίμηση \(x^{(0)}\), η οποία δεν ικανοποιεί απόλυτα την εξίσωση. Υπό προϋποθέσεις μπορούμε να πάρουμε μια καλύτερη εκτίμηση \(x^{(1)}\), αν επιλύσουμε κάθε εξίσωση του συστήματος θεωρώντας γνωστά όλα τα \(x_i\) εκτός από ένα.

\[x_i^{\left(t\right)}=\frac{1}{a_{i,i}}\left[b_i-\sum_{j=1}^{j=i-1}a_{i,j}x^{(t-1)}_j -\sum_{j=i+1}^{j=n} a_{i,j} x^{(t-1)}_j \right]\]

Επαναλαμβάνοντας πολλές φορές την διαδικασία, οι νέες τιμές \(x^{\left(t\right)}\) έρχονται πολύ κοντά στις προηγούμενες \(x^{(t-1)}\), σημαίνει ότι η λύση συγκλίνει και μπορούμε να σταματήσουμε αναλόγως το αποδεκτό επίπεδο ανοχής (tolerance). Η παραπάνω μέθοδος φέρει την ονομασία Jacobi και σημαντικό στοιχείο της είναι ότι το \(x_i^{\left(t\right)}\) υπολογίζεται μόνο με τις εκτιμήσεις του προηγούμενου βήματος \(t-1\).

4.5.2. Gauss-Seidel#

Επιτάχυνση της σύγκλισης μπορεί να επιτευχθεί με την μέθοδο Gauss-Seidel, κατά την οποία χρησιμοποιούνται οι πιο επικαιροποιημένες τιμές, δηλαδή τιμές του τρέχοντος βήματος \(t\) για \(j<i\) και τιμές του προηγούμενου βήματος \(t-1\) για \(j>i\):

\[x_i^{\left(t\right)}=\frac{1}{a_{i,i}}\left[b_i-\sum_{j=1}^{j=i-1}a_{i,j}x^{\left(t\right)}_j -\sum_{j=i+1}^{j=n} a_{i,j} x^{\left(t-1\right)}_j \right]\]

4.5.3. Κριτήριο σύγκλισης#

Η ικανή, αλλά όχι αναγκαία, συνθήκη για της σύγκλιση επαναληπτικών μεθόδων είναι η διαγώνια κυριαρχία[1] του πίνακα \(A\), δηλαδή:

\[|a_{i,i}|>\sum_{j=1,j\neq i}^{j=n}|a_{i,j}|, \;\text{για}\; i=1,n\]

4.6. Ειδικές περιπτώσεις#

Στην πράξη συναντάμε συχνά αλγεβρικά συστήματα με ιδιαίτερα χαρακτηριστικά όπως:

συμμετρία
αραιά διατεταγμένα στοιχεία (sparse matrix)
- σε συγκεκριμένες διαγώνιους ή
- σε συγκεκριμένες περιοχές (bands) ή
- σε τυχαίες θέσεις.

Για αυτές τις περιπτώσεις έχουν αναπτυχθεί ειδικές μέθοδοι που επιταχύνουν σημαντικά την επίλυση χρησιμοποιώντας και λιγότερη μνήμη. Ιδιαίτερα χρήσιμος είναι ο αλγόριθμός Thomas που χρησιμοποιείται σε τριδιαγώνια συστήματα.

Ειδικά όταν το σύστημα που επιλύουμε είναι αρκετά μεγάλο, αξίζει να αναζητήσουμε εξειδικευμένες λύσεις. Μια καλή αρχή για την αναζήτηση είναι η numpy.linalg και η πιο αναλυτική scipy.linalg.