SymPy-demo 1: Matrixalgebra og determinanter#
Demo af Jakob Lemvig, Christian Mikkelstrup, Hans Henrik Hermansen, Karl Johan Måstrup Kristiansen og Magnus Troen. Revideret 22-10-26 af shsp.
from sympy import *
init_printing()
Matricer og ligningssystemer#
Et lineært ligningssystem med én løsning#
Vi vil gerne løse det lineære ligningssystem
Dette kan gøres på adskillige måder - lad os gennemgå nogle metoder herunder.
Metode 0 - Definér alle ligningerne#
En mulighed er simpelthen at definere alle ligningerne og bede SymPy om at løse dem.
x1,x2,x3 = symbols('x1:4')
eq1 = Eq(-x2 + x3, 2)
eq2 = Eq(2*x1 + 4*x2 - 2*x3, 2)
eq3 = Eq(3*x1 + 4*x2 + x3, 9)
eq1, eq2, eq3
Lineære ligningssystemer kan i SymPy løses med linsolve()-kommandoen:
linsolve((eq1,eq2,eq3), (x1,x2,x3))
Systemet har altså præcis en løsning:
Metode 1 - Konvertér systemet til matrixform#
En anden tilgang til at løse et ligningssystem er, når du har systemets koefficientmatrix \(\mathbf A\) og højreside \(\mathbf b\). Da vi allerede har defineret ligningerne ovenfor, kan kommandoen linear_eq_to_matrix() benyttes til at opstille denne matrix og vektor.
A,b = linear_eq_to_matrix([eq1,eq2,eq3],[x1,x2,x3])
A,b
Definér matricer#
Du vil ofte skulle definere din matrix manuelt. Dette gøres med kommandoen Matrix():
A = Matrix([
[0,-1,1],
[2,4,-2],
[3,4,1]])
A
Heri indtastes hver række som en liste med komma mellem listerne. Bemærk, at der også skal tilføjes en omkransende [...] om alle listerne.
En alternativ måde at indtaste elementerne i Matrix()-kommandoen på er følgende, hvor du udfylder det angivne antal rækker og kolonner rækkevis fra venstre mod højre:
Matrix(3,2,[1,2,3,4,5,6])
Definér vektorer#
En vektor defineres med samme kommando. Her udelades blot den omkransende [...], og kun én liste indtastes. Denne liste vil blive fortolket som en søjlevektor:
b = Matrix([2,2,9])
A,b
Altså, vi definerer matricer ved at angive en liste af lister, hvor hver indre liste angiver en række, mens vi definerer søjlevektorer ved kun at angive én liste og ikke andet.
Løs ved brug af matrixformen#
Kommandoen linsolve() fra Metode 0 ovenfor kan også tage matricer som input:
linsolve((A,b))
Metode 2 - Gauss-elimination og reduceret trappeform#
Gauss-Jordan solve#
Har man opskrevet koefficientmatricen og systemets højreside, kan SymPy også løse systemet ved Gauss-Jordan-metoden.
A.gauss_jordan_solve(b)
Denne metode er generelt meget overskuelig fra et brugersynspunkt. Men den resulterer i en fejl, hvis ligningssystemet ikke har nogen løsninger.
A2 = Matrix(2,2, [1,2,0,0])
b_uden_losning = Matrix([1,2])
try:
A2.gauss_jordan_solve(b_uden_losning)
except Exception as e:
print(e)
Linear system has no solution
Reduceret trappeform#
Alternativt kan man opskrive totalmatrixen, som man ville gøre det i hånden:
Det gøres med .hstack(), der sammensætter to matricer/vektorer vandret.
T = Matrix.hstack(A,b)
T
Herefter kan den reducerede trappeform (reduced row-echelon form) findes med .rref().
T.rref(pivots = False)
Igen ser vi, at løsningen er \(x_1 = 4\), \(x_2 =-1\) og \(x_3 =1 \).
pivots = False angiver blot, at vi ikke er interesserede i, hvilke søjler der indeholder pivotelementer. Uden dette argument fås følgende output, som indeholder en smule mere information:
T.rref()
Metode 3 - Simuleret håndregning#
Har man brug for at se hvert trin undervejs i en matrix-reduktion, kan man i stedet benytte SymPys .elementary_row_op() til at udføre rækkeoperationer semi-manuelt.
Herunder ombyttes række 1 og 2, ganges den nye række 1 med \(1/2\) og tilføjes \(-3\) gange række 1 til række 3:
T1 = T.elementary_row_op('n<->m', 0, 1)
T2 = T1.elementary_row_op('n->kn', 0, S(1)/2)
T3 = T2.elementary_row_op('n->n+km',2,-3,0)
T1,T2,T3
Første søjle er hermed færdig, og vi fortsætter med de andre to:
T4 = T3.elementary_row_op('n->kn',1,-1)
T5 = T4.elementary_row_op('n->n+km',2,2,1)
T6 = T5.elementary_row_op('n->n+km',0,-2,1)
T4,T5,T6
T7 = T6.elementary_row_op('n->kn',2,S(1)/2)
T8 = T7.elementary_row_op('n->n+km',0,-1,2)
T9 = T8.elementary_row_op('n->n+km',1,1,2)
T7,T8,T9
Af denne reducerede trappeform aflæser vi den samme løsning som fra de andre metoder ovenfor.
Denne semi-manuelle metode kan være nødvendig, når et system indeholder en ubekendt. Fx, hvis en rækkeoperation undervejs dividerer med et udtryk indeholdende den ubekendte, risikeres division med nul for visse værdier af denne ubekendte. Disse værdier skal håndteres særskilt og kan give helt andre trappeformer og løsninger, men med SymPys .rref() ignoreres de!
Et lineært ligningssystem med flere løsninger#
Vi ønsker at løse det inhomogene ligningssystem, hvis koefficientmatrix er givet ved
A = Matrix([[1,3,2,4,5],[2,6,4,3,5],[3,8,6,7,6],[4,14,8,10,22]])
A
og hvis højreside er givet ved
b = Matrix([9,3,5,32])
b
Metode 1 - Løs direkte i SymPy#
Vi har nu set flere forskellige metoder til løsning af et ligningssystem, metoder som selvfølgelig også kan anvendes på inhomogene systemer. Lad os prøve at benytte nogle af dem.
A.gauss_jordan_solve(b)
linsolve((A,b))
Løsningen indeholder frie parametre (der er uendeligt mange løsninger). SymPy anvender symbolerne \(\tau_0\) og \(\tau_1\) som de frie parametre. De opfylder generelt \(\tau_0, \tau_1 \in \mathbb{F}\).
Ovenstående SymPy-outputs er ikke letlæselige, så lad os opskrive løsninger med LaTex for en pæn præsentation:
Bemærk, at Python anvender nul-indeksering, så en liste af frie parametre som her vil blive nummereret startende fra \(0\).
Metode 2 - Løs vha. rref()#
Vi opskriver først totalmatricen ved:
T = A.row_join(b)
T
Dernæst finder vi den reducerede trappeform:
T.rref()
Vi kan nemt aflæse løsningen til at være:
Vi udelod argumentet pivots=False fra kommandoen denne gang for klart at kunne se, hvilke søjler der indeholder pivotter. Der er ingen pivotter i søjle 3 og 5, hvorfor variablerne \(x_3\) og \(x_5\) bliver vores frie parametre \(\tau_0\) og \(\tau_1\) i dette eksempel.
Tjek#
Lad os tjekke efter, at antallet af frie parametre er korrekt. Vi fandt frem til to i det ovenstående.
Antallet af ubekendte er \(n=5\), så antallet af frie parametre skal gerne være lig med \(n-\rho(\mathbf A)\), hvor \(\rho\) er matricens rang.
A.rank()
5-A.rank()
Som forventet.
Matrixalgebra#
Grundlæggende matrixalgebra og -aritmetik#
Vi vil benytte følgende matricer til eksemplerne i dette afsnit:
A = Matrix([[2,1],[3,0],[7,11]])
B = Matrix([[1,1],[9,3],[-7,-1]])
C = Matrix([[2,1,3],[-6,5,8]])
A,B,C
Vi kan gange matricen \(\mathbf A\) med en konstant ved:
k = symbols('k')
k*A
Da \(\mathbf A\) og \(\mathbf B\) er af samme type, \(\mathbb{R}^{3\times 2}\), er matrixsummen \(\mathbf A+\mathbf B\) defineret.
A+B
Vi kan også udregne linearkombinationer af matricer, for eksempel \(3\mathbf A-5\mathbf B\):
3*A-5*B
Da antallet af søjler i \(\mathbf A\) passer med antallet af rækker i \(\mathbf C\), kan vi bestemme matrixproduktet \(\mathbf A \cdot \mathbf C\).
A*C
Bemærk, at matrixprodukter ikke nødvendigvis kommuterer (vi kan ikke regne med, at \(\mathbf A\mathbf C\) giver det samme som \(\mathbf C\mathbf A\)), hvilket det følgende viser tydeligt:
A*C, C*A
Den transponerede matrix kan findes på flere måder:
A.T, A.transpose(), transpose(A)
Rangen af en matrix kan findes ved:
A.rank()
Hvis en matrix er invertibel, kan dens inverse matrix bestemmes med .inv():
M = Matrix([[2,0,1],[1,2,0],[-2,1,4]])
M.inv(), M*M.inv()
Indeksering i matricer#
Elementer i en matrix kan tilgås ved, at matricen opfattes som et 2D array. Kommandoen M[n,m], hvor M er en matrix, vil give os elementet placeret i række \(n\) og kolonne \(m\). Her er det særligt vigtigt at huske, at Python bruger nulindeksering.
Lad os se på matricen:
M = Matrix([[2,0,1],[1,2,0],[-2,1,4]])
M
Elementet i række 3 og i søjle 3 kan altså fås ved:
M[2,2]
Skal vi have fat i en hel søjle benyttes .col(indeks)-kommandoen. Det følgende vil give dig søjle 2:
M.col(1)
For rækker har vi den tilsvarende .row(indeks).
M.row(0)
Vi kan være endnu mere sofistikerede i vores håndtering af indholdet i en matrix. Hvis vi ønsker at fjerne den første række og den sidste kolonne, så fortæller vi blot SymPy, at vi kun vil have de sidste to rækker og de første to kolonner. Dette gør vi ved at bruge slices som følger:
M[1:, :2]
Læs mere om slicing i Python.
SymPy kan udtrække en snitmatrix fra en matrix med kommandoen .minorMatrix(i,j), hvor i og j er indekserne på den række og søjle, der skal fjernes.
M.minorMatrix(0,-1)
Determinanter#
Vi betragter følgende matrix:
For at finde determinanten her, benytter vi SymPys indbyggede kommando:
A = Matrix([[0,2,3,4],[2,0,4,3],[3,4,0,2],[4,3,2,0]])
A.det()
eller tilsvarende:
det(A)
Determinanten af en snitmatrix kan selvfølgelig findes ved at kombinere de to kommandoer: det(A.minorMatrix(i,j)), men der findes også en direkte kommando til formålet: A.minor(i,j).
det(A.minorMatrix(0,-1)), A.minor(0,-1)
Anvendelser af determinanter#
Eksistens af invers matrix#
Betragt matrixen
Vi ønsker at bestemme \(a \in \mathbb{R}\), således at \(\mathbf A\) er invertibel. Ifølge lærebogen er \(\mathbf A\) invertibel, hvis og kun hvis \(\operatorname{det}(\mathbf A) \neq 0\). Derfor vil vi løse
for \(a\). Dette vil give os alle \(a\)-værdier, for hvilke \(\mathbf A\) ikke er invertibel. Vi udfører dette som følger:
a = symbols('a', real = True)
A = Matrix([[1,2,3],[2,4,1],[3,a,7]])
A
solveset(Eq(det(A), 0), a, S.Reals)
Altså er \(\mathbf A\) invertibel for \(a \in \mathbb{R}\setminus\{6\}\).
Lineær (u)afhængighed#
Betragt følgende fire vektorer i \(\mathbb{C}^4\):
v1 = Matrix([1+I, 3, 0 , 7*I])
v2 = Matrix([2, 4-I, 2*I, 8-I])
v3 = Matrix([3+I, 7-I, 2*I, 8+6*I])
v4 = Matrix([3,-1-I, 7*I, 0])
Det ønskes at undersøge, hvorvidt \(\mathbf v_1,\mathbf v_2,\mathbf v_3,\mathbf v_4\) er lineært uafhængige. Først opstiller vi matricen
A = Matrix.hstack(v1,v2,v3,v4)
A
Vi udregner nu determinanten af denne matrix:
det(A)
Da \(\det(\mathbf A) = 0\), kan vi konkludere, at vektorerne \(\mathbf v_1,\mathbf v_2,\mathbf v_3,\mathbf v_4\) er lineært afhængige (se den relevante sætning i kursuslærebogen).
Lad os nu kun fokusere på \(\mathbf v_1,\mathbf v_2,\mathbf v_4\). Er denne mængde af blot tre af vektorerne lineært uafhængig? Da \(\mathbf B = [\mathbf v_1,\mathbf v_2,\mathbf v_4] \in \mathbb{C}^{4 \times 3}\) ikke er en kvadratisk matrix, kan vi ikke udregne en determinant, og derfor kan vi ikke benytte metoden ovenfor. I stedet vil vi se på deres linearkombination og undersøge, om der findes koefficienter, der ikke er nul, som lader linearkombinationen blive lige nul (se kursuslærebogen for mere).
Lad \(c_1,c_2,c_3 \in \mathbb{C}\) være vilkårlige skalarer, og antag
Dette svarer til et ligningssystem, hvis koefficientmatrix udgøres af de tre vektorer som søjler. Vi løser dette system:
cs = symbols('c1:4')
B = Matrix.hstack(v1,v2,v4)
linsolve((B,zeros(4,1)), cs)
Den eneste løsning er \(c_1 = c_2 = c_3 = 0\). Derfor er \(\mathbf v_1,\mathbf v_2,\mathbf v_4\) lineært uafhængige.
Tips og tricks til matrixkonstruktion#
Nogle foretrækker at opskrive matricer i SymPy som følger:
Matrix(3,3,[1,2,3,4,5,6,7,8,9])
Her angiver (3,3) dimensionen, og rækkerne i matricen indtastes nu som en enkelt liste, [1,2,3,...], fremfor at indtaste en liste per række.
Enhedsmatricer laves nemt i SymPy ved eye(dimension):
eye(3), eye(3,2)
Diagonalmatricer kan defineres ved Matrix.diag(liste-af-værdier):
Matrix.diag([1,2,3,4,5])
Fulde 0-matricer eller 1-matricer kan opnås med zeros(m, n) og ones(m, n):
zeros(2), ones(3,2)
ones(m,n) kan anvendes til at opstille en matrix, hvori alle elementer er ens:
k = 7
k * ones(3,2)
Mere avancerede matrixkonstruktioner kan gøres med Matrix(m, n, lambda i,j: udtryk), en kommando der afsluttes med argumentet udtryk, som angiver, hvad der skal ske med hvert element. For eksempel:
Matrix(3,2, lambda i,j: 3*i + j)
Vi kan fx opstille enhedsmatricen således:
Matrix(3,3, lambda i,j: 1 if i == j else 0)