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

4.1. Διατύπωση του προβλήματος

Ένα σύστημα \(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}\]

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

• \(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}\]

ή ισοδύναμα:

\[Α x =b \]

4.2. Απαλοιφή Gauss

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

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

Ορισμός

Η γραμμή ή εξίσωση, που χρησιμοποιούμε για να μηδενίσουμε τις τιμές κάτω από την διαγώνιο, ονομάζεται οδηγός (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} \)

Solution to Exercise 4.1

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

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

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

\[ \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}\]

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

\[ \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}\]

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

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

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

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

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

Exercise 4.2

Γράψτε μια συνάρτηση που επιλύει ένα γραμμικό σύστημα \(n\) αγνώστων με την μέθοδο Gauss χωρίς οδήγηση. Η συνάρτηση θα δέχεται ως παραμέτρους τους πίνακες Α και b και θα επιστρέφει τον πίνακα x. Χρησιμοποιήστε την συνάρτηση για την επίλυση του γραμμικού συστήματος της Άσκησης 4.1.

Solution to Exercise 4.2

import numpy as np


def solveGauss(A: np.ndarray, b: np.ndarray) -> np.ndarray:
    """
    Solve Ax=b with Gauss elimination.
    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.]

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

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

\[\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}&=&\underbrace{\frac{n^2}{2}}_\text{Σημαντικός όρος}+\underbrace{\mathcal{O}(n)}_\text{Αμελητέοι όροι}\\ \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{Σημαντικός όρος}+\underbrace{\mathcal{O}(n^2)}_\text{Αμελητέοι όροι} \end{align} \end{split}\]

Important

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

Προσέξτε ότι, όταν αλλάξει το άνω όριο του αθροίσματος σε \(n-2\), δεν επηρεάζεται ο σημαντικός όρος:

\[ \sum_{i=1}^{n-2} i^2 = 1 + 2^2 + 3^2 + \dots + (n-2)^2 = \frac{(n-2)(n-1)(2(n-2)+1)}{6}=\frac{n^3}{3}+\mathcal{O}(n^2) \]

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

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

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

Exercise 4.3

Ένας υπολογιστής χρειάστηκε 0.05 ms για την επίλυση ενός γραμμικού συστήματος 1000 αγνώστων με την μέθοδο Gauss. Εκτιμήστε πόσο χρόνο (ms) θα χρειαστεί για την επίλυση ενός συστήματος 2000 αγνώστων με την ίδια μέθοδο. Μπορείτε να αγνοήσετε εξαρτήσεις μικρότερης τάξης για απλοποίηση των πράξεων. Θεωρήστε ότι ο χρόνος υπολογισμού αυξάνει γραμμικά με το πλήθος των πράξεων.

Solution to Exercise 4.3

\[\begin{split} t_1=0.05 \text{ms}\\ n_1=1000\\ n_2=2000 \end{split}\]

Το πλήθος των πράξεων μεγαλώνει με την τρίτη δύναμη του \(n\) στην μέθοδο Gauss:

\[\frac{N_2}{N_1}=\left(\frac{n_2}{n_1}\right)^3\]

και ο χρόνος υπολογισμού αυξάνει γραμμικά με το πλήθος των πράξεων:

\[\frac{t_2}{t_1}=\frac{N_2}{N_1}\]

επομένως:

\[ t_2={t_1}\frac{N_2}{N_1}={t_1}\left(\frac{n_2}{n_1}\right)^3=0.05\left(\frac{2000}{1000}\right)^3= 0.40 \text{ms} \]

4.4. Απαλοιφή 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:

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

2. Μετά την διαδικασία απαλοιφής των στοιχείων κάτω από την διαγώνιο, ακολουθεί η διαδικασία απαλοιφής των στοιχείων πάνω από την διαγώνιο.

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

Exercise 4.4

Επιλύστε το γραμμικό σύστημα της Άσκησης 4.1 με την απαλοιφή Gauss-Jordan.

Solution to Exercise 4.4

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

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

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

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

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

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

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

\[ \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}\]

2. Τα στοιχεία πάνω από την διαγώνιο στην τελευταία στήλη μηδενίζονται όπως στην απαλοιφή 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}\]

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

\[ \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.5. Παραγοντοποίηση LU

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

\[Α=LU\]

όπου

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

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

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

\[Lm=b\]

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

\[Ux=m\]

Το μόνο που λείπει πλέον είναι να βρεθούν οι πίνακες \(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\) έχει μονάδες στην διαγώνιο, ενώ τα στοιχεία κάτω από την διαγώνιο είναι οι συντελεστές που χρησιμοποιήθηκαν στην απαλοιφή. Για παράδειγμα το στοιχείο \(a_{2,1}\) απαλείφθηκε με την εξίσωση:

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

επομένως \(l_{2,1}=1/4\). Ο πλήρης πίνακας \(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 Παραγοντοποιεί τον πίνακα Α σε P, L και U, όπου P o πίνακας μεταθέσεων.

• scipy.linalg.lu_factor Παραγοντοποιεί τον Α και επιστρέφει τους πίνακες lu και piv (πακεταρισμένη μορφή των P,L,U).

• scipy.linalg.lu_solve Επιλύει το σύστημα με τους πίνακες lu, piv και b.

Exercise 4.5

Επιλύστε το γραμμικό σύστημα της Άσκησης 4.1 με παραγοντοποίηση LU χρησιμοποιώντας:

1. την συνάρτηση numpy.linalg.solve.

2. την συνάρτηση scipy.linalg.lu_solve και όποια άλλη συνάρτηση από την ίδια βιβλιοθήκη χρειάζεται.

Solution to Exercise 4.5

from scipy.linalg import lu, lu_factor, lu_solve
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)

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

# P,L,U are packed in lu (combination of L and U) and piv (pivot table)
lu, piv = lu_factor(A)
# print(lu)
# print(piv)
x = lu_solve((lu, piv), b)
print(x)

[ 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.]]
[ 1. -0.  2.]

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

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

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

4.6.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.6.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.6.3. Κριτήριο σύγκλισης

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

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

Exercise 4.6

Ελέγξτε αν ο πίνακας \(Α\) της Άσκησης 4.1 έχει διαγώνια κυριαρχία.

Solution to Exercise 4.6

• Για \(i=1\), \(|4|>|2|+|1|\), ισχύει.

• Για \(i=2\), \(|-1|>|1|+|2|\), δεν ισχύει.

• …

Επομένως ο πίνακας \(Α\) δεν έχει διαγώνια κυριαρχία.

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

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

1. συμμετρία

2. αραιά διατεταγμένα στοιχεία (sparse matrix)

   • σε συγκεκριμένες διαγώνιους (bands) ή

   • σε συγκεκριμένες περιοχές ή

   • σε τυχαίες θέσεις.

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

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

---

[1]

Μια διαφορετική διατύπωση της διαγώνιας κυριαρχίας δίνεται με το κριτήριο του Scarborough.