Python, en met name het Symbolic Python pakket afgekort met SymPy en te importeren als sympy, heeft voldoende ingebouwde faciliteiten om stelsel lineaire vergelijkingen numeriek en algebraisch op te lossen. We illustreren ze aan de hand van het stelsel

$$\left\{\begin{array}{rrrrrl} x & + & 2y & + & 3z & =\;2\\ & & y & + & 2z & =\;1\\ 3x & + & y & + & z & =\;3 \end{array}\right.$$ één oplossing, namelijk $$\cv{x\\y\\z}= \left(\!\!\begin{array}{r} 1\\-1\\1\end{array}\right)$$

We zullen SymPy hier niet volledig behandelen, maar je kunt hier computer algebra mee bedrijven, dat wil zeggen, niet alleen rekenen met getallen, maar ook met symbolen. Door zelf met de pakketten te stoeien en op de officiële website www.sympy.org de documentatie te raadplegen kom je al een heel eind.

We beginnen met het SymPy pakket te importeren met een afgekorte naam en het stelsel van vergelijking in te toetsen in een Python sessie. In sympy schrijf je een vergelijking niet op met = of ==, maar met de uitdrukking Eq(x,y) voor $x=y$. Handiger is om de vergelijking te schrijven als $x-y = 0$ en dan gewoon x-y te gebruiken. Je moet ook (net als in MATLAB) expliciet aangeven welke symbolen je wilt gebruiken. De bij een stelsel lineaire vergelijkingen passende matrixvorm te bepalen met het commando linear_eq_to_matrix. Let op dat het pakket sympy zijn eigen datastructuur voor matrices heeft. Voor grote stelsels is het beter om het numpy pakket te gebruiken, maar dan moet je de lineaire algebra theorie wel kennen om te weten wat je doet en vlot te kunnen werken.

>>> import sympy as sy
>>> x, y, z = sy.symbols('x y z')
>>> var = x,y,z
>>> sys = [x + 2*y + 3*z - 2, y + 2*z - 1, 3*x + y + z - 3]
>>> A, b = sy.linear_eq_to_matrix(sys, var)
>>> A
Matrix([
[1, 2, 3],
[0, 1, 2],
[3, 1, 1]])
>>> b
Matrix([
[2],
[1],
[3]])

We kunnen nu linsolve gebruiken om de vector $\vec{v}$ te vinden zodanig dat $A\vec{v}=\vec{b}$

>>> opl = sy.linsolve((A,b), var); opl
{(1, -1, 1)}

De oplossing van het oorspronkelijke stelsel hadden we ook kunnen vinden hadden we ook via solve kunnen vinden zonder op de matrixvorm over te gaan.

>>> opl = sy.solve(sys,var); opl
{x: 1, y: -1, z: 1}
>>> opl[x], opl[y], opl[z]
(1, -1, 1)

Voordeel van de matrixvorm is dat je de theorie voor het oplossen van vergelijkingen beter in de hand hebt. Zo kun je het stelsel oplossen bijvoorbeeld Gauss-eliminatie toepassen om het stelsel op te lossen. Dit kan impliciet via de functie solve_linear_system of expliciet via de methode rref, waarvan de naam een acroniem is voor row-reduced echelon form, iets wat wij de gereduceerde trapvorm noemen.

>>> Ab = A.row_join(b); Ab   # de gerande matrix
Matrix([
[1, 2, 3, 2],
[0, 1, 2, 1],
[3, 1, 1, 3]])
>>> sy.solve_linear_system(Ab, x, y, z)
{x: 1, y: -1, z: 1}
>>> GT = Ab.rref(); GT[0]   # berekening van de gereduceerde trapvorm
Matrix([
[1, 0, 0, 1],
[0, 1, 0, -1],
[0, 0, 1, 1]])
>>> GT[0][:, 3]  % selectie van de vierde kolom
Matrix([
[ 1],
[-1],
[ 1]])

$\phantom{x}$

Zoals gezegd kun je bij numerieke berekeningen beter het NumPy pakket gebruiken, maar dan moet je zelf alles in matrixberekeningen organiseren. Het voorbeeld kan er dan als volgt uitzien. Merk op dat het NumPy pakket, net als zijn grote broer SciPy (een afkorting voor Scientific Python), geen ingebouwde commando's heeft voor het berekenen van de gereduceerde trapvorm en dat we daarom switchen tussen dit pakket en SymPy om deze klus te klaren.

>>> import numpy as np
>>> A = np.array([[1, 2, 3], [0,1,2], [3, 1, 1]]); print(A)
[[1 2 3]
 [0 1 2]
 [3 1 1]]
>>> b = np.array([[2],[1],[3]]); print(b)
[[2]
[1]
[3]]
>>> opl = np.linalg.solve(A,b); print(opl)
[[ 1.]
[-1.]
[ 1.]]
>>> Ab = np.hstack((A,b))  # gerande matrix
>>> print(Ab)
[[1 2 3 2]
[0 1 2 1]
[3 1 1 3]]
>>> import sympy as sy
>>> Abmat = sy.Matrix(Ab)
>>> GT = np.array(Abmat.rref()[0]); print(GT)  # gereduceerde trapvorm via sympy
[[1 0 0 1]
[0 1 0 -1]
[0 0 1 1]]