Lineaire Algebra – Matrices – Numerieke methoden

Dit onderwerp staat niet meer in Lay, maar stond wel in de eerste edities.

We bekijken twee methoden om oplossingen van een stelsel vergelijkingen iteratief te benaderen. We gaan uit van een stelsel vergelijkingen van de vorm \(A\mathbf{x}=\mathbf{b}\) dat uniek oplosbaar is (\(A\) is dus inverteerbaar). Beide methoden zijn gebaseerd op een recursief algoritme van de vorm

\[M\mathbf{x}_{k+1}=N\mathbf{x}_k+\mathbf{b},\quad k=0,1,2,\ldots,\]

waarbij \(A=M-N\). Als de rij \(\{\mathbf{x}_k\}_{k=0}^{\infty}\) convergeert naar een vector \(\mathbf{x}\), dan geldt:

\[M\mathbf{x}=N\mathbf{x}+\mathbf{b}\quad\Longleftrightarrow\quad(M-N)\mathbf{x}=\mathbf{b}\quad\Longleftrightarrow\quad A\mathbf{x}=\mathbf{b}.\]

Dus: \(\mathbf{x}\) is de unieke oplossing van \(A\mathbf{x}=\mathbf{b}\).

De splitsing \(A=M-N\) kan op meerdere manieren worden uitgevoerd. Bij de twee methoden die wij bekijken nemen we aan dat de diagonaalelementen van \(A\) allemaal ongelijk aan nul zijn. Bij de methode van Jacobi is \(M\) een diagonaalmatrix met op de hoofddiagonaal de diagonaalelementen van \(A\) en bij de methode van Gauss-Seidel is \(M\) een benedendriehoeksmatrix waarvan het benedendriehoeksdeel gelijk is aan dat van \(A\).

Problemen die zich hierbij voordoen:

1. De methoden werken niet altijd.


2. De processen convergeren soms tergend langzaam of helemaal niet.


3. Hoe moet men de startvector \(\mathbf{x}_0\) kiezen om het proces (zo snel mogelijk) te laten convergeren?


4. Als men kan kiezen, welke methode werkt dan beter?

We bekijken een voorbeeld en gaan daarna kort in op enkele van deze problemen.

Voorbeeld: Beschouw \(A\mathbf{x}=\mathbf{b}\) met \(A=\begin{pmatrix}5&2&-1\\-1&5&2\\2&1&5\end{pmatrix}\) en \(\mathbf{b}=\begin{pmatrix}20\\5\\3\end{pmatrix}\).

De methode van Jacobi: \(M=\begin{pmatrix}5&0&0\\0&5&0\\0&0&5\end{pmatrix}\) en \(N=\begin{pmatrix}0&-2&1\\1&0&-2\\-2&-1&0\end{pmatrix}\). Dus:

