1. No category

HH/ITE/BN
Matriser och Mathematica
1
Något om Matriser och Mathematica
Bertil Nilsson
2015-08-15
Denna väntan varje jul  xi

1
1 1
 x mod 1
1 2 i
x0
x1
x2
x3
x4
x5
x6
x7
x8
x9
x10
x11
x12
x13
x14
x15
x16
x17
x18
x19
x20
x21
x22
x23
x24
x25
x26
x27
x28
x29
x30
x31
x32
x33
x34
x35
x36
x37
x38
x39
x40
x41
x42
x43
x44
x45
x46
x47
x48
x49
x50
x51
x52
x53
x54
x55
x56
x57
x58
x59
x60
x61
x62
x63
x64
x65
x66
x67
x68
x69
x70
x71
x72
x73
x74
x75
x76
x77
x78
x79
x80
x81
x82
x83
x84
x85
x86
x87
x88
x89
x90
x91
x92
x93
x94
x95
x96
x97
x98
x99
x100
x101
x102
x103
x104
x105
x106
x107
x108
x109
x110
x111
x112
x113
x114
x115
x116
x117
x118
x119
x120
x121
x122
x123
x124
x125
x126
x127
x128
x129
x130
x131
x132
x133
x134
x135
x136
x137
x138
x139
x140
x141
x142
x143
x144
x145
x146
x147
x148
x149
x150
x151
x152
x153
x154
x155
x156
x157
x158
x159
x160
x161
x162
x163
x164
x165
x166
x167
x168
x169
x170
x171
x172
x173
x174
2
Matriser och Mathematica
HH/ITE/BN
ť Förord
På följande sidor presenteras en elementär "streetwise guide" till matriser med flitig användning av Mathematica. Framställningen är
fåordig, fri från pedanteri men i någon mening fullständig. Det man väsentligen behöver veta om begrepp, terminologi, beteckningar
och teori för att modellera och lösa problem i framtida kurser och yrkesliv som ingenjör, naturvetare eller lärare klarläggs och
typiska exempel ges.
ť Matrisbegreppet
En uppsättning av element, oftast vanliga tal, ordnade "rektangulärt" i rader och kolonner kallas för en matris. En sådan betecknas
nästan alltid med stora "feta" bokstäver, t.ex. . Om en matris innehåller m rader och n kolonner säger vi att den är av "typen m
gånger n" och skriver typ
m n. Inte sällan brukar man dekorera matrisens namn med ett index som förtydligar dess typ, m n .
Elementen i matrisen anges sedan med motvarande "lilla" bokstav och ett index som talar om dess position. Exempelvis har vi en
matris med m rader och n kolonner
m n
aij m
a11 a12
a21 a22
a1n
a2n
am1 am2
amn
n
Elementet aij är elementet i i:te raden och j:te kolonnen räknat uppifrån respektive från vänster. Man säger att element aij befinner
sig på plats i, j . Ofta ser man att elementet också refereras med skrivsättet
ij
vilket är inspirerat av att matrisen ibland skrivs ut i
full form med hakparenteser [] istället för runda ().
Exempel: För matrisen
2
1
7
3 2
gäller att typ
5
3
2
3 2 eftersom den har 3 rader och 2 kolonner. Vidare är exempelvis elementen a21
Definition. Två matriser
nA nB . Detta skrivs typ
och för vilka typ
= typ .
mA nA och typ
Exempel: Antag att och är matriser sådana att typ
är inte av samma typ, det vill säga typ
typ .
1 och a32
2.
mB nB , säges vara av samma typ om mA
4 2 och typ
2 4. Då har båda matriserna 4 2
mB och
8 element, men de
Om en matris har lika många rader som kolonner, det vill säga n n , kallas den av naturliga skäl kvadratisk. Sådana är mycket
vanliga i tillämpningar, vilket vi får anledning att återkomma till. En matris som består av endast en rad kallas av uppenbara skäl för
en radmatris eller radvektor och en matris som består endast av en kolonn för en kolonnmatris eller kolonnvektor. Slutligen
faller det sig då naturligt att ibland tolka en matris av typen 1 1 som en skalär eller ett tal. Denna dualitet med lite liberalare syn på
typerna gör att räkningarna "flyter" på lite enklare i löpande text. Men när man verkligen genomför beräkningar måste man vara på
sin vakt.
I Mathematica skrivs feta bokstäver, t.ex. A eller , med "fet" font respektive dsA , där ds står för double-struck. När det gäller
textceller i Mathematica brukar författaren blanda efter dagsnotering medan det i inputceller alltid är den senare som gäller. Det går
naturligtvis lika bra att använda vanliga tecken! Internt i Mathematica representeras en matris som en vektor med vektorer men
presenteras på den form vi känner igen från ovan. För att mata in en matris föreslås att man börjar med ny matris i paletten varefter
man mekar om till rätt typ och fyller i sina element enligt nedan. Samma gäller naturligtvis för redigering av redan befintlig matris,
se nedan. Om man inte får snygga 
 i en output cell utan
,
,
,
kan det vara läge att ändra i Edit>Preferences...>Eval-
uation>Format type of output cells>TraditionalForm.
Ny matris från palette: 

Ny rad under markörens läge:
Ny kolonn till höger om markörens läge:
Hoppa till nästa :
͓ ͔
͓,
Ny rad ett:
ställ markören till vänster om a11 och
Ny kolonn ett:
ställ markören till vänster om a11 och
Ta bort hel rad kolonn: svärta med markören, sedan
͓ ͔
͓,
HH/ITE/BN
Matriser och Mathematica
3
Så när som vissa funktionsnamn ansluter sig Mathematica väldigt smidigt till definitionerna i matrisalgebra.
4 3 1
;
1 5 2
Dimensions
2, 3
Att välja ut ett eller flera element ur matrisen brukar kallas att indexera. Detta anges med dubbla hakparenteser [[...]]. De två
hakparenteserna kan göras lite kompaktare med
och
, vilket använts i exemplen nedan.
2, 3
Hämta element a23
2, All
Hämta rad 2
2
1, 5, 2
All, 3
Hämta kolonn 3
1, 2
2, 3, 1
Hämta delmatris
a23 a21
i denna ordning
2, 1
1, 2 , 3, 1

1
2
Hämta delmatris 
a13 a11

a23 a21
4

1
1, 3
7
Tilldela element a13 nytt värde
7

4 3 7

1 5 2
100 200
300 400
1, 2 , 3, 2


Tilldela delmatrisen 
a12 a13
 nya värden
a22 a23
100 200

300 400
4 200 100

1 400 300
ť Likhet för matriser
Definition. Två matriser och säges vara lika, vilket skrivs
, om de är av samma typ och elementen på motsvarande
platser i de båda matriserna är lika, det vill säga aij bij för alla förekommande i och j. Om dessa villkor inte är uppfyllda säger
vi att matriserna är olika och skriver
.
Exempel: Matriserna
2 3
2 5
1 3
uppfyller uppenbarligen kraven för likhet så vi kan skriva
2 5
1 3
True
7
2
2 5
1 3
7
2
7
2
och
2 3
2 5
1 3
7
2
. I Mathematica används naturligtvis två likhetstecken!
4
Matriser och Mathematica
HH/ITE/BN
ť Nollmatris
Definition. Om samtliga element i en matris är 0 kallas matrisen för nollmatris och reserverar namnet
finns oändligt med nollmatriser, en för varje typ.
för denna. Notera att det
Exempel: Matriserna
3 2
är exempel på nollmatriser, men
3 2
2 2
0 0
0 0
0 0
och
0 0
0 0
2 2
eftersom de är av olika typ.
ť Addition och subtraktion av matriser
Definition. Antag att
för vilken typ
typ
och är två matriser av samma typ. Med summan av
typ och cij aij bij för alla förekommande i och j.
och , betecknat
, menas den matris
Addition sker alltså elementvis vilket kräver att matriserna måste vara av samma typ. Så med
a11 a12
a21 a22
a1n
a2n
och
m n
am1 am2
b11 b12
b21 b22
b1n
b2n
bm1 bm2
bmn
m n
amn
har vi
a11
a21
m n
m n
b11
b21
a12
a22
b12
b22
a1n
a2n
b1n
b2n
m n
m n
am1
bm1 am2
bm2
amn
m n
bmn
Exempel: Följande kalkyl går tydligen bra. Inte fel att använda stödlinjer när man "räknar"...
2 5 4
1 3 0
2 5 4
1 3 0

1 2
11 8
1 2
11 8
2
2
7
1
1
5
11 3
2 4
8
2
0
7
3 7 2
10 11 7
2
7
3 7 2

10 11 7
Räknelagar: Om
1
2
3
, ,
och
är matriser av samma typ gäller
kommutativa lagen
associativa lagen
Subtraktion av matriser får vi nu odramatiskt genom att snegla på addition ovan och byta alla "+" till " ".
Definition. Antag att och är två matriser av samma typ. Med subtraktion mellan och , betecknat
, menas den
matris för vilken typ
typ
typ och cij aij bij för alla förekommande i och j. Som vanligt är subtraktion inte kommutativ, det vill säga
.
a11
a21
m n
m n
b11
b21
a12
a22
b12
b22
a1n
a2n
b1n
b2n
m n
m n
am1
bm1 am2
bm2
amn
bmn
m n
HH/ITE/BN
Matriser och Mathematica
5
Exempel: Följande kalkyl går tydligen bra
2 5 4
1 3 0
2 5 4
1 3 0

1
12
3
5
1 2
11 8
1 2
11 8
2
2
7
1
1
5
11 3
2 4
8
2
0
1
12
7
3
5
6
7
2
7
6

7
ť Multiplikation av matris med ett tal
Definition. Med produkten av talet s och matrisen
s ij s ij för alla förekommande i och j.
Detta betyder helt enkelt att varje element i
m n
s
m n
, betecknat s , menas den matris för vilken typ s
multipliceras med s och resultatet ärver typen från
s
a11 a12
a21 a22
a1n
a2n
sa11 sa12
sa21 sa22
sa1n
sa2n
am1 am2
amn
sam1 sam2
samn
typ
och
.
m ns
Exempel: Dé blir bara svårare och svårare
3
3
2 5 4
1 3 0

6 15 12

3 9 0
Räknelagar: Om
1 s
s
2 s t
s
3 st
ts
,
s
t
st
2 5 4
1 3 0
3 2
3
1
3 5 3 4
3 3 3 0
6 15 12
3 9 0
är matriser av samma typ och s och t godtyckliga tal gäller
distributiva lagen
distributiva lagen
associativa lagen
ť Multiplikation av matriser
Matrismultiplikation är något mer komplicerad än de nästan självklara operationer vi har sett hittills. Den är dock vald med omsorg
för att vara så användbar som möjligt i diverse tillämpningar.
Definition. Multiplikation av och , skrivet
utan multiplikationstecken, är tillåtet endast då antalet kolonner i är lika
med antalet rader i . Resultatet blir då en ny matris som ärver antalet rader från och antalet kolonner från . Alltså matrismultiplikationen
p
m p
q n
m n
är utförbar endast om p
q. Elementet cij
aik bkj
k 1
Elementet cij fås alltså som skalärprodukten mellan den i:te raden i
och den j:te kolonnen i , vilket förklarar att antalet kolonner i
måste vara lika med antalet rader i . Matrismultiplikation är i allmänhet inte kommutativ, det vill säga
nen
säger vi att eftermultipliceras av och omvänt att förmultipliceras av .
. I konstruktio-
Exempelvis bestämning av elementet c32 vid eftermultiplikation av en 4 3-matris med en 3 2-matris. Lägg märke till typerna för
såväl och som resultatmatrisen .
6
Matriser och Mathematica
b12
b22
b32
a31 a32 a33
4 3
HH/ITE/BN
a31 b12
a32 b22
a33 b32
3 2
4 2
Exempel: Låt
6 4
3 5
2 3
1
2
och
7
5
3
3 2
4
8
6
så har vi
2 2
2 3
6 4
3 5
3 2
6 7
1
2
2
4 5
3 7
5 5
7
5
3
3
1 3
4
8
6 32
6 4 4
2 3
3 4
8
1 6
5
8
2 6
2 2
59
10
14
40
3 3
30
54
0
48
20
42
2 2
eller
3 3
3 2
7
5
3
2 3
4
6 4 1
8
3 5 2
6 32
7 6 4 3
7 4
5 6
3 6
Vi ser att
8
6
3
3
5 4
3 4
2 3
4 5
7
8 5 5
6 5
1
4 2
1
6 2
1
3
8 2
1
21
9
3 3
trots att båda multiplikationerna i just detta fall är tillåtna. De får ju inte ens samma typ.
Matrismultiplikation i Mathematica anges med vanlig punkt likt skalärprodukt mellan vektorer (som ju faktiskt är lite inblandad).
Observera att punkten inte får utelämnas eftersom detta tolkas som elementvis listmultiplikation, det vill säga cij aij bij . Vi provkör
med matriserna i föregående exempel.
6 4
3 5
1
;
2
7
5
3
4
8 ;
6
.

