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

Gebruik makend van deze functie willen we het lineair stelsel $$Ax = b$$ oplossen. Hiertoe bouwen we de matrix die ontstaat door aan de matrix $$A$$ één kolom toe te voegen, namelijk de kolommatrix $$b$$. Via de functie $$\verb!gauss_piv()!$$ krijgen we dan een nieuw stelsel, waarbij de coëfficiëntenmatrix een bovendriehoeksmatrix is. We noteren dit nieuw stelsel (met zelfde oplossing als het originele) als:

$$ Cx = d $$

Met $$C$$ dus een bovendriehoeksmatrix. Hieruit kunnen we de oplossing $$x$$ vinden door het proces van achterwaartse substitutie. Er geldt immers:

$$ C_{n-1,n-1}x_{n-1} = d_{n-1} \\ C_{n-2,n-2}x_{n-2} = d_{n-2} - C_{n-2,n-1}x_{n-1} $$

en meer algemeen dus:

$$ C_{i,i} x_i = d_i - \sum_{j = i+1}^{N-1} C_{ij} x_j $$

Startend met de bepaling van $$x_{n-1}$$ kunnen we dus achtereenvolgens $$x_{n-2}$$, $$x_{n-3}$$ t.e.m. $$x_0$$ vinden. Merk op dat dit slechts kan op voorwaarde dat de rang van de uitgebreide matrix gelijk is aan $$M$$. Is dit niet het geval, dan heeft het stelsel geen unieke oplossing.

Schrijf een functie $$\verb!solve(A, b)!$$ die de unieke oplossing oplevert van het stelsel $$Ax = b$$. Het argument $$A$$ is hierbij een NumPy-tabel met $$M>0$$ rijen en $$M>0$$ kolommen. Het argument $$b$$ is eveneens een NumPy-tabel, met $$M$$ rijen en 1 kolom. Het resultaat is opnieuw een kolomvector, namelijk:

• indien het stelsel een unieke oplossing heeft, dan is die kolomvector de unieke oplossing
• indien het stelsel geen unieke oplossing heeft, dan is die kolomvector de nul-vector, waarij het aantal elementen van de nul-vector gelijk is aan de rang van de uitgebreide matrix.

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], [8], [6]], dtype=float)
solve(a, b) = 
[[-14.63414634]
 [  5.95121951]
 [  4.24390244]]
 
a = np.array([[1, 1, 1], [2, 2, 2], [3, 3, 3]], dtype=float)
b = np.array([[10], [8], [6]], dtype=float)
solve(a, b) = 
[[ 0.]
 [ 0.]]