Voor een reële, vierkante $$N \times N$$ matrix $$\mathbf{A}$$ definiëren we dat de reële, vierkante $$N \times N$$ matrix $$\mathbf{C}$$ een vierkantswortel van $$\mathbf{A}$$ is indien $$\mathbf{C}^2 = \mathbf{A}$$. In deze opdracht gaan we op zoek naar een dergelijk vierkantswortel, waarbij de matrix $$\mathbf{A}$$ diagonaliseerbaar is en bovendien enkel strikt positieve, reële eigenwaarden bezit.

Een matrix $$\mathbf{A}$$ die diagonaliseerbaar is, kunnen we schrijven als:
$$ \mathbf{A} = \mathbf{Q}\mathbf{D}\mathbf{Q^{-1}} $$
waarin

• $$\mathbf{Q}$$ een reële, vierkante $$N \times N$$ ($$N \gt 0$$) matrix is, waarvan de kolommen de eigenvectoren van $$\mathbf{A}$$ voorstellen
• $$\mathbf{D}$$ een reële, vierkante $$N \times N$$ diagonaalmatrix is, waarbij de elementen op de hoofddiagonaal de eigenwaarden van $$\mathbf{A}$$ voorstellen (hierbij hoort de eigenwaarde op kolom $$i$$ van $$\mathbf{D}$$ bij de eigenvector op kolom $$i$$ uit $$\mathbf{Q}$$

Je kan gemakkelijk nagaan dat een wortel uit $$\mathbf{A}$$ dan gegeven wordt door:
$$\mathbf{C} = \mathbf{Q}\mathbf{D}^{1/2}\mathbf{Q}^{-1} $$
Hierin is $$\mathbf{D}^{1/2}$$ de diagonaalmatrix die ontstaat door op de hoofddiagonaal de positieve vierkantswortel te noteren van het element op dezelfde positie uit $$\mathbf{D}$$.

Schrijf een functie wortel() met als enig argument een 2-dimensionale NumPy-tabel (dus van het type array), die een matrix voorstelt. Deze matrix is reëel, vierkant en heeft enkel strikt positieve eigenwaarden. De functie levert als resultaat de vierkantswortel uit die matrix, die je vindt door het bovengeschetste algoritme te volgen. Het resultaat is opnieuw een NumPy-rij (dus opnieuw van het type array).

TIPS:

• je kan een Python-lijst met lengte $$N$$ naar een $$N \times N$$ tabel omzetten via de functie np.diagflat()
• maak voor deze oefening o.a. gebruik van de functies uit de module np.linalg

Voorbeeld

A = np.array([[ 2.80149813,  0.54307116, -0.20599251, -0.23220974],
 [ 1.47565543,  3.42509363, -0.40262172, -0.72659176],
 [ 1.20599251, -0.32771536,  2.67602996,  1.07116105],
 [-1.98127341,  0.78838951,  0.42509363,  1.09737828]])
C = wortel(A)
print(C)
#[[ 1.647659    0.15887501 -0.0551021  -0.06483627]
# [ 0.38213817  1.850458   -0.0996205  -0.2298639 ]
# [ 0.4877458  -0.15476029  1.62008461  0.40300282]
# [-0.81561688  0.32728355  0.15586579  1.02806276]]