59
10
14

40
.
30
54
0
48
20
42
1
21
9
Avslutningsvis poängteras ännu en gång att faktorernas ordning i produkten
är väsentlig, det vill säga
är något helt annat än
. Eftersom multiplikation av matriser är definierad endast om antalet kolonner i den första matrisen är lika med antalet rader i den
andra så är det inte säkert att produkten
är definierad även om
är det. Om både
och
är definierade behöver de ändå
inte vara lika. Ett fall då båda produkterna
och
alltid är definierade är då och båda är kvadratiska matriser av samma
typ.
behöver dock inte heller i detta fall vara lika med
.
Räknelagar: Om ,
1
C
2
3
och
är matriser gäller under förutsättning att ingående operationer är tillåtna
distributiva lagen
distributiva lagen
associativa lagen
HH/ITE/BN
Matriser och Mathematica
7
ť Matristransponering
Definition. Den matris som bildas av en given matris om man låter rader och kolonner byta roller (det vill säga raderna görs
till kolonner och tvärtom), kallas transponatet till och betecknas med , vilket utläses "A-transponat". Om typ
m n, så
är alltså typ
n m.
I Mathematica finns Transpose, vilket också kan skrivas lite kompaktare som
" " som resultat. Detta ska inte förväxlas med " upphöjt till T"; T .
tr
direkt efter matrisen med ett litet dekorativt
2 5 4
1 3 0
2
5
4
1
3
0
I ingenjörslitteratur ser man ofta skalärprodukt mellan två vektorer skriven som
. Detta grundar sig på att vektorer då ömsom
betraktas som kolonnmatriser. Exempelvis får vi med god kontroll över ingående typer
2
1
3
3 1
1
1
2
2 1 3
3 1
1 3
1
1
2
2 1
1
1
3 2
1 1
7
1 1
3 1
För att det ska bli just en skalärprodukt är det nu naturligt att betrakta den sista matrisen med en rad och en kolonn som en skalär, det
vill säga ett tal. Det finns uppenbara risker med denna liberala syn på typbegreppet och det finns datorprogram som tillåter det.
Mathematica tillhör inte dessa utan gör alltid en mycket rigorös typkontroll inför varje operation. För den seriöse användaren är
detta inget problem utan en trygghet. Använd indicering för att hämta ut önskat tal.
2
1
3
aTb
1
1
2
.
7
tal
aTb 1, 1
7
För att transponera en vektor, vilket egentligen är omöjligt eftersom den bara har en dimension, måste den i Mathematica först göras
till en radmatris, det vill säga en matris med en rad. Eftersom en matris representeras internt som en vektor med vektorer är det bara
att "stänga" in vektorn inom {}.
InputForm
1 2 3

3 4 5
1, 2, 3, 3, 4, 5
1, 2, 3
1 2 3
1, 2, 3
1
2
3
Räknelagar: Om
1
2
och
är matriser och s ett tal gäller
3 s
s
4
Notera ordningen i högerledet
ť Kvadratisk matris
Om en matris har lika många rader som kolonner, det vill säga n n , kallas den inte oväntat kvadratisk. För sådana matriser är det
lite "tårta på tårta" att säga typ n n, så inte sällen säger man istället ordning n. Då och då brukar man också höra ordning som
synonym för typ, vilket inte brukar vålla någon förvirring. Kvadratiska matriser är mycket vanligt förekommande i tillämpningar,
vilket vi får anledning att återkomma till.
8
Matriser och Mathematica
HH/ITE/BN
I en kvadratisk matris kallas elementen aii för huvuddiagonalen i . Om alla element under huvuddiagonalen är noll kallas den
högertriangulär eller övertriangulär och vänstertriangulär eller undertriangulär om alla element över huvuddiagonalen är noll.
2 5 9
Exempel: Matrisen 0 7 1 är högertriangulär och
0 0 3
7
1
8
1
0
5
9
3
0
0
4
5
0
0
är vänstertriangulär.
0
3
Om det för en kvadratisk matris gäller att
kallas matrisen symmetrisk och antisymmetrisk om
. Diagonaleleaii
2aii 0
aii 0. Varje kvadratisk matris kan skrivas som summan
menten i en antisymmetrisk matris är alltid noll, ty aii
av en symmetrisk och en antisymmetrisk matris. Produkten av en godtycklig matris och dess transponat är symmetrisk, ty
.
Exempel: För godycklig matris
inte ens blir av samma typ.
2 5 4
.
1 3 0

är alltid
och
symmetriska matriser, dock inte nödvändigtvis lika eftersom de vanligtvis
2 5 4
1 3 0
45 13

13 10
2 5 4
1 3 0
.
2 5 4
1 3 0
5 7 8
7 34 20
8 20 16
Summan av talen i huvuddiagonalen kallas för spåret (eng. trace). Tydligen har vi för
i senaste Exemplet.
Tr
55
Huvuddiagonalen själv hämtas med en liten variant av samma funktion.
Tr
, List
5, 34, 16
ť Diagonalmatris och enhetsmatris
Om en matris är kvadratisk och alla element utom de på huvuddiagonalen är noll kallas matrisen för en diagonalmatris. Ett
.
vanligt sätt att ange en sådan är diag a11 , a22 ,
DiagonalMatrix

3, 6
3 0

0 6
Om en diagonalmatris har enbart ettor på huvuddiagonalen kallas den enhetsmatris. För denna brukar man reservera namnen
(identity matrix) eller (enhetsmatris). Namnet kommer sig av att den tjänar ungefär som en etta gör vid "vanlig" räkning. Notera
liksom för nollmatris att det finns oändligt med sådana här matriser, en för varje typ.
IdentityMatrix 2

1 0

0 1
Exempel: Effekten av att förmultiplicera en matris
element i .
.

1 1 1 1
1 1 1 1
3 3 3 3

6 6 6 6
med en diagonalmatris
är att multiplicera varje rad i
med motsvarande
HH/ITE/BN
Matriser och Mathematica
Effekten av att eftermultiplicera en matris
element i
1
1
1
1
3
3
3
3
med en diagonalmatris
är att multiplicera varje kolumn i
9
med motsvarande
1
1
.
1
1
6
6
6
6
ť Potenser
Kvadratiska matriser är de enda matriser man kan utsätta för upprepad multiplikation med sig själv. Om vi låter vara enhetsmatrisen definierar vi naturligt potenser som 0 , 1
, 2
, 3
och så vidare. Eller mera formell rekursiv definition
om n
n
n 1
0
om n är ett positivt heltal
Negativa eller brutna exponenter befattar vi oss inte med. Naturligtvis gäller de vanliga potenslagarna
m n
mn
.
m
n
m n
och
Exempel: När det gäller potenser av matriser så skriv inte upphöjt till i Mathematica. Detta tolkas då som att varje element i
matrisen skall upphöjas till angiven exponent.
2
4

2
4
3
5
3

5
2

4 9

16 25
.

16
28
21

37
Då exponenten n blir stor finns funktionen MatrixPower som tröst.
. . . . . . . . . . . . . .

2 549 464 720 604
4 482 738 450 156
MatrixPower

3 362 053 837 617

5 911 518 558 221
, 15
2 549 464 720 604
4 482 738 450 156
3 362 053 837 617

5 911 518 558 221
ť Matrisvärd funktion av en variabel
I tidigare kurs har vi stiftat bekantskap med en funktion av en reell variabel. Vi har sett att detta begrepp kan enkelt spilla över till en
vektorvärd funktion av en reell variabel. Då är steget litet till en matrisvärd funktion av en reell variabel. Vi låter helt enkelt varje
element i matrisen vara en funktion av samma variabel, säg
x
a11 x
a21 x
a12 x
a22 x
Definition av derivation och integration gör knappast någon förvånad.
10
Matriser och Mathematica
a11 x
x
a21 x
x
x
x
a12 x
x
a22 x
x
x
x
HH/ITE/BN
a11 x
x
a12 x
x
a21 x
x
a22 x
x
Naturligtvis kan man ha bestämd integral om så önskas.
Exempel: Mathematica är naturligtvis kunnig som vanligt
x2
2x
D
Sin x
Log x
Cos x
ArcTan x
;
,x
2x
cos x
2x
2
sin x
1
1
x
x2 1
x

