We wensen het stelsel
$$
\mathbf{Ax} = \mathbf{b}
$$
met $$\mathbf{A}$$ een reële, inverteerbare $$n \times n$$ matrix,
$$
\left[
\begin{array}{ccccc}
a_{00} & a_{01} & a_{02} & \cdots &a_{0(n-1)}\\
a_{10} & a_{11} & a_{12} & \cdots &a_{1(n-1)}\\
a_{20} & a_{21} & a_{22} & \cdots &a_{2(n-1)}\\
\cdots & \cdots & \cdots & \cdots & \cdots\\
a_{(n-1)0} & a_{(n-1)1} & a_{(n-1)2} & \cdots &a_{(n-1)(n-1)}
\end{array}
\right]
$$
$$\mathbf{b}$$ een reële $$n \times 1$$ kolomvector, en $$\mathbf{x}$$ een te zoeken $$n \times 1$$ reële
kolomvector, op te lossen.
Voor diagnonaaldominante matrices $$\mathbf{A}$$, waarvoor geldt dat
$$
|a_{ii}| \ge \sum_{j \ne i} |a_{ij}|
$$
kunnen we dit stelsel oplossen via hieronder geschetste aanpak.
We schrijven $$\mathbf{A}$$ als de som van een diagonaalmatrix $$\mathbf{D}$$ en een matrix $$\mathbf{R}$$
met enkel nullen op de diagonaal, m.a.w. $$\mathbf{A} = \mathbf{D} + \mathbf{R}$$, met
$$ \mathbf{D} =
\left[
\begin{array}{ccccc}
a_{00} & 0 & 0 & \cdots &0\\
0 & a_{11} & 0 & \cdots &0\\
0 & 0 & a_{22} & \cdots &0\\
\cdots & \cdots & \cdots & \cdots & \cdots\\
0 & 0 & 0 & \cdots &a_{(n-1)(n-1)}
\end{array}
\right]
$$
en
$$ \mathbf{R} =
\left[
\begin{array}{ccccc}
0 & a_{01} & a_{02} & \cdots &a_{0(n-1)}\\
a_{10} & 0 & a_{12} & \cdots &a_{1(n-1)}\\
a_{20} & a_{21} & 0 & \cdots &a_{2(n-1)}\\
\cdots & \cdots & \cdots & \cdots & \cdots\\
a_{(n-1)0} & a_{(n-1)1} & a_{(n-1)2} & \cdots & 0
\end{array}
\right]
$$
We bepalen de oplossing van het stelsel iteratief via de constructie van een rij van kolomvectoren
$$\mathbf{x_i}$$, via
$$
\mathbf{x}_{k+1} = \mathbf{D}^{-1}(\mathbf{b} - \mathbf{R}\mathbf{x}_k)
$$
De inverse van de diagonaalmatrix $$\mathbf{D}$$ is uiteraard makkelijk te bepalen!
We starten de iteratie met de nulvector (dus $$\mathbf{x}_0 = \mathbf{0}$$). De rij $${\mathbf{x}_k}$$ convergeert naar
de oplossing van het stelsel $$\mathbf{Ax} = \mathbf{b}$$ indien de matrix $$\mathbf{A}$$ diagonaaldominant is.
Schrijf een functie
solve_diagdom() die het stelstel $$\mathbf{Ax}=\mathbf{b}$$ oplost.
Je mag ervan uit gaan dat de matrix $$\mathbf{A}$$ inverteerbaar maar niet noodzakelijk diagonaaldominant is.
De functie heeft drie argumenten, namelijk:
- de matrix $$\mathbf{A}$$ als een NumPy-tabel van reële getallen ($$n$$ rijen en $$n$$ kolommen)
- de kolomvector $$\mathbf{b}$$ als een NumPy-tabel van reële getallen ($$n$$ rijen, 1 kolom)
- een strikt positief reëel getal $$f$$. Dit getal bepaalt wanneer de iteratie stopt. We stoppen de iteratie zodra geldt dat
$$
|\mathbf{A}\mathbf{x}_k - \mathbf{b}| < f
$$
Hier stelt $$|\mathbf{M}|$$ de norm van de reële matrix $$\mathbf{M}$$ voor, gedefinieerd als de
vierkantswortel uit de som van de kwadraten van alle elementen van $$\mathbf{M}$$.
Het resultaat van de functie is:
- een NumPy-tabel van 1 rij en 1 kolom, die als enig element het getal 0.0 bevat, indien de argumentmatrix $$\mathbf{A}$$ NIET diagonaaldominant is
- de gevonden benadering voor de oplossing van het gegeven stelsel in het andere geval (dus een NumPy-rij met $$n$$ reële getallen)
Merk op: het resultaat van je functie wordt, ten behoeve van de evaluatie in Dodona, naar een lijst omgezet waarbij afgekapt wordt op 6 decimalen.
Het resultaat van je functie MOET wel degelijk een NumPy-tabel zijn met $$n$$ rijen en 1 kolom.
Voorbeeld
A = np.array([[11., 1., 1.], [2., 25., 2.], [30., 3., 52.]])
b = np.array([[1.], [2.], [3.]])
x = solve_diagdom(A, b, 0.0001)
[[ 0.083434]
[ 0.072544]
[ 0.004599]]
#resultaat op 6 decimalen weergegeven)