Los deze oefening pas op na het oplossen van "Oplossen van Lineaire Stelsels".

Gegeven is de functie $$\verb!gauss_piv(a)!$$ die de NumPy-tabel ($$M$$ rijen en $$N$$ kolommen) $$\verb!a!$$ via Gauss-eliminatie naar een bovendriehoeksmatrix reduceert (via het Gauss-Jordan algoritme met pivoteren). De functie levert als resultaat de getransformeerde matrix, alsook de rang van de originele matrix.

def gauss_piv(a):
    rank = 0
    m, n = a.shape
    i = 0
    for j in range(n):
        p = np.argmax(abs(a[i:m,j]))
        if p>0 :
                a[[i,p+i]] = a[[p+i, i]]
        if a[i,j] != 0:
                rank += 1
                for r in range(i + 1, m):
                    a[r,j:] -= (a[r,j]/a[i,j])*a[i,j:]
                i += 1
        if i >= m:
                break
    return a, rank

Hiervan gebruik makend, kunnen we lineaire stelsels van de vorm $$Ax = b$$ oplossen, met $$A$$ een $$M \times M$$ NumPy-tabel en $$b$$ een $$M \times 1$$ NumPy-tabel. Wanneer we nu stelsels met dezelfde coëfficiëntenmatrix $$A$$ en verschillende kolomvectoren $$b$$ (namelijk $$b_0$$, $$b_1$$, ...), moeten oplossen, kunnen we de reductie van $$A$$ hergebruiken bij het oplossen van de stelsels $$Ax = b_0$$, $$Ax = b_1$$, ... Hiertoe bouwen we een nieuwe matrix $$B$$ die bestaat uit de kolommen $$b_0$$ t.e.m. $$b_{N-1}$$. Het proces van achterwaartse substitutie voeren we dan $$N$$ keer uit.

Schrijf een functie $$\verb!solve_multi(A, B)!$$ met als argumenten de $$M \times M$$ coëfficiëntentabel $$A$$ en de $$M \times N$$ NumPy-tabel $$B$$. De functie levert als resultaat de $$M \times N$$-tabel $$X$$ waarbij de kolom met rangnummer $$i$$ de unieke oplossing voorstelt van het stelsel $$Ax = b_i$$. Je mag hierbij aannemen dat al deze stelsels een unieke oplossing bezitten.
Schrijf een tweede functie $$\verb!inv(A)!$$ die de inverse berekent van een vierkante $$M \times M$$ argumentmatrix $$A$$. Merk hierbij op dat dit gemakkelijk kan door $$\verb!solve_multi()!$$ op te roepen met als 2de argument de eenheidsmatrix. Je mag aannemen dat de matrix $$A$$ inverteerbaar is!

TIPs:

• om de uitgebreide matrix te construeren, kan je gebruik maken van de NumPy-functie $$\verb!np.append()!$$. Raadpleeg de help-functie voor een correct gebruik!
• gebruik de functie $$\verb!np.dot()!$$ om de sommaties bij de achterwaartse substitutie efficiënt te berekenen

Voorbeeld

a = np.array([[1, 2, 3], [7, 10, 12], [9, 21, 3]], dtype=float)
b = np.array([[10, 5], [8, 4], [6, 2]], dtype=float)
solve_multi(a, b) = 
[[-14.63414634  -7.26829268]
 [  5.95121951   2.90243902]
 [  4.24390244   2.15447154]]

a = np.array([[1, 2, 3], [7, 10, 12], [9, 21, 3]], dtype=float)
inv(a) = 
[[-1.80487805  0.46341463 -0.04878049]
 [ 0.70731707 -0.19512195  0.07317073]
 [ 0.46341463 -0.02439024 -0.03252033]]