x3
cos x
3
sin x
2x
x log x
2
x x tan
1
x
1
2
log x2
1
ť Determinant
Varje kvadratisk matris kan tilldelas ett bestämt tal, som kallas matrisens determinant. Detta tal kan beräknas med en formel som
innehåller alla matrisens element. Determinanter har mest teoretiskt värde och kommer till användning i diverse kvalitativa överläggningar, bland annat vid lösning av linjära ekvationssystem som vi kommer att studera i ett senare avsnitt.
Låt
beteckna en kvadratisk matris av ordning n
a11 a12
a21 a22
a1n
a2n
an1 an2
ann
n n
Till denna matris skall vi definiera ett enda tal som vi kallar determinanten till matrisen . Detta tal betecknar vi med
och skrivs ut med samma element som i fast med raka parenteser istället för runda.
det
a11 a12
a21 a22
a1n
a2n
an1 an2
ann
eller det
Benämningarna element, rad, kolonn, ordning och huvuddiagonal används på samma sätt för determinanter som för matriser. Hur
räknas då detta tal ut? Vi börjar från början
Definition.
(1) Om är av ordning 1 så är
(2) Om
är av ordning 2 så är
a11
a11 .
a11 a12
a21 a22
a11 a22
a12 a21 .
Determinater av högre ordning måste utvecklas enligt en visst mönster till en följd av determinanter av ordning 2. Vi beskriver här
hur denna utveckling av determinanter går till.
Antag att determinanten är av ordning n. Vi kan välja att utveckla den längs vilken rad eller kolonn som helst. Säg att vi väljer rad i.
För varje element i denna rad aij skapar vi en ny determinant genom att stryka rad i och kolonn j. Det som återstår är en ny determinant av ordning n
Cij
1
i j
ij
1 som kallas för underdeterminanten
till plats ij. Determinanten till
ij
med avseende på plats ij. Bilda komplementet (eng. cofactor)
är nu slutligen skalärprodukten mellan elementen i rad i och deras komplement.
Om nu n 1 2 måste alla dessa n 1 determinanter utvecklas enligt samma recept till n 1 st n 2 determinanter å igen å igen
ända till dess man kommit ner till (duktigt många) determinanter av ordning 2 som beräknas enligt ovan. Färdig
HH/ITE/BN
Talet
1
i j
Matriser och Mathematica
11
kan anta två värden 1 och 1 och uppträder i ett mönster som liknar ett "schackbräde"
Ŝ
Ŝ
Ŝ
ś ś ś
Ş
Detta hjälper oss att genomföra kalkylen som klarnar efter ett
2
4
1
Exempel: Bestäm
1
7
6
5
3 .
8
Lösningsförslag: Eftersom vi har en determinant av ordning 3 måste vi utveckla. Då väljer man en rad eller kolonn som innehåller
många nollor. Detta reducerar arbetsinsatsen dramatiskt. Men inga nollor i sikte här så vi väljer lite slumpartat att utveckla längs rad
2. Tecknen framför termerna i skalärprodukten mellan elementen och underdeterminanterna hämtar vi från "schackbrädet" ovan.
Redan efter en utveckling är vi framme vid 3 stycken determinanter av ordning 2 som vi kan "räkna" ut direkt. Vi får
Ʋ
Ŝ
Ŝ
Ŝ
ś ś ś
2
4
1
1
7
6
5
3
8
4
4
4
1 8
2
4
1
Ş
1 5
2 1 5
7 3
7 4 7 3
3
6 8
1 6 8
1 5
2 5
2 1
7
3
6 8
1 8
1 6
5 6
7 2 8 5 1
3 2 6
4 38 7 11 3 13 114
2
4
1
1
7
6
5
3
8
1 1
Detta sätt att beräkna en determinant är mycket kostsamt vid stora n. Om vi som illustration väljer att bestämma determinanten för
en, kan man tycka måttligt stor, determinant av ordning 25 med hjälp av en dator som klarar 1 miljon multiplikationer per sekund (1
Mflops) så får vi vänta en stund vid datorn; ca 1017 sekunder, vilket är ungefär universums ålder ! Vid all seriös numerisk analys
med dator är metodval mycket viktigt och kanske helt avgörande för att lyckas. I ett senare avsnitt ska vi se hur vi kan få svaret på
bråkdelar av en sekund.
I Mathematica används en funktion med självdokumenterande namnet Det som väljer olika algoritmer beroende på determinantens
storlek. Lägg märke till att den ska matas med en matris. Här senaste exemplet.
2
Det 4
1
1
7
6
5
3 
8
114
Räknelagar: Om
1
2
3
4
5
6
7
och
är matriser av ordning n och s ett tal gäller
Om alla element i en viss rad kolonn multipliceras med samma faktor, multipliceras
a11 sa12
a11 a12
determinanten med samma faktor, t.ex.
s
.
a21 sa22
a21 a22
s
sn
Determinantens värde ändras inte om man gör en linjärkombination av rader eller kolonner,
a11 a12
a11 sa21 a12 sa22
t.ex.
.
a21 a22
a21
a22
Om två rader eller kolonner byter plats, byter determinanten tecken.
Om två rader eller kolonner är lika, är determinanten noll.
12
Matriser och Mathematica
HH/ITE/BN
ť Invers matris
Om
är en kvadratisk matris och det existerar en (kvadratisk) matris
matrisen
kallas inversen till
sådan att
1
. För denna använder vi beteckningen
1
säger vi att
är inverterbar och
. Alltså
1
Man kan visa att om 1 existerar så är den unik, det vill säga det finns bara en enda. Ty om är en invers till och även gör
0.
anspråk på att vara det, har vi
, vilket visar att faktiskt
. Om är inverterbar så är
I någon mening kan man säga att invers matris fungerar som inverterade värdet gör för vanliga tal, vilket också skrivs "upphöjt till
1",
7
7
1
7
7
7 7
1
1. Därmed är det kanske inte helt fel att säga att vi har matrisdivision även om man aldrig får skriva
Med hjälp av inversen kan vi skriva ner den symboliska lösningen till matrisekvationer, t.ex.
söks. Lösningen får vi genom att förmultiplicera båda sidor med inversen till
1
1
1
1
, där
och
1
.
är givna och
.
I Mathematica beräknas inversen med hjälp av funktionen Inverse.
3 2
Inverse 4 6
2 7
17
13
16
71
2
71
11
71
19
71
16
71
17
71
10
71
71
71
1
5 
3
En matris av typ 2 2 ska man kunna invertera för hand,
1

a11 a12

a21 a22
1
1

a22
a21
a12

a11
1
a11 a22 a12 a21

a22
a21
a12

a11
En enkel kalkyl verifierar detta
1

a11 a22

a12 a21
a22
a21
a12
a
a
. 11 12 
a11
a21 a22
1 0

0 1
Exempel: Lös matrisekvationen
2
2
2 1
och
1 0
då
4 6 .
Lösningsförslag: Stuva om och lös ut . Notera speciellt efter andra ekvivalenspilen att det vid faktoriseringen är nödvändigt att
skjuta in en enhetsmatris av passande typ.
2
2
2
2
2
2
2
2
1
.
Nu låter vi Mathematica göra räknandet.
2 1
;
1 0
Inverse
4
;
6
.
2 IdentityMatrix 2
.
16
7
10
7
En ängslig koll
. .
2
True
En symmetrisk matris kallas ortogonal om
. För en sådan matris är alltså 1
. Ett nödvändigt och tillräckligt villkor
för att skall vara ortogonal är att :s kolonnvektorer är normerade och parvis ortogonala.
HH/ITE/BN
Matriser och Mathematica
13
Liksom determinanter förekommer inversen nästan aldrig i praktiken utan används mest vid diverse teoretiska överläggningar. Om
man mot förmodan behöver beräkna den finns effektiva algoritmer. Detta återkommer vi till i ett senare avsnitt.
Räknelagar: Om
vektorer gäller
1
1  1
1
2
3
4
5

1
1
och
är matriser av samma ordning,
och
ortogonala matriser av samma ordning och
6

1
Notera ordningen i högerledet
speciellt
1
ortogonal
8
9
1
passande
ortogonal
7
1
och
Enhetsmatrisen är ortogonal
10
ť Linjära ekvationssystem
System av linjära ekvationer utgör grunden för den linjära algebran på ungefär samma sätt som gränsvärden och derivator utgör
grunden för analysen. I princip alla modeller inom tillämpad matematik, allt från hållfasthetsberäkningar med finita elementmetoden
via ekonomisk resursteori till flödesproblem i en blodåder eller runt ett flygplan utmynnar alltid i ett stort linjärt ekvationssystem.
Att utveckla effektiva lösningsalgoritmer för sådana problem har därför sysselsatt matematiker i 70 år och kommer att bli än mer
aktualiserat i framtiden. Men vi börjar från början.
Med en linjär ekvation i två obekanta x och y menar vi en ekvation av typen
ax
by
c
där a, b och c är givna reella tal. Ordet linjär syftar på att de obekanta endast får förekomma på formen x, y. Sammansättningar och
funktionsuttryck är inte tillåtet. Som exempel har vi en linjär ekvation
6x
2y
10
medan
6x2
2y
10 och x
2
y
1 och x y
y
3
inte är linjära. Vi ser att t.ex. x 1, y 2 satisfierar ekvationen 6x 2 y 10 ovan och man säger att detta är en lösning. En annan är
x 0, y 5. Man kan hitta fler, i själva verket oändligt många eftersom för varje x har vi ett tillhörande y 5 3x. Detta beror på
att vi har fler obekanta än ekvationer.
Med ett linjärt ekvationssytem menar vi ett antal linjära ekvationer som gäller samtidigt. Även här kan man ha fler obekanta än
ekvationer, färre eller lika många. Man talar då om underbestämt, överbestämt respektive kvadratiskt linjärt ekvationssystem.
Alla fallen är mycket vanliga i diverse tillämpningar, speciellt det sistnämnda. Här ett exempel med tre ekvationer och tre obekanta.
2x
3x
4x
3y
4y
2y
9z
5z
2z
5
1
6
3
4
2
9
5
2
x
y
z
Det är brukligt att skriva linjära ekvationssystem på matrisform
2x
3x
4x
3y
4y
2y
9z
5z
2z
5
1
6
2
3
4
5
1
6
där kallas koefficientmatris, högerled och lösningsvektor. Om högerledet är nollvektorn säger vi att ekvationssystemet är
homogent. För stunden ska vi koncentrera oss på kvadratiska system och eftersom dessa dominerar branchen brukar de också lägga
beslag på namnet linjära ekvationssystem. Vill man prata om under- eller överbestämt linjärt ekvationssystem får man lov att nämna
detta explicit. Senare kommer vi att studera överbestämda system i samband med minsta kvadratmetoden och underbestämda vid
linjärprogrammering.
Det ser sunt ut med lika många ekvationer som obekanta, men vi kan faktiskt ha allt från entydig lösning, ingen lösning till oändligt
och . En geometrisk betraktelse av situationen blir speciellt åskådlig för ett
många lösningar! Vilket som inträffar beror på
linjärt ekvationssystem med två ekvationer och två obekanta x, y. Varje ekvation blir en rät linje i planet.
14
Matriser och Mathematica
x
x
y
y
3
x
x
2
y
y
y
2
HH/ITE/BN
3
2
x
x
y
y
y
y
5
5
5
4
4
4
3
3
2
2
1
1
1
1
2
x
2
3
2
1
0
unik lösning
3
3
1
2
x
2
0
ingen lösning
1
0
1
2
x
0
oändligt med lösningar
Vilket fall som inträffar beror på
Om
0 säger vi att ekvationssystemet är illa konditionerat. Genom att studera figurerna ovan förstår vi att små störningar av
elementen i eller ger en dramatisk ändring av lösningsvektorn . Sunda ekvationslösare, exempelvis Solve i Mathematica,
signalerar en varning när så är fallet. Vanligtvis beror detta på svagheter i modelleringen av det problem som man vill lösa.
Exempel: Betrakta de båda ekvationssystemen
x
x
Det första systemet har lösningen x, y
1
2
a
y
y
1, a
a
1
och
x
y
x
999
1000
1 och det andra x, y
a
y
1
1000
999a, 1000 a
1
. För a
1 har båda
systemen samma lösning, x, y
1, 0 . Men om vi nu ändrar a något lite till a 1.05, så ändras det första systemets lösning
obetydligt till x, y
1.025, 0.025 , medan det andra ändras våldsamt till x, y
48.95, 50 . Det första systemet är därför
väl konditionerat, det vill säga okänsligt för små störningar, medan det andra systemet är illa konditionerat.
Vi fortsätter med några exempel på modellering som leder till linjära ekvationssystem. Vi koncentrerar oss på detta för den praktiserande ingenjören viktigaste momentet och överlåter åt Mathematica att uttala sig om lösningsmängder. Dock ska vi i ett senare
avsnitt orientera om den vanligaste datoranpassade lösningsalgoritmen.
Exempel: AB Len&Fin tillverkar en kräm som enligt reklamen sägs motverka rynkor. Denna kräver tre olika råvaror. Inköpspriset
per gram råvara är 1 kr, 1.50 kr respektive 2 kr. Fraktkostnaderna per gram råvara är 2 kr, 1 kr respektive 1.50 kr. Till kund levereras
burkar med kräm som väger 50 gram och betingar ett pris av 80 kr i råvarukostnad och 70 kr i fraktkostnad. Hur många gram av de
olika råvarorna går det åt för att tillverka en burk?
Lösningsförslag: Antag att råvarornas storlekar är x, y respektive z. Informationen räcker för att möblera ett ekvationssystem med
tre ekvationer, en för vikten, en för råvarukostnad och slutligen en för fraktkostnaden.
x y z 50
1x 1.5 y 2z
2x 1 y 1.5z
80
70
Inte helt oväntat är det vår gamle arbetshäst Solve i Mathematica som gör jobbet. Fördelen är bland annat att den lägger sig väldigt
nära modelleringen. Inga riskabla omlastningar av ekvationerna till formen
vilket många datorprogram kräver, t.ex. Matlab.
Solve
x
x y z
50,
1 x 1.5 y 2 z
2 x 1 y 1.5 z
10., y
20., z
80,
70
20.
Vilket är svaret på den brännande frågan angående antal ingående gram av de olika råvarorna.
Exempel: Tre vätskor blandas i volymförhållande 1:2:4 varvid blandningen får densiteten 0.9. Om de blandas i volymförhållande
2:3:1 blir densiteten 0.8 och om de blandas 3:1:5 blir den 1.1. Ställ upp och lös det ekvationssystem som bestämmer densiteterna för
de i blandningen ingående tre vätskorna.
Lösningsförslag: Antag att densiteterna är Ρ1 , Ρ2 respektive Ρ3 . Informationen räcker för att meka ihop ett ekvationssystem med tre
ekvationer, en för varje soppa
Solve
1 Ρ1 2 Ρ2 4 Ρ3
2 Ρ1 3 Ρ2 1 Ρ3
1 2 4 0.9,
2 3 1 0.8,
HH/ITE/BN
Matriser och Mathematica
3 Ρ1
Ρ1
1 Ρ2
1.41429, Ρ2
5 Ρ3
0.3, Ρ3
3
1
15
5 1.1
1.07143
Exempel: För att bestämma temperaturfördelningen i en kvadratisk platta indelas denna i sexton mindre
kvadratiska plattor arrangerade i fyra rader och fyra kolonner. Temperaturen i plattorna längs randen är
kända. Se fig Bestäm temperaturen i de fyra inre plattorna om man antar att temperaturen i en sådan
platta är medelvärdet av temperaturen i de fyra kantgrannarna. Metoden kallas finita differensmetoden
FDM och differentialekvationen vi tillämpat den på heter Laplace differentialekvation och är
mycket vanlig i fysik.
38 46 20 19
35 T1 T2 15
10 T3 T4
5
8
1
5
0
C
Lösningsförslag: Vi får en ekvation för var och en av de fyra inre plattorna enligt meningen "Bestäm "
1
T1
4
SolveT1
1
46
T2
T3
35 , T2
20
4
1
T3
T1
T4
5
10 , T4
T2
4
T1
28, T2
18, T3
15
T4
T1 ,
4
1
13, T4
5
0
T3 
First
4
9
38 46 20 19
35 T1 T2 15
ListContourPlotReverse
 . T1 4 , InterpolationOrder