\[M\mathbf{y}=N\mathbf{x}+\mathbf{b}\quad\Longleftrightarrow\quad\left\{\begin{array}{rcrrrrr}5y_1&=&&-2x_2&+x_3&+20\\[2.5mm] 5y_2&=&x_1&&-2x_3&+5\\[2.5mm]5y_3&=&-2x_1&-x_2&&+3.\end{array}\right.\]

Startend met \(\mathbf{x}_0=\mathbf{0}\) vinden we nu: \(\mathbf{x}_1=\begin{pmatrix}4\\1\\3/5\end{pmatrix}\), \(\mathbf{x}_2=\begin{pmatrix}93/25\\39/25\\-6/5\end{pmatrix}\) enzovoort. In decimalen:

\[\mathbf{x}_1=\begin{pmatrix}4.0\\1.0\\0.6\end{pmatrix},\quad\mathbf{x}_2=\begin{pmatrix}3.72\\1.56\\-1.20\end{pmatrix},\quad \mathbf{x}_3=\begin{pmatrix}3.136\\2.224\\-1.200\end{pmatrix},\quad\mathbf{x}_4=\begin{pmatrix}2.8704\\2.1072\\-1.0992\end{pmatrix}, \quad\mathbf{x}_5=\begin{pmatrix}2.93728\\2.01376\\-0.96960\end{pmatrix},\quad\ldots.\]

De methode van Gauss-Seidel: \(M=\begin{pmatrix}5&0&0\\-1&5&0\\2&1&5\end{pmatrix}\) en \(N=\begin{pmatrix}0&-2&1\\0&0&-2\\0&0&0\end{pmatrix}\). Dus:

\[M\mathbf{y}=N\mathbf{x}+\mathbf{b}\quad\Longleftrightarrow\quad\left\{\begin{array}{rrrcrrrrr}5y_1&&&=&&-2x_2&+x_3&+20\\[2.5mm] -y_1&+5y_2&&=&&&-2x_3&+5\\[2.5mm]2y_1&+y_2&+5y_3&=&&&&3.\end{array}\right.\]

Startend met \(\mathbf{x}_0=\mathbf{0}\) vinden we nu: \(\mathbf{x}_1=\begin{pmatrix}4\\9/5\\-34/25\end{pmatrix}\), \(\mathbf{x}_2=\begin{pmatrix}376/125\\1341/625\\-3226/3125\end{pmatrix}\) enzovoort. In decimalen:

\[\mathbf{x}_1=\begin{pmatrix}4.00\\1.80\\-1.36\end{pmatrix},\quad\mathbf{x}_2=\begin{pmatrix}3.00800\\2.14560\\-1.03232\end{pmatrix}, \quad\mathbf{x}_3=\begin{pmatrix}2.93529600\\1.99998720\\-0.97411584\end{pmatrix}, \quad\mathbf{x}_4=\begin{pmatrix}3.00518195200\\1.99068272640\\-1.00020932608\end{pmatrix},\quad\ldots.\]

De exacte oplossing is \(\mathbf{x}=\begin{pmatrix}3\\2\\-1\end{pmatrix}\). We zien dat de methode van Gauss-Seidel (behoorlijk) sneller convergeert dan de methode van Jacobi. Dat is vaak het geval, maar daar staat tegenover dat de methode van Gauss-Seidel wat ingewikkelder is dan die van Jacobi. Er zijn echter ook voorbeelden te construeren waarbij de methode van Jacobi juist sneller gaat. Er zijn ook voorbeelden te construeren waarin één van de methoden niet convergeert of zelfs beide niet convergeren.

Er bestaat echter ook een eenvoudige conditie die garandeert dat beide methoden convergeren. Deze conditie is niet noodzakelijk (ook in andere gevallen kan convergentie optreden), maar in de praktijk kan vaak aan deze voorwaarde worden voldaan. In dat geval is convergentie voor beide methoden dus gegarandeerd. Voor deze conditie is het noodzakelijk dat de matrix \(A\) vierkant is. Een vierkante matrix \(A\) heet strikt diagonaal dominant als de absolute waarde van elk diagonaalelement van \(A\) groter is dan de som van de absolute waarden van de andere elementen in dezelfde rij. Men kan aantonen dat in dat geval de matrix \(A\) inverteerbaar is (\(A\mathbf{x}=\mathbf{b}\) is dus uniek oplosbaar) en dat zowel de methode van Jacobi als de methode van Gauss-Seidel convergeert naar de unieke oplossing voor iedere keuze van de startvector \(\mathbf{x}_0\). Het bewijs van deze bewering laten we hier buiten beschouwing.

Merk op dat de matrix \(\begin{pmatrix}5&2&-1\\-1&5&2\\2&1&5\end{pmatrix}\) strikt diagonaal dominant is, want: \(|5| > |2|+|\pm1|=3\).

Voorbeeld: Beschouw \(A\mathbf{x}=\mathbf{b}\) met \(A=\begin{pmatrix}10&3&2\\3&10&2\\3&2&10\end{pmatrix}\) en \(\mathbf{b}=\begin{pmatrix}15\\15\\15\end{pmatrix}\).

Gemakkelijk is te zien dat \(A\) strikt diagonaal dominant is (want: \(|10| > |3| + |2| = 5\)) en dat de unieke oplossing \(\mathbf{x}=\begin{pmatrix}1\\1\\1\end{pmatrix}\) is.

De methode van Jacobi: \(M=\begin{pmatrix}10&0&0\\0&10&0\\0&0&10\end{pmatrix}\) en \(N=\begin{pmatrix}0&-3&-2\\-3&0&-2\\-3&-2&0\end{pmatrix}\). Dus:

\[M\mathbf{y}=N\mathbf{x}+\mathbf{b}\quad\Longleftrightarrow\quad\left\{\begin{array}{rcrrrrr}10y_1&=&&-3x_2&-2x_3&+15\\[2.5mm] 10y_2&=&-3x_1&&-2x_3&+15\\[2.5mm]10y_3&=&-3x_1&-2x_2&&+15.\end{array}\right.\]

Startend met \(\mathbf{x}_0=\mathbf{0}\) vinden we nu:

\[\mathbf{x}_1=\begin{pmatrix}1.5\\1.5\\1.5\end{pmatrix},\quad\mathbf{x}_2=\begin{pmatrix}0.75\\0.75\\0.75\end{pmatrix},\quad \mathbf{x}_3=\begin{pmatrix}1.125\\1.125\\1.125\end{pmatrix},\quad\mathbf{x}_4=\begin{pmatrix}0.9375\\0.9375\\0.9735\end{pmatrix}, \quad\mathbf{x}_5=\begin{pmatrix}1.03125\\1.03125\\1.03125\end{pmatrix},\quad\ldots.\]

De methode van Gauss-Seidel: \(M=\begin{pmatrix}10&0&0\\3&10&0\\3&2&10\end{pmatrix}\) en \(N=\begin{pmatrix}0&-3&-2\\0&0&-2\\0&0&0\end{pmatrix}\). Dus:

\[M\mathbf{y}=N\mathbf{x}+\mathbf{b}\quad\Longleftrightarrow\quad\left\{\begin{array}{rrrcrrrrr}10y_1&&&=&&-3x_2&-2x_3&+15\\[2.5mm] 3y_1&+10y_2&&=&&&-2x_3&+15\\[2.5mm]3y_1&+2y_2&+10y_3&=&&&&15.\end{array}\right.\]

Startend met \(\mathbf{x}_0=\mathbf{0}\) vinden we nu:

\[\mathbf{x}_1=\begin{pmatrix}1.50\\1.05\\0.84\end{pmatrix},\quad\mathbf{x}_2=\begin{pmatrix}1.0170\\1.0269\\0.98952\end{pmatrix}, \quad\mathbf{x}_3=\begin{pmatrix}0.99402600\\1.00388820\\1.00101456\end{pmatrix},\quad \mathbf{x}_4=\begin{pmatrix}0.998630628000\\1.00020789960\\1.00036923168\end{pmatrix},\quad\ldots.\]

---

Laatst gewijzigd op 12 april 2021