10 T3 T4 5
8 5 0 1
2, ContourLabels
True
Exempel: En tunn rektangulär glasskiva med massan m hålls i horisontellt läge
av tre personer enligt figur, där koordinatsystem och diverse mått är angivna.
a Sök ortsvektorerna A , B , C till personerna och G till skivans tyngdpunkt.
b Frilägg glasskivan med informationen att personerna endast lyfter i z–led.
c Bestäm med Solve de tre lyftkrafternas z–komponenter ur mekanikens
jämviktssamband
och
.
i
i
i
i
Lösningsförslag: Det är bara att följa receptet. Först a) Mekanik räknas alltid i 3D.
A
B
C
G
360, 0, 0 ;
720, 480 720, 0 ;
720 1680, 480, 0 ;
1
720 1680, 480 720, 0 ;
2
b) Personernas sökta lyftkrafter, enligt uppgift verkande endast i z–led. Slutligen tyngdkraften som verkar på glasskivan i negativ z–
led.
A
B
C
G
0,
0,
0,
0,
0, FAz ;
0, FBz ;
0, FCz ;
0, mg ;
c) Först jämviktssambandet med två vektorekvationer
i
och
i
i
i
.
16
Matriser och Mathematica
jämvikt

A
B
0
1200 FBz
600 mg
G,
C
A
A
B
B
C
C
G
360 FAz
mg
720 FBz
2400 FCz
Första raden i matrisen är kraftjämvikt i de tre koordinatriktningarna
0,
motsvarande för momentjämvikt Mx
2 3
My
0
G
0
480 FCz 1200 mg
HH/ITE/BN
0 och
Mz
Fx
FBz
FAz
FCz
0
0,
Fy
0 och
Fz

0
0. Andra raden är
0. De två vektorekvationerna expanderar alltså till
6 skalära ekvationer. Så blir det alltid i mekanik. Nu är det bara att lösa ut de sökta kraftkomponenterna.
Solve jämvikt, FAz , FBz , FCz
23 mg
FAz
57 mg
, FBz
79
158
55 mg
, FCz
2
Exempel: Lös matrisekvationen

158
2
då
2 1
och
1 0
4 6 .
Lösningsförslag: Detta är en repris på ett exempel under avsnittet "Invers matris". Jämför med tidigare lösningsmetod när vi denna
gång istället betraktar matrisekvationen som ett linjärt ekvationssystem och låter Solve direkt bestämma de obekanta elementen i
. Det enda vi behöver göra är att ansätta rätt. En titt i högerledet övertygar oss om att typ
typ , så smidigt och lätt
2 1
;
1 0
ekv
. .
5 x11 2 x21
2 x11 x21
4
;
6

x11
;
x21
2
2 x11
2 x21
. Solve ekv
4
6
First
16
7
10
7
I Mathematica finns väsentligen två funktioner som löser ekvationssystem. Vilken man väljer beror lite på hur problemet är formulerat. Om vi har koefficientmatris och högerled färdigmöblerade passar LinearSolve bra. Exempelvis exemplet med AB
Len&Fin
x y z 50
1x 1.5 y 2z
2x 1 y 1.5z
LinearSolve
1 1
1
1 1.5 2 ,
2 1 1.5
80
70
50
80 
70
10.
20.
20.
LinearSolve
1 1
1
1 1.5 2 , 50, 80, 70 
2 1 1.5
10., 20., 20.
Annars blir oftast valet den gamla välbekanta arbetshästen Solve. Den har fördelen att den ligger nära "modelleringen" och ger ett
självdokumenterat svar på regelform, vilket underlättar fortsatt användning. Dessutom är den smart nog att dra nytta av svagt
kopplade ekvationer, det vill säga många nollor i koefficientmatrisen, och sparar systemet i så kallat glest format. Detta sänker
lösningstiden dramatiskt om man har väldigt många obekanta. I många tillämpningar är det vardag att lösa ekvationssystem med
flera miljoner obekanta. Vi ser att även då det är bäddat för LinearSolve är kanske ändå Solve att föredra.
Solve
x
10., y
1 1
1
x
1 1.5 2 . y
z
2 1 1.5
20., z
50
80 
70
20.
När man matar Solve gäller att båda sidor om
tecknet ska ha identisk struktur eller att den ena sidan är en skalär.
HH/ITE/BN
Matriser och Mathematica
17
ť Rotation (vridning) av vektor kring en koordinataxel
Låt vara en tvådimensionell vektor och studera en rotation av denna kring origo. Positiv rotation
en vinkel Θ definierar vi som moturs rotation. I figuren återges situationen före och efter rotationen.
Eftersom vektorns längd bevaras får vi med summaformler för sinus och cosinus följande kalkyl
Θ , sin
Θ
vΘx , vΘy
Θ
Θ cos
cos cos Θ
cos Θ
Θ cos
cos Θ
sin Θ
sin Θ cos Θ
Θ
sin Θ , sin cos Θ cos sin Θ
sin Θ , Θ sin cos Θ
Θ sin
Θ cos
cos
Θ
sin Θ
cos
Θ
Θ
sin
sin Θ cos Θ
Θ
y
vΘ
v
sin
Θ
sin Θ
vx
vy
x
Nu generaliserar vi till det tredimensionella fallet och inser efter en stunds eftertanke att vi just hanterat positiv vridning kring zaxeln. Eftersom vz bevaras i detta fall, gäller det bara att expandera matrisen ovan till typen 3 3 och meka in en 1:a på rätt ställe. På
samma sätt kan analysen genomföras för rotation kring x- och y-axlarna. Vi lämnar detta som övning. Så, till slut får vi de tre
fundamentala rotationsmatriserna kring respektive koordinataxel. Notera att det inte är feltryck i y Θ
x
1
0
0 cos Θ
0 sin Θ
Θ
0
sin Θ ,
cos Θ
y
cos Θ
0
sin Θ
Θ
0 sin Θ
1
0
,
0 cos Θ
z
cos Θ
sin Θ
0
Θ
sin Θ
cos Θ
0
0
0
1
Exempel: Bestäm matrisen för den avbildning som erhålls om man först roterar 30 kring x-axeln och sedan 60 kring y-axeln.
Lösningsförslag: Notera att ordningen vid matrismultiplikation är viktigt som vanligt! Så den sammansatta avbildningen blir
, där
a
y 60
x 30
Cos 60
0
Sin 60
0 Sin 60
1
0
0 Cos 60
1
3
3
2
4
4
3
0
1
0
. 0 Cos 30
0 Sin 30
0
Sin 30
Cos 30
1
2
2
3
1
3
2
4
4
Om vi tar rotationerna i omvänd ordning får vi något helt annat...
1
0
0 Cos 30
0 Sin 30
1
2
0
0
Sin 30
Cos 30
.
Cos 60
0
Sin 60
0 Sin 60
1
0
0 Cos 60
3
2
3
3
4
2
1
4
3
1
3
4
2
4
ť Basbyte och koordinattransformation
Låt vara en vektor och låt , , och , , vara två olika baser i rummet. Vi söker
ett samband mellan :s komponenter vx , v y , vz respektive vx , v y , vz i de två
baserna. Baser i planet kommer som vanligt ut som ett specialfall genom att bara
betrakta de två första basvektorerna, vilket också är lättare att åskådliggöra. Eftersom
en vektor, som ju ofta har en fysikaliskt tolkning, måste vara densamma oberoende
av "vem" som mäter den med sin "linjal", får vi sambandet mellan :s komponenter i
de två baserna
vx
vy
För att komma vidare måste vi veta hur basvektorerna , ,
vz
och , ,
vx
vy
y
y
j
j
x
i
i
vz
ligger i förhållande till varandra. Antag att
x
18
Matriser och Mathematica
där tij är komponenterna för basvektorerna , ,
t11
t21
t31
t12
t22
t32
t13
t23
t33
HH/ITE/BN
med avseende på basvektorerna , , . Sätter vi in detta i uttrycket ovan för :s
representation i de två systemen och identifierar koefficienter framför , och
i höger och vänster led får vi
vx
t11 vx
t12 v y
t13 vz
vy
t21 vx
t22 v y
t23 vz
vz
t31 vx
t32 v y
t33 vz
Detta är det sökta sambanden mellan komponenterna vx , v y , vz och vx , v y , vz i de två baserna. Detta kan skrivas på matrisform
Matrisen
vx
t11 vx
t12 v y
t13 vz
vy
t21 vx
t22 v y
t23 vz
vx
vy
vz
t31 vx
t32 v y
t33 vz
vz
t11 t12 t13
t21 t22 t23
t31 t32 t33
vx
vy
vx
vy
vx
vy
vz
vz
vz
kallas transformationsmatrisen mellan de två systemen. Den talar alltså om hur basvektorerna ,
och
är riktade i
förhållande till basvektorerna , och . Observera att elementen i :s första kolonn utgör :s komponenter med avseende på basen ,
och . Analogt bildas andra och tredje kolonnen i med hjälp av komponenterna för respektive med avseende på basen ,
och . Om båda baserna är ortonormerade, alltså ON-system, så är kolonnerna i matrisen parvis ortogonala och normerade, det vill
säga är en ortogonal matris, så 1
. Därmed är det inte så "kostsamt" att "gå" mellan de olika representationerna. Även de tre
rotationsmatriserna x Θ , y Θ och z Θ i föregående avsnitt är ortogonala.
Vi ska nu använda detta för att bestämma ett samband mellan koordinaterna
x, y, z och x, y, z för en punkt P i koordinatsystemen O, , ,
y
O, , , . Om origo O för det senare systemet har koordinaterna xO , yO , zO 
med avseende på systemet O, , ,
x
xO , y
yO , z
j
y
så har vektorn OP komponenterna
zO  med avseende på basen , ,
j
O
och komponenterna
i
där
xO
y
yO
z
zO
x
O
kan vi så äntligen teckna sambandet mellan koordinaterna x, y, z och x, y, z
x
x
i
x, y, z med avseende på basen , , . Med hjälp av teorin ovan för basbyte
för en punkt P i koordinatsystemen O, , ,
.P
respektive
respektive O, , , 
x
y
z
x
y
z
xO
x
yO
y
z
zO
x
x
y
z
y
z
xO , yO , zO  avslutningsvis är ortsvektorn för O med avseende på basen , , .
Exempel: I ett tvådimensionellt fall är koordinaterna
y
för punkten P givna i systemet O, ,  som är vridet
30 moturs i förhållande till O, , . Sök koordinaterna
för motsvarande punkt P i O, , .
y
.P
j
i
O
i
Lösningsförslag: Vi har origo och transformationsmatrisen
3, 1 ;
3
1
2
2
1
3
2
2
Cos 30
Sin 30
Sin 30
Cos 30
x
30
j
O
1.3,0.9
3,1
x
 ; , det vill säga och är kolonner.
HH/ITE/BN
Matriser och Mathematica
19
så motsvarande punkt P i O, ,
P
. 1.3, 0.9
3.67583, 2.42942
En sista ängslig kontroll att vi kommer tillbaka till P i O, ,  med vetskapen att
.P
är ortogonal.

1.3, 0.9
ť Egenvärden och egenvektorer
y
Låt vara en given kvadratisk matris. Om det finns en vektor och ett reellt eller komplext
tal Λ så att
Λ
kallas Λ för ett egenvärde till och en till Λ hörande egenvektor. I dessa sammanhang
bortser vi från den ointressanta, så kallade triviala lösningen att är nollvektorn. Vidare ska
vi betrakta komplexa egenvärden lite styvmoderligt eftersom Λ då inte ligger i samma rum
som . Med ett språkbruk från vektorernas värld kan man alltså säga att en egenvektor till
Ax Λx
x
x
är en vektor som efter transformation
är parallell med sig själv, det vill säga
Λ . Egenvärdesanalys används
mycket flitigt inom de flesta tillämpningar. Här några enkla observationer av lineariteten som ofta kommer till användning
1
2
Om 1 och 2 är egenvektorer till för samma egenvärde Λ, så är också 1
2 egenvektor till
för samma egenvärde Λ, ty
Λ
Λ
Λ
1
2
1
2
1
2
1
2 .
Om är egenvektor till för egenvärde Λ och s ett reellt tal så är också s egenvektor till
för samma egenvärde Λ, ty s
s
sΛ
Λ s . Egenvektorer är alltså bara bestämda så när som på en syftlinje.
Vi får efter omstuvning ett homogent ekvationssystem
Λ
Λ
Ett nödvändigt och tillräckligt villkor för att detta system ska ha icke-triviala lösningar är att koefficientmatrisens determinant är
noll. Detta kallas för sekularekvationen och är en polynomekvation av grad n i Λ om är av ordning n.
Λ
0
Denna ekvation har alltid n rötter om dessa räknas med multiplicitet. Mängden av alla egenvärden till
Exempel: Bestäm egenvärden och motsvarande egenvektorer till
kallas
:s spektrum.
3 2
.
1 2
Lösningsförslag: Vi får sekularekvationen
3 2
1 2
Λ
1 0
0 1
Λ2
Egenvektor till Λ1
1. Ansätt
e1

3
0
x
 så får vi
y
Λ
2
1
5Λ
3
x
4
2
0
0
Λ
Λ1
3
1 och Λ2
1 x 2y 0
2 1 y 0
x
x
y
y
Λ 2
Λ
2 1
0
4
0
. Ekvationerna är parallella som sig bör. Annars
0
är det felräknat! Vi kan nu välja exempelvis x 1 och y
1. Egenvektorerna är bara bestämda med avseende på syftlinje! Man
brukar välja så man får "smidiga" tal eller så är man konsekvent och normerar egenvektorerna. Så första paret med egenvärde och
tillhörande egenvektor Λ1 1, e1 1, 1 .
Egenvektor till Λ2
4. Ansätt
e2

x
 så får vi
y
3
x
4 x 2y 0
2 4 y 0
andra paret med egenvärde och tillhörande egenvektor Λ2
4,
e2
x
x
2, 1
2y 0
. Ok! Välj exempelvis x
2y 0
2 och y
.
I Mathematica gör Eigensystem hela jobbet eller Eigenvalues och Eigenvectors för portionsvis leverans
Eigensystem
3 2

1 2
1, varav
20
Matriser och Mathematica

4
2, 1
HH/ITE/BN
1

1, 1
Egenvärdena till en reell symmetrisk matris är alla reella och egenvektorerna är parvis vinkelräta. Reella symmetriska matriser är
mycket vanliga i diverse tillämpningar.
1 2
.
2 1
Exempel: Bestäm egenvärden och motsvarande egenvektorer till
Lösningsförslag: Vi får sekularekvationen
1 2
2 1
1 0
0 1
Λ
Λ2
Egenvektor till Λ1
1. Ansätt xe1

1
0
2
2Λ
3
1
2x
x
 så får vi
y
Λ
2
1
0
0
Λ
Λ1
1
Λ 1
1 och Λ2
3
0
0
y
y
1 x 2y
1
1 y
x
x
Λ
2 2
0
0
. Ekvationerna är parallella som sig bör.
0
Annars är det felräknat! Vi kan nu välja exempelvis x 1 och y
1. Egenvektorerna är bara bestämda med avseende på syftlinje!
Man brukar välja så man får "smidiga" tal eller så är man konsekvent och normerar egenvektorerna. Så första paret med egenvärde
1, xe1 1, 1 .
och tillhörande egenvektor Λ1
Egenvektor till Λ2
3. Ansätt xe2

1 3 x 2y
2x 1 3 y
x
 så får vi
y
andra paret med egenvärde och tillhörande egenvektor Λ2
3, xe2
0
0
x
x
1, 1
y
y
0
0
. Ok! Välj exempelvis x
1 och y
1, varav
.
Vi låter Mathematica kontrollera riktigheten i vår kalkyl
Eigensystem

3
1, 1
1 2

2 1
1

1, 1
I många tillämpningar finner vi ofta system bestående av flera delar, vilka är ömsesidigt beroende av varandra genom regelbundet
utbyte av förnödenheter och tjänster. Det kan vara ett biologiskt system, växter och djur, som livnär sig på varandra eller på annat
sätt är beroende av varandra. Det kan vara ett ekonomiskt system med ömsesidigt utbyte av varor, tjänster och kapital. Ett sådant
system kan åtminstone på ett ytligt plan beskrivas genom att man anger utbytet av förnödenheter vid varje tidsperiod. Detta leder
naturligt till så kallade utbytesdiagram och dess övergångsmatris. Vi exemplifierar först med en äng, följt av ugglor och råttor i
skogen för att avsluta med Leslies djurhållningsmodell.
Exempel: På en äng som lämnats åt sitt öde finns bara gräs och buskar. Man har observerat att
a Under en växtsäsong övertar buskar 30 av den mark som var bevuxen med gräs.
b Under en växtsäsong dör 20 av buskarna och marken återgår till att vara gräsbevuxen.
Hur länge dröjer det innan halva ängen är bevuxen med buskar om det var 10 från början ?
Finns det något jämviktstillstånd? Hur ser i så fall fördelningen av gräs och buskar ut ?
Lösningsförslag: Låt xg vara andelen gräs och xb andelen buskar på ängen, det vill säga xg
ett år och
n 1
xb
1, så är
n
året därpå. De två observationerna ger oss då
xg 0.3xg 0.2xb
xg
xb
a
n 1
xg
b
xb 0.2xb 0.3xg
b
a
xb
n 1
0.7 0.2
0.3 0.8
xg
xb
n 1
n
n
n
Så tillståndet från ett år till nästa ges av överföringsmatrisen .
1
10
7 2
;
3 8
Nu är det bara att starta med 10% buskar år 0 och se hur det ser ut efter ett år efter två år och efter tre år.
0.9, 0.1 ,
. 0.9, 0.1 ,
. . 0.9, 0.1 ,
. . . 0.9, 0.1
xg , xb
tillståndet
HH/ITE/BN

Matriser och Mathematica
21
0.9 0.65 0.525 0.4625

0.1 0.35 0.475 0.5375
Kanske går det mot ett jämviktstillstånd? En bild över utvecklingen under de första 20 åren ger stöd åt vår förmodan.
år 20; p NestList . &, 0.9, 0.1 , år ;
ListPlot p, PlotRange
0, 1.01 , 0, 1 , Joined True, PlotStyle Red, AxesLabel
GridLines Automatic, GridLinesStyle Directive GrayLevel 0.8 ,
Epilog
Red, MapThread Text 1, 2, Background White &, Range 0, år , p
xg , x b ,
xb
1.0
0.8
0.6
17
16
15
14
13
20
19
18
12
11
10
987654
3
2
0.4
1
0.2
0
xg
0.0
0.0
0.2
0.4
0.6
0.8
1.0
Verkar hopa sig vid 0.4, 0.6 . Vi gör en ängslig test några år senare tiotusen år
MatrixPower
, 10 000 . 0.9, 0.1
0.4, 0.6
Ja, så här brukar modellering med dator gå till. Oftast räcker det och man kommer kanske inte längre teoretiskt. Men här kan vi dock
späka oss lite till
Typiskt är att vi har en följd av vektorer 0 , 1 , 2 ,
och, bortsett från "startvektorn" 0 , så gäller det att varje vektor kan fås ur
föregående genom multiplikation med en konstant matris enligt n 1
n för n 0, 1, 2, 3,
n
Vi vill undersöka vad som händer med följden n , när n
. Vi har redan insett att upprepad användning av formeln ger n
0,
så vårt problem skulle också kunna sägas vara att undersöka potenserna av matrisen . Vi ska se att vi kan skaffa oss en formel för
n ! Vår kvadratiska matris har naturligtvis egenvärden och tillhörande egenvektorer
Eigensystem
1
1
2
2
 , 1
1, 1
3
Låt oss skriva vår startvektor 0 som en linjärkombination av egenvektorerna
bestämmer tal c1 och c2 sådana att 0 c1 e1 c2 e2
9
1
Solve
,
10
3
c1
5
2
c1  , 1
3

10
0
2
1
5 1
1, 1 
c2
2
3
5
1n
12

e1
e1
1

2
e2 
n:
3
5
e1
1 1

2 2
e2 
5 1
1

2

e1
3


e1
1

2
1 n
2
e2
1 2
2
3
5

1

2
3
3
5
e2 
e1
12
1 n 1
e1
2
e2
e2

e1
3
5
1
1 1

2 2

3
5
e2 
1

2
för alla n
2
 , 1
5 3

1
1 2
2
0, 1, 2,
n 1
1, 1
2
0.9 0.65 0.525 0.4625 0.4

0.1 0.35 0.475 0.5375 0.6
.n
0, 1, 2, 3, 10 000
1 1

2 2

e1
Vi testar på några tillstånd vi känner sedan tidigare
3
det vill säga vi gör ett basbyte och

1
3
5
e2 ,
2
3
5
n
och
1
, c2
Nu är det lätt att meka ihop en formel för
3
e1
N
e2
12
e1

e2 
3
5
13
1

2
e1
1 2
2

e2
1

2
1 3
2
e2
22
Matriser och Mathematica
Om matrisen
n
n
har egenvektorerna
c1 Λn1 e1 c2 Λn2
0
och
e1
e2
HH/ITE/BN
med egenvärden Λ1 respektive Λ2 och
0
c1
e1
c2
e2 ,
så gäller allmänt att
e2
Så här ser det ut om är en 2 2-matris. Om är en 3 3-matris har man 3 egenvektorer, 3 egenvärden och en summa av tre termer
och så vidare. Sedan finns det tyvärr "undantagsmatriser", för vilka egenvektorerna inte räcker till för att varje 0 skall kunna skrivas
som en linjärkombination av egenvektorerna. Sagan är sann även då egenvärdena är komplexa.
Nu drar vi oss till minnet att
Λn 0, när n
, om Λ
1
. Av detta kan vi inse följande
n
Λ
, när n
, om Λ
1
Om alla egenvärden har absolutbelopp 1, så n 0, när n
, varmed menas att alla komponenter i vektorerna n 0. Detta
är en slags stabilitet.
Om något Λ har absolutbelopp 1, har vi termer i summan som
(utom i undantagsfallet då motsvarande koefficient c 0).
Detta är instabilitet.
. Stabilitet.
Om ett av egenvärdena, säg Λ1 1 och övriga är till absolutbeloppet 1, så får vi att n c1 e1 , när n
I vårt fall gäller tydligen det sistnämnda, så
0?
olika starttillstånd
n
c1
e1
3 2
 ,
5 3
1
2
3
 5 , 5  vilket stämmer bra! Kvarstår bara en stilla oro angående
2
3
Vi kan starta i vilket tillstånd som helst, alla går mot  5 , 5 ! Det finns gott om tillämpningar där man vill att
överföringsmatrisen ska ha egenvärden enligt det sista alternativet ovan, exempelvis sluten ekonomi och stabil djurpopulation. Det
vi har studerat är ett så kallat diskret dynamiskt system, (vi återkommer till kontinuerliga i samband med differentialekvationer).
kallas bana eller trajektoria och en stabil punkt för attraktor. Figuren ovan är exempel på en
Följden av vektorer 0 , 1 , 2 ,
trajektoria med en attraktor.
Exempel: I storskogen är populationen av ugglor U och råttor R , räknat i tusental,
kopplade från månad k till nästföljande månad k 1 av ekvationerna
Uk 1
0.5Uk 0.4Rk
Rk 1
pUk 1.1Rk
Förklara utförligt innebörden av de fyra ingående koefficienterna Plantera nu lite U
och R i skogen och undersök utvecklingen under några månader för lite olika värden
på p. Förklara Bestäm sedan p så att det blir en stabil uppsättning djur över tiden.
Om det då finns 20 ugglor i skogen, hur många är råttorna då ?
Lösningsförslag: Vi börjar med de fyra koefficienterna. 0.5Uk innebär att vid avsaknad av råttor i storskogen så kommer hälften av
ugglorna att dö under månaden. På motsvarande sätt innebär 1.1Rk att vid avsaknad av ugglor så ökar antalet råtter med 10% under
en månad. Om nu båda arterna är närvarande så ökar ugglorna sin population med 0.4Rk medan pUk innebär att råttornas antal
minskar beroende på de jagande ugglorna, i snitt blir alltså 1000 p råttor uppätna av en uggla varje månad. Sedan har vi parametern
p. Av föregående exempel har det framgått att det vore en idé att titta på egenvärdena till överföringsmatrisen för olika värden på p.
0.5 0.4
;Λ
p 1.1
0.2 4.
2.25
Eigenvalues
10. p , 0.2 
2.25
10. p
4.
Eftersom vi inte här befattar oss med komplexa egenvärden så måste 0 p 0.225. Undre gränsen är naturlig eftersom ugglornas
närvaro knappast bidrar till att öka antalet råttor. Studera nu de två egenvärdena Λi som funktion av p i en graf
Plot Λ, p, 0, 0.225 , PlotStyle
Blue, Red , AxesLabel
"p", "Λ1 ,Λ2 "
Λ1 ,Λ2
1.1
1.0
0.9
0.8
0.7
0.6
p
0.05
Vi ser att 0
Λ1
0.10
0.15
0.20
1 och för något kritiskt värde pc har vi
Λ2 1, för p pc
. Så följande fall utkristalliserar sig
0 Λ2 1, för p pc
HH/ITE/BN
Matriser och Mathematica
23
För p pc brakar antal ugglor/råttor mot oändligheten (instabilitet). Modellen gäller inte så länge eftersom det varken finns mat
eller utrymme här på jorden till "så" många.
För p pc dör såväl ugglor som råttor ut. Med tanke på den biologiska mångfalden så är detta en något tråkig form av stabilitet.
För p pc har vi Λ2 1 och 0 Λ1 1 och vi går mot ett stabilt jämviktstillstånd i storskogen, oberoende av hur många djur vi
planterar ut i "begynnelsen". Detta kritiska värde är
pc
Solve Λ 2
p
1
First
0.125
med tillhörande egenvärden
Λ . pc
0.6, 1.
Kontrollera ovanstående förutsägelser genom att "simulera systemet", det vill säga räkna fram antalen ugglor/råttor ett antal tidsenheter framåt för olika värden på p. Skriv ut var tolfte månad. Det går bra (nästan bättre) att räkna för godtyckliga startvärden U och
R.
Vi börjar med p
pc och ser att både ugglor och råttor ökar till antalet med tiden.
f x
Nest .
.p
NestList f, U, R , 5
0.05 &, x, 12 ;
Simplify
U
1.60134 R 0.141263 U 2.26074 R
3.39401 R 0.30058 U 4.79043 R
7.19165 R 0.636908 U 10.1506 R
15.2386 R 1.34956 U 21.5083 R
32.2894 R 2.85962 U 45.5745 R
Sedan p
R
0.200167 U
0.424251 U
0.898956 U
1.90482 U
4.03618 U
pc och ser att både ugglor och råttor dör ut.
f x
Nest .
.p
NestList f, U, R , 5
0.2 &, x, 12 ;
Simplify
U
R
0.537176 R 0.254747 U
0.551018 R 0.268588 U
0.15915 R 0.0793833 U
0.159341 R 0.0795749 U
0.0450515 R 0.0225231 U
0.0450541 R 0.0225257 U
0.0127253 R 0.00636261 U 0.0127253 R 0.00636265 U
0.00359402 R 0.00179701 U 0.00359402 R 0.00179701 U
Sist det intressant stabila p pc som säger oss att oberoende vad vi startar med U, R kommer vi så småningom att få ett stabilt
läge R 0.25U, 1.25R 0.3125U .
f x
Nest .
. pc &, x, 12 ;
NestList f, U, R , 4
Simplify
U
R
0.997823 R 0.247279 U 1.24946 R 0.31182 U
0.999995 R 0.249994 U 1.25 R 0.312499 U
1. R 0.25 U
1.25 R 0.3125 U
1. R 0.25 U
1.25 R 0.3125 U
Avslutningsvis har vi frågan om det jämviktstillstånd som motsvarar 20 ugglor.
Eigensystem
Λ,

1.
0.624695, 0.780869
. pc
0.6

0.970143, 0.242536
U och R ska alltså förhålla sig som
1
1, 1
1., 1.25
Så 20 ugglor motsvarar 20 1.25 1000
25 000 råttor i storskogen.
24
Matriser och Mathematica
HH/ITE/BN
Exempel: Leslies populationsmodell P. H. Leslie 1945 är ett försök att beskriva utvecklingen hos
en population med hänsyn tagen till åldersstrukturen. Modellen räknar vanligen endast med hondjuren,
eftersom de är avgörande för reproduktionen. Man tänker sig således att populationen handjur hela tiden
kan hållas på en tillräckligt hög nivå för att modellens antaganden skall vara rimliga. Vi beskriver modellen
med tre olika åldersgrupper. Generellt tillåter modellen ett godtyckligt antal åldersgrupper. Låt hondjuren
delas in i tre åldersgrupper G0 de nyfödda , G1 och G2 . Vi skall sätta upp en modell för antalet djur i dessa
tre olika grupper vid tidpunkterna 0, 1, 2, 3, . Alla djur som vid en viss tidpunkt befinner sig i den sista
åldersgruppen G2 kommer vid nästa tidpunkt av vara döda. En viss andel pi av de hondjur som befinner
sig i gruppen Gi kommer att överleva till nästa tidpunkt. Vidare kommer varje hondjur i åldersgruppen Gi
i genomsnitt föda fi hondjur. Om det nu vid tidpunkten n finns xn 0 , xn 1 , xn 2 hondjur i de tre ålders–
grupperna får vi följande samband mellan kullarna vid tidpunkt n och nästföljande tidpunkt n 1.
xn
xn
xn
1
1
1
0
1
2
f0 xn 0
p0 xn 0
p1 xn 1
Om vi nu samlar åldersgrupperna i n xn 0 , xn 1 , xn 2
trisen som alltid är kvadratisk. Vi får direkt i vårt fall
f0 f1 f2
p0 0 0
0 p1 0
f1 xn 1
f2 xn 2
kan sambandet skrivas
n 1
och inser att den vid fler åldersgrupper generaliseras till
n,
där
är den så kallade Lesliema-
f0 f1 f2
p0 0 0
0 p1 0
0 0 p2
f3
0
0
0
Vi provkör med lite siffervärden på parametrarna
p0
p1
f0
f1
f2
0.5 50 av de nyfödda överlever
0.8 80 av ettåringarna överlever
0 inga hondjur i G0 reproducerar sig
1.6 1.6 honor föds i genomsnitt av varje ettåring
0.5 0.5 honor föds i genomsnitt av varje tvååring
0 1.6 0.5
0.5 0 0
0 0.8 0
Vi antar att åldersstrukturen vid tidpunkten 0 är känd. På grund av en miljökatastrof finns inga djur i grupperna G0 och G2 . Antalet
djur i G1 är 100. Det betyder att 0 0, 100, 0 . Följande tabell visar hur åldersstrukturen förändras under de första tjugo åren.
0 1.6 0.5
; NestList Round
0.5 0
0
0 0.8 0
.
&, 0, 100, 0 , 20
0 160 40 128 64 110 77 101 83 95 87 94 89 93 89 91 89 91 89 91 89
100 0 80 20 64 32 55 38 50 42 48 44 47 44 46 44 46 44 46 44 46
0 80 0 64 16 51 26 44 30 40 34 38 35 38 35 37 35 37 35 37 35
Frågan är nu om åldersfördelningen konvergerar, det vill säga med tiden närmar sig en viss stabil jämviktsfördelning. Ser så ut efter
lite svajig inledning. Då måste det gälla att
och fördelningen ges av den egenvektor som motsvarar egenvärdet ett.
Λ, e

Eigensystem
1.
0.842152, 0.421076, 0.336861
0.723607
0.696556, 0.481308, 0.532121
Så den sökta stabila jämviktsfördelningen blir
e 1
Total e 1
0.526316, 0.263158, 0.210526
0.276393

0.177642, 0.321358, 0.930146
HH/ITE/BN
Matriser och Mathematica
25
Exempel: I mekanik används egenvektorerna för att peka ut ett naturligt
koordinatsystem för en kropp. Vi känner igen begreppet masströghets
xx
yx
m
yy
y2 m
m
r m. Detta kan utvidgas
m
yx m
m
xy m
x2 m
y
till att man bestämmer den så kallade masströghetsmatrisen med analog
m
m
2
moment för en kropp runt en given axel I
uppbyggnad i 3D. Integralerna
xy
y
xy m kallas deviationsmoment. Samma
j
j
struktur och egenvärdesanalys finner man exempelvis i hållfasthetslära,
statistik och i sökmotorn Google som använder egenvektorer för att ta
x
i
genvägar när den söker genom enorma datamängder. Även egenvärdena
kan ges en fysikalisk tolkning. Men först modellen. En smal stång med
längden L och massan m är placerad i vårt xy koordinatsystem enligt
vidstående figur. Sök masströghetsmatrisen .
i
30
x
Lösningsförslag: Börja med att bestämma stångens "funktion"
y
Tan 30
x;
Sedan tröghetsmatrisen. Vi har m
L
Cos30 
2

L
Cos30 
2
L2 m
xy
x2
Ρ
y '2
1
m
1
D y, x
x i färskt minne
2
x
L
L2 m
48
16
L2 m
L2 m
16
y2
xy
Ρ s
3
16
3
Eigensystem
L2 m
12
1

0
, 1  3 , 1
3
Vi ser att egenvektorerna precis pekar ut det naturliga x y-systemet för kroppen! Dessa kallas för kroppens huvudtröghetsriktningar.
Vidare är egenvärdena lika med masströghetsmomenten (huvudtröghetsmomenten) runt respektive riktning, det vill säga de berättar
om kroppens utsträckning i respektive riktning.
Lite mer om egenvektorer och egenvärden
1
2
Om
Om
har en egenvektor med egenvärdet Λ så har n en egenvektor med egenvärdet Λn .
inverterbar
alla egenvärden till är skilda från noll.
3
Om
inverterbar och egenvektor till
4
Om alla element i är positiva har minst ett positivt egenvärde och alla komponenter i motsvarande
egenvektor är ickenegativa. Perrons sats .
5
Låt Ci
till
6
z : z
aii
i j
med egenvärdet Λ så är också egenvektor till
1
1
med egenvärdet Λ .
aij  vara cirklar i det komplexa talplanet. Då ligger samtliga egenvärden
i unionen av Ci , det vill säga i
Ci . Gerschgorins sats .
satisfierar sin sekularekvation. Cayley–Hamiltons sats . Detta märkliga resultat betyder alltså om
sekularekvationen är Λn cn 1 Λn 1
c1 Λ c0 0 så gäller att n cn 1 n 1
c1
c0
2
Exempel: Undersök om Gerschgorin fungerar på matrisen 3
5
3 1
4 5 .
8 7
z: z 2
5 8
Lösningsförslag: Gerschgorins tre cirklar får vi med receptet (5) ovan. C1
C2 z : z
4
3 5
z: z
4
8 och till sist C3 z : z 7
dena med sekularekvationen
2
3
5
3 1
4 5
8 7
Λ
1 0 0
Λ 0 1 0
0 0 1
1 Λ2
4Λ
2
80
Λ
3
5
0
0
Λ1
.
21
3
1
4 Λ
5
8
7 Λ
0
21 , Λ2
21
Λ3
5Λ2
3 1
z: z 7
76Λ
21 och Λ3
80
1
z: z 2
4,
13 . Sedan egenvär-
0
26
Matriser och Mathematica
HH/ITE/BN
En liten bild över situationen som bekräftar Gerschgorin kan inte skada. Vi låter Mathematica göra grovjobbet
2
3
5
2 1
3 1
4 5 ;Λ
8 7
21 , 2 1
Eigenvalues
21 , 1
Graphics Table c
i,
Blue, Green, Orange
Red, PointSize 0.025
i, 3
, AspectRatio
i , 0 ; r Plus
i , Circle c, r ,
, Point Λ i , 0
Automatic, Axes
Abs
i
Abs
i, i ;
Text Ci , c r 0.72 1, 1 ,
1, 0 ,
, Text Subscript "Λ", i , Λ i , 1.3
True, AxesLabel
x, y
,
y
10
C3
C2
5
Λ2
10
Λ3
5
C1
Λ1
5
10
15
20
1,
, n till matrisen
x
5
10
Vi ser att egenvärdena Λ j , j
3
i 1Ci .
ligger i unionen
ť Gauss eliminationsmetod för lösning av linjära ekvationssystem
Eftersom lösning av ekvationssystem är vardagsmat för alla världens datorer är det viktigt att ha stabila och effektiva algoritmer, inte
minst då antalet obekanta kanske uppgår till flera millioner. Man skiljer på iterativa och direkta ekvationslösare. Exempel på en
iterativ metod hittar vi i ett senare avsnitt.
Här ska vi studera den vanligaste direkta metoden som också varit förhärskande under många decennier. Men innan vi går vidare ska
vi betrakta ett litet ekvationssystem med två ekvationer och två obekanta
2x
4x
6y
7y
3
8
1
2
och göra följande mer eller mindre triviala observationer. Lösningen till ett ekvationssystem ändras inte om vi
1. Byter plats på två rader.
2. Multiplicerar en rad med ett tal skilt från noll.
3. Adderar en rad till en annan.
Detta kallas för elementära radoperationer. Ekvationssystem som kan omformas i varandra med elementära radoperationer kallas
radekvivalenta och har alltså samma lösning.
Egentligen borde man kanske prata om ekvationer i stället för rader men denna nomenklatur hör samman med att vi lite längre fram
kommer att genomföra arbetet i en matris. Punkt 1 är självklar och punkt 2 känner vi igen som "man ska göra samma sak på båda
sidor", det vill säga i exempelvis ekvation (2): 123 4x 7 y 123 8. Punkt 3 vilar på samma argument, så ekvation (2) igen:
4x 7 y 2 3 8 2 3. Men enligt ekvation (1) är 2x 6 y 3, vilket vi nu utnyttjar för att byta ut 3:an i vänsterledet i senaste
versionen av ekvation (2): 4x 7 y 2 2x 6 y 8 2 3
5 y 2. Genom att på ett sådant här smart sätt utnyttja väsentligen
punkterna 2 och 3 kan man reducera ett ekvationssystem till en enkel form som är lätt att lösa. Denna systematiska metod kallas
Gauss eliminationsmetod och man kan på grund dess enkelhet lätt dressera en dator att härma den. Häng mé!
2x
4x
Om vi multiplicerar ekvation (1) med
4
2
6y
7y
3
8
1
2
och subtraherar den från ekvation (2) får vi
HH/ITE/BN
Matriser och Mathematica
2x
4x
6y
3
7y
4
2
27
1
2x
6y
4
2
8
3
2
vilket hyfsas till
2x
6y
5y
3
2
1
2
Det är nu uppenbart att det "smarta" utnyttjandet av punkt 2 och 3 bestod i valet av just
4
2
så att x elimineras från ekvation (2) efter
multiplikation och subtraktion av ekvation (1). Strategin är att med elementära radoperationer arbeta sig systematiskt igenom
kolonnerna i med syfte att göra den högertriangulär, det vill säga så att alla element under huvuddiagonalen blir noll. Detta är
själva eliminationssteget (eng. elimination). Elementen på huvuddiagonalen i som används flitigt i nämnaren på den lämpliga
faktorn kallas pivotelement (eng. pivot element). De obekanta kan nu enkelt bestämmas genom successiv bestämning "underifrån",
så kallad bakåtsubstitution (eng. back substitution).
2 :
5y
2
2
5
y
1 : 2x
6
2

5
3
x
27
10
Gauss eliminationsmetod är därmed ett systematiskt sätt att bestämma lösningen till ekvationssystemet
utan att använda
sig av inversen till . Det brukar vara praktskt att utföra räkningarna i totalmatrisen eller utökade matrisen (eng. augmented
matrix), bestående av med en extra kolonn , vilket skrivs som kolonnkatenering
. Efter eliminationssteget säger man att
totalmatrisen har trappform (eng. row-echelon form). Som "biprodukt" efter eliminationssteget har vi även determinanten till
som produkten av talen i huvuddiagonalen.
Man kan visa att Gauss eliminationsmetod är den som kräver minst antal operationer (multiplikationer och additioner) för att lösa ett
givet ekvationssystem, lösningstiden är proportionell mot n3 , det vill säga T kn3 där n är antalet obekanta och k en konstant som
beror på datorn. Därför är Gauss en typisk "arbetshäst" i numerisk analys. Vi illustrerar metoden med ytterligare ett exempel, denna
gång med tre ekvationer och tre obekanta.
Exempel: Sök lösningen
till ekvationssystemet
2 4
4 6
6 8
, där
1
5 ,
4
2
3 och
6
x1
x2 . Använd Gauss eliminationsx3
metod.
Lösningsförslag: Först tittar vi i facit
2 4
Solve 4 6
6 8
1
x1
5 . x2
x3
4
13
x1
7
1
, x2
2
2
3 
6
2
, x3

7
2 4
Sedan raskt över till Gauss med totalmatrisen 4 6
6 8
1 2
5 3 . Använd systematiskt elementära radoperationer för att göra
4 6
högertriangulär. Med k menas talen i rad k i totalmatrisen och med
Kolonn 1: 2
2 4
4 6
6 8
4
2
2
1 2
5 3
4 6
2
0
0
4
2
4
1
7
7
4
2
6
2
4
6
4
2
Kolonn 2: 3
6
2
1 , sedan 3
2
1
0
1 :
4
2 6
2 8
4
2
6
2
1
4 5
4 4
4
2
6
2
2
1
3
1
6
4
2
6
2
2
0
0
2
2
4
2
4
1
7
7
2
1
0
2
0
0
4
2
0
1
7
7
7
28. Undrar ändå
2 :
2
0
0
4
2
4
4
2
1
7
2
7
4
2
2
1
7 0
Vi är nu färdiga med eliminationssteget och noterar att
2 4
Det 4 6
6 8
menas att systemen är radekvivalenta.
1
5 
4
4
2
1
2
2
2
1
2
28
Matriser och Mathematica
HH/ITE/BN
28
Slutligen bestämmer vi
Ekv 3 :
7x3
2
Ekv 2 :
2x2 7
genom bakåtsubstitution
2

7
1

2
Ekv 1 : 2x1 4
2
7
x3
1
1
x2
2

7
2
13
7
1
2
2
7
1
2
13
7
x1
Om vi skapar totalmatrisen
2 4
4 6
6 8
M
1 2
5 3 ;
4 6
kan vi sedan betrakta skådespelet som följande programsnutt i Mathematica spelar upp.
n Length M ;
Print "
Elimination
"
DoPrint "j ", j, ": ", M ;
M i, j
;M i
Doq
qM j ;
M j, j
a10 i
j
a10 j
j
Print"j ", j, ", i ", i, ", ",
, i, j
, " ", q, "
1, n ;
, j, 1, n 
Print "
Bakåtsubstitution
Table 0, n 1 ;
Do
i
M i, n 1
M i .
Print xi , "
",
i ;
, i, n, 1, 1
Elimination
2 4
j 1: 4 6
6 8
j 1, i 2,
a21
1 2
5 3
4 6
2
2
0
6
4
2
8
1
7
4
2
1
6
3
2
0
0
4
2
4
1
7
7
2
1
0
4
2
0
1
7
7
2
1
2
a11
j 1, i 3,
a31
a11
2
j 2: 0
0
j 2, i 3,
4
2
4
a32
1
7
7
2
a22
2
j 3: 0
0
4
2
0
1
7
7
2
1
0
2
0
0
2
1
2
Bakåtsubstitution
2
x3
7
1
x2
2
13
x1
7
"
M i, i ;
", M;
HH/ITE/BN
Matriser och Mathematica
29
I allmänhet erhåller man det bästa numeriska resultatet genom att välja så stort pivotelement som möjligt och för att åstadkomma
detta är det ibland nödvändigt att sortera om raderna enligt elementära radoperationer punkt 1. Detta kallas pivotering. En god
implementering av Gauss eliminationsmetod, exempelvis Solve i Mathematica, tar hänsyn till detta.
ť Jacobis metod för beräkning av invers matris
1
I praktiken används inversen nästan aldrig explicit. Ofta dyker den upp i sällskap med multiplikation av vektor,
. Då väljer
man istället att med hjälp av Gauss eliminationsmetod lösa ekvationssystemet
. Anledningen till detta är att det krävs
väsentligt mindre arbete i en dator.
Om vi trots allt är intresserade av inversen, så kan vi använda en till Gauss eliminationsmetod snarlik metod, nämligen Jacobis
bildar vi totalmatrisen
. Med samma elementära radoperationer som vid Gauss omformar vi sedan
metod. Istället för
1
denna till
1
. Att detta blir rätt inses av följande kalkyl där
är den sökta inversen till
;
1
1
.
3 2
4 6
2 7
Exempel: Bestäm inversen till
1
5 . Använd Jacobis metod.
3
Lösningsförslag: Först tittar vi i facit
3 2
Inverse 4 6
2 7
17
13
16
71
2
71
11
71
19
71
16
71
17
71
10
71
71
71
1
5 
3
3 2
Sedan raskt över till Jacobi. Vi får totalmatrisen 4 6
2 7
1 1 0 0
5 0 1 0 . Använd nu elementära radoperationer för att systematiskt
3 0 0 1
arbeta igenom de tre första kolonnerna så att dessa bildar en enhetsmatris. Med k menas talen i rad k i totalmatrisen och med
menas att systemen är radekvivalenta.
Kolonn 1: Dividera 1 med 3, sedan 2
3 2
4 6
2 7
1 1 0 0
5 0 1 0
3 0 0 1
1
1
0 1
0
17
3
1
3
19
10
11
3
1
3
2
5
2
3
1 0
0 1
0 0
1
3
5
2
5
16
71
5
3
10
,
3
1 0
3
10
0 1
0
0
0 1
1
5
3
10
17
71
1
3
0 0
1
0 1 0
0 0 1
0
sedan 1
0 0
71
,
10
17
3
8
5
19
10
11
3
2
3
3
5
2
5
2
3
1
5
3
10
0
8

5
0
1 0 0
0
0 1
17
71
2
5
16
71
10
71
0 0 1
Vi är nu färdiga och kan direkt avläsa
Om vi skapar totalmatrisen
1
2 1 :
1
3
19
3
1
3
4
3
0 0
1
1 0
0
3
0
0 1
0
17
3
2 , slutligen 3

19
10
2
3
10
3
2 7
sedan 1
Kolonn 3: Dividera 3 med
8
5
19
10
1
3
4 6
2 7
Kolonn 2: Dividera 2 med
2
3
2
3
4 1 , slutligen 3
17
71
2
71
16
71
0
1 0
0
0 1
1
0 0
3 , slutligen 2
13
71
3
10
17
71
16
71
13
71
11
71
17
71
3
5
2
5
8
5
19
10
0 1 0
10
71
0 0 1
16
71
19
71
10
71
.
1
3
19
3
11
3
2 :
8
5
19
10
71
10
1 0 0
0
2
3
10
3
17
3
1
5
3
10
17
10
0
0
1
3 :
17
71
2
71
16
71
13
71
11
71
17
71
16
71
19
71
10
71
1
3
4
3
2
3
0 0
1 0
0 1
30
Matriser och Mathematica
3 2
4 6
2 7
M
HH/ITE/BN
1 1 0 0
5 0 1 0 ;
3 0 0 1
kan vi sedan betrakta skådespelet som följande programsnutt i Mathematica spelar upp.
n
Length M ;
DoPrint "j ", j, ": ", M ;
M j
p
;
M j, j ; M j
p
Print "j ", j, ", ", a10 j j , " ", p, "
", M ;
Do If i j, q M i, j ;
M i
qM j ;
Print "j ", j, ", i ", i, ", ", a10 i j , " ", q, "
, i, 1, n
, j, 1, n 
3 2
j 1: 4 6
2 7
1 1 0 0
5 0 1 0
3 0 0 1
1
j 1, a11 3
2
1
1
3
3
3
4 6
2 7
5
3
1
j 1, i 2, a21 4
0
2
1
0
j 1, i 3, a31 2
0
1
j 2: 0
0
2
1
3
10
3
19
1
3
4
3
17
3
11
3
2
3
3
3
0
0
2
1
3
10
3
19
3
3
7
3
1
3
19
3
4
3
17
3
11
3
2
3
3
3
1
3
3
19
3
2
17
10
11
5
2
3
3
3
0 1
8
0
1
10
17
5
2
3
3
3
3
1
5
19
5
2
5
3
10
71
5
8
10
17
10
5
10
0 0
0
10
10
11
8
0 1
0 1
0
3
5
3
3
10
0
1
1
1 0
1
5
2
0
71
j 3, a33
1 0
3
0
0 0
0 0
0 0
5
19
1
17
j 3: 0 1
1
0 1
3
8
0 1
3
10
1 0
1 0
1 0
3
0
1 0
0
j 2, i 3, a32
0 0
3
4
2
2
j 2, i 1, a12
1
0 0
1
3
0 1 0
0 0 1
2
1
10
j 2, a22
0 0
0
0
0
1
3
1
5
19
5
2
5
3
10
71
5
8
10
17
10
5
10
0
0
1
0
0
1
8
3
1
5
19
5
2
5
3
10
5
16
10
17
71
71
1
0
0
10
71
", M
;
HH/ITE/BN
Matriser och Mathematica
1 0
8
j 3, i 1, a13
0 1
5
0 0
0
19
13
71
2
71
3
5
16
10
17
71
71
10
1
1 0 0
19
j 3, i 2, a23
17
0 1 0
10
0 0 1
31
16
71
0
10
71
17
13
16
71
2
71
11
71
19
71
16
71
17
71
10
71
71
71
ť En iterativ metod för lösning av linjära ekvationssystem
Fördelen med Gauss eliminationsmetod är att den i ett ändligt antal operationer ger rätt resultat på given tid T kn3 . Emellertid
hamnar man i stora svårigheter då problemstorleken ökar och inte länge får plats i datorns minne utan måste delas upp i delmatriser (partioneras) och lagras på disk. Detta ökar lösningstiden drastiskt till mycket längre än formeln ovan antyder. Vidare kan man
ju fråga sig om det är rimligt att bestämma en lösning med full maskinnoggrannhet, kanske 15-20 digitala siffror, om man betänker
att vid matematisk modellering gör man ju en modell av verkligheten, kanske full av diverse idealiseringar? Då kan det vara dax att
fundera på en iterativ lösare.
Iterativa metoder baserar sig väsentligen på att man minimerar kvadratsumman
. Detta kan göras mycket snålt
eftersom algoritmen vanligtvis endast kräver vektorn
explicit. Vid mycket stora problem är detta en dramatisk reduktion av
minnesbehovet i datorn jämfört med själv, n jämfört med n2 tal. Denna situation har man exempelvis vid hållfasthetsberäkningar
med den så kallade finita elementmetoden där
är lättare att beräkna än .
Den bakomliggande idén är enkel. Gissa en lösning 0 , mata algoritmen och ut kommer 1 som är en bättre lösning till ekvationssystemet. Denna skickar man naturligtvis in igen å igen ända till dess ingen förändring sker. Detta kallas att iterera (efter engelskans
iterate = upprepa). Man får alltså en sekvens av allt bättre i . Den stora fördelen är den tidigare nämnda att endast
behöver
levereras från problembeskrivningen samt att man själv kan bestämma lösningens noggrannhet genom att avbryta itererandet när
Ε, där Ε är ett på förhand bestämt litet tal.
detta är uppnått. Vanligtvis testas noggrannheten efter iteration i som
i
i 1
När(/om) detta inträffar har vi konvergens.
Naturligtvis finns det nackdelar. Den största är ju frågan hur många iterationer som behövs innan önskad noggrannhet Ε uppnås?
Detta avgörs av det så kallade konditionstalet för , som i sin tur beror på hur utspridda egenvärdena till är. Genom så kallad
förkonditionering med en matris som ''liknar'' 1 ökar prestandan rejält. Stor möda ägnas att finna dylika ''billiga'' förkonditioneringar. Detta sökande tillsammans med utveckling av effektiva iterativa metoder har varit mycket framgångsrikt på senare år. Redan
för små problem har de blivit konkurrenskraftiga och för de riktigt stora har de varit den enda framkomliga vägen. En enkel medlem
i denna klass av metoder brukar kallas steepest descent och kan enkelt beskrivas med en analogi från verkliga livet.
Antag att du befinner dig på en alpsluttning och vill komma ned i dalen (minimering!). På grund av dimma är sikten begränsad till
några meter. En strategi är då att gå i den riktning som lutar som mest nedåt. Fortsätt att gå rakt i denna riktning ända till det börjar
luta uppåt. Tag ut ny riktning nedåt. Upprepa (iterera!) till dess att du kommer till en punkt där det är uppförsbacke i alla riktningar.
Du har då kommit till dalens lägsta punkt!
Nu klär vi denna strategi i något mera precis matematisk dräkt. Metoden med steepest descent är som sagt helt enkelt en metod att
söka lokalt minimum till en funktion. Strategin är att från aktuell punkt i 1 förflytta sig i negativa gradientens riktning
. Man kan med hjälp av figurerna nedan dra sig till minnes att detta innebär i en riktning vinkelrät mot den
i 1
i 1
nivåkurva man befinner sig på. Det optimala steget i denna riktning är ju att gå så länge funktionen avtar. Detta ges av Αi i 1 . Från
denna nya punkt i
Αi i 1 tar man så ut en ny riktning och så vidare till dess man kommit till en punkt där ingen negativ
i 1
gradientriktning kan hittas. Vi har då kommit till minpunkten. Metoden fungerar speciellt bra på kvadratiska funktioner, vilket det är
frågan om här. Steglängden Αi , vars uppbyggnad vi inte fördjupar oss i för tillfället, är utformad för att vara optimal i detta fall, det
vill säga vi tar hela den aktuella nerförsbacken i ett steg. Vetskapen att en kvadratisk funktion bara har ett enda lokalt och därmed
globalt minimum ökar tryggheten ytterligare.
brukar kallas residualen (felet) eller residualvektorn eftersom den är just skillnaden mellan vänsterled och
Vektorn i 1
i 1
2
högerled i ekvationssystemet. Så den utgör därmed också ett mått på lösningens kvalitet eftersom
0 då vi har konvergens. Detta illustreras i nästa exempel.
32
Matriser och Mathematica
HH/ITE/BN
Nivåkurvor
5
30
30
4
50
3
y
58
61
2
61.6
6
1 20
55
10
0
0
1
30
2
42
3
4
5
x
Exempel: Dressera Mathematica till att lösa ekvationssystemet
4 2
2 7
x
y
20
med steepest descent. Figurerna ovan vilar på
22
detta exempel där man tydligt ser vägen ned i dalen med optimal steglängd.
Lösningsförslag: Vi kör igång direkt, starta i
4 2
;
2 7
20, 22 ;
0
0, 0
häng mé
0.0, 0.0 ; track
;
1, 1 ;
10 6 ,
While .
. ;
.
;
. .
Α ;
Print
, Α, , Α , .
AppendTo track,
;
Α
;

2.62004, 2.88204 , 0.131002, 20., 22. , 2.62004, 2.88204 , 884.
3.73566, 1.86783 , 0.297043, 3.75578, 3.41434 , 1.11563, 1.01421 , 25.7636
3.90881, 2.05829 , 0.131002, 1.32168, 1.45385 , 0.173143, 0.190457 , 3.86052
3.98253, 1.99127 , 0.297043, 0.248197, 0.225634 , 0.0737253, 0.067023 , 0.112513
3.99397, 2.00385 , 0.131002, 0.0873422, 0.0960764 , 0.011442, 0.0125862 , 0.0168593
3.99885, 1.99942 , 0.297043, 0.0164019, 0.0149108 , 0.00487207, 0.00442915 , 0.000491354
3.9996, 2.00025 , 0.131002, 0.00577193, 0.00634912 , 0.000756133, 0.000831746 , 0.0000736265
 3.99992, 1.99996 , 0.297043, 0.0010839, 0.000985368 , 0.000321966, 0.000292697 , 2.1458 10 6 
 3.99997, 2.00002 , 0.131002, 0.000381433, 0.000419576 , 0.0000499684, 0.0000549652 , 3.21535 10 7 
Vilken väg tog vi?
track

0. 2.62004 3.73566 3.90881 3.98253 3.99397 3.99885 3.9996 3.99992 3.99997

0. 2.88204 1.86783 2.05829 1.99127 2.00385 1.99942 2.00025 1.99996 2.00002
Stegen är markerade med röda punkter i figuren ovan. Lägg speciellt märke till att vi väsentligen är framme redan efter några
iterationer. Till slut stillar vi vår oro
Solve
x
4, y
. x, y
2
Avslutningsvis nöjer vi oss med att nämna att liknande iterativa algoritmer kan konstrueras för att bestämma inverser samt egenvärden och egenvektorer.