1.4. Linj¨ ara ekvationssystem

(1)

1.3. Matrisr¨ akning

I kursens förra del beskrevs hur man löser ett lineärt ekvationssystem i matrisform: Ax = b. Vi kommer nu att studera problemet i detalj, men först skall vi behandla multiplikation av en matris och en vektor, samt multiplikation av tv˚a matriser. I princip finns dessa operationer färdigt implementerade i MATLAB, men för att vi bättre skall lära oss att handskas med matriser, s˚a lönar det sig att studera uppställningen av matrisproblem och operationer med matriser mera i detalj.

Elementen i en matris kan konstrueras p˚a olika sätt. Det finns t.ex. speciella matriser, som kan genereras utg˚aende fr˚an ett matematiskt uttryck för matriselementen. Ett exempel är Hilberts matris:

a_ij = 1

i + j − 1, som kan genereras med programmet

A = zeros(n,n);

for i = 1:n for j = 1:n

A(i,j) = 1/(i+j-1);

end end

där matrisen först blivit nollställd med zeros(n,n).

(2)

Detta program utnyttjar inte det faktum, att Hilberts matris är symmetrisk: aij = a_ji ∀(i, j). Om vi beaktar detta, räcker det med att beräkna elementen p˚a diagonalen och (t.ex.) elementen ovanför diagonalen, resten följer av symmetrivillkoret:

A = zeros(n,n);

for i = 1:n for j = i:n

A(i,j) = 1/(i+j-1);

A(j,i) = A(i,j);

end end

Det finns faktiskt en inbyggd funktion hilb(n) i MATLAB, som alstrar en n–radig Hilbert-matris, som ocks˚a kan anv¨andas ( type hilb visar en listning av programmet).

Ett annat exempel ¨ar en matris, som uppbyggs av binomialkoefficienterna (Pascals triangel)

n k

=

( n!

k!(n−k)! om 0 ≤ k ≤ n 0 i annat fall

(3)

I detta fall kan vi konstruera matriselementen t.ex. genom formeln p_ij =

i − 1 j − 1

.

Att räkna ut binomialkoefficienterna direkt via deras definition är arbetsdrygt. Istället lönar det sig att använda rekursionslikheten

p_ij = p_i−1,j−1 + p_i−1,j.

Ett program, som utnyttjar denna relation ser vi nedan (jfr pascal i MATLAB):

function p=binomial(n) p = zeros(n,n)

p(:,1) = ones(n,1);

for i=2:n for j=2:i

p(i,j) = p(i-1,j-1) + p(i-1,j);

end end

(4)

Vandermondes matris har studerats tidigare i samband med interpolation:

V = 2 66 4

1 x₁ x²₁ x³₁ 1 x₂ x²₂ x³₂ 1 x₃ x²₃ x³₃ 1 x₄ x²₄ x²₄

3 77 5

Den kan enklast konstrueras kolumnvis p˚a följande sätt (jfr även vander i MATLAB):

n = length(x);

V(:,1) = ones(n,1);

for j = 2:n

V(:,j) = x.*V(:,j-1);

end

I en cirkulant–matris skiftas raderna cykliskt:

C = 2 66 4

a₁ a₂ a₃ a₄ a₄ a₁ a₂ a₃ a₃ a₄ a₁ a₂ a₂ a₃ a₄ a₁

3 77 5

(5)

En s˚adan matris kan vi t.ex. generera ur definitionen

c_ij = a[(n−i+j) mod n]+1,

eller utnyttja det faktum att c(i, :) f˚as genom att skifta c(i − 1, :) ett steg mot h¨oger. Vi kan l¨att skriva tv˚a Matlab–program som utnyttjar dessa metoder:

function C=circ1(a)

% Invariabel:

% a : en radvektor

% Utvariabel:

% C : en cirkulant, d¨ar C(1,:) = a n = length(a);

C = zeros(n,n);

for i = 1:n for j = 1:n

C(i,j) = a (rem(n-i+j,n)+1);

end end

(6)

function C=circ2(a)

% Invariabel:

% a : en radvektor

% Utvariabel:

% C : en cirkulant, d¨ar C(1,:) = a n = length(a);

C = zeros(n,n);

C(1,:) = a;

for i = 2:n

C(i,:) = [C(i-1,n) C(i-1,1:n-1)];

end

Det andra alternativet ¨ar det kortaste, och samtidigt ocks˚a det snabbaste programmet.

M˚anga av de matriser, som man stöter p˚a i fysiken inneh˚aller matriselement, som är lika med 0. Förekommer de rikligt, brukar man tala om glesa (eng.: sparse) matriser. En triangulär matris har element som är olika noll endast p˚a diagonalen, och under (eller ovanför) densamma (jfr tril, och triu i MATLAB).

(7)

En undre triangul¨ar matris ser t.ex. ut s˚a h¨ar:

L = 2 66 4

a 0 0 0

b c 0 0

d e f 0

g h p q

3 77 5

Tridiagonala matriser har element, som ¨ar olika noll endast p˚a diagonalen, och de tv˚a omgivande subdia- gonalerna:

T = 2 66 4

a b 0 0

c d e 0

0 f g h

0 0 h p

3 77 5

Detta slag av matriser kallas ocks˚a bandmatriser, eftersom de element, som är olika 0, ing˚ar i band, som omger diagonalen. Ett specialfall av en bandmatris är en diagonal matris, som kan konstrueras fr˚an en vektor med hjälp av MATLAB–kommandot diag.

(8)

Om t.ex. d = [1 2 3 4], s˚a ger kommandot D = diag(d) upphov till matrisen

D = 2 66 4

1 0 0 0

0 2 0 0

0 0 3 0

0 0 0 4

3 77 5

Vi skall nu studera matrisoperationer, och börja med multiplikation av en matris och en vektor. Antag att vi vill beräkna y = Ax, där A är en m × n–matris, och x en vektor i ett n–dimensionellt euklidiskt rum. Det normala sättet att utföra denna operation är att beräkna elementen yi, i = 1 : m

y_i = Xn

j=1

a_ijx_j

som skal¨arprodukter av matrisens rader och vektorn x.

(9)

Produkten kan ber¨aknas med MATLAB–programmet

[m,n] = size(A);

y = zeros(m,1);

for i = 1:m for j = 1:n

y(i) = y(i) + A(i,j)*x(j);

end end

Naturligtvis f˚ar man samma resultat med kommandot y = A*x, s˚a vi behöver allts˚a inte beräkna produkten av en matris och en vektor som en dubbelslinga. Men det är änd˚a lärorikt att studera programmet. Vi kan t.ex. vektorisera den inre slingan, som bildar produkten av den i:te raden i A (A(i,:)) och vektorn x:

function y = matvekr(A,x)

% Invariabler:

% A: mxn-matris

% x: kolumnvektor med n element

% Utvariabel:

% y: A*x (radorienterad metod) [m,n] = size(A); y = zeros(m,1);

for i = 1:m

y(i) = A(i,:)*x;

end

(10)

Denna metod att beräkna produkten är radorienterad. Vi kan ocks˚a använda en kolumnorienterad metod.

Principen framg˚ar av f¨oljande exempel:

2

41 2

3 4

5 6

3 5

7 8

= 2

47 · 1 + 8 · 2 7 · 3 + 8 · 4 7 · 5 + 8 · 6

3 5 =

2

47 · 1 7 · 3 7 · 5

3 5 +

2

48 · 2 8 · 4 8 · 6

3

5 = 7 2 41

3 5

3

5 + 8 2 42

4 6

3 5 ,

dvs resultatvektorn är en lineärkombination av kolumnerna i A, där koefficienterna är xj. Detta leder till följande program:

function y = matveks(A,x)

% Invariabler:

% A: mxn-matris

% x: vektor med n element

% Utvariabel:

% y: A*x (kolumnorienterad metod) [m,n] = size(A);

y = zeros(m,1);

for j = 1:n

y = y + A(:,j)*x(j);

end

Denna funktion ¨ar ekvivalent med det ursprungliga programmet, d˚a b˚ada slingorna blivit omkastade.

(11)

Den inre slingan beskriver nu en operation av formen vektor ← vektor · skalär + vektor (y=A(:,j)*x(j)+y), som kallas för saxpy (”Scalar Alpha X Plus Y”) operationen. Denna benämning kommer fr˚an namnet p˚a en subrutin, som ing˚ar i BLAS (Basic Linear Algebra Subroutines) biblioteket. I vektorform kan saxpy operationen i matveks uttryckas

2 66 4

y(1) y(2)...

y(m) 3 77 5 =

2 66 4

A(1, j) A(2, j)

...

A(m, j) 3 77

5 x(j) + 2 66 4

y(1) y(2)...

y(m) 3 77 5

Trots att vartdera programmet använder 2mn flyttalsoperationer, s˚a är kolumnversionen ofta snabbare än radversionen, p˚a grund av att matriser lagras i minnet kolumnvis, dvs kolumnerna följer efter varandra.

Vi skall nu behandla multiplikation av tv˚a matriser. Produkten av en m×p–matris A och en p×n–matris B ¨ar en m × n–matris C = AB, vars element kan ber¨aknas ur formeln

c_ij = Xp

k=1

a_ikb_kj,

d¨ar i och j uppfyller villkoren 1 ≤ i ≤ m och 1 ≤ j ≤ n.

(12)

Detta innebär, att varje element i C är skalärprodukten av en rad i A och en kolumn i B. Matrisen C kan s˚alunda beräknas med MATLAB–programmet

C = zeros(m,n);

for j = 1:n for i = 1:m

for k = 1:p

C(i,j) = C(i,j) + A(i,k)*B(k,j);

end end end

Alternativt kan samma beräkning utföras med kommandot C = A*B (som i allmänhet är den snabbaste metoden). Men vi kan ocks˚a utnyttja vektorisering vid utförandet.

Om vi vektoriserar den innersta slingan i programmet, som beräknar skalärprodukten av den i:te raden i A och den j:te kolumnen i B (allts˚a k ersätts med ’:’), s˚a f˚ar vi funktionen

(13)

function C = matmulp(A,B)

% Invariabler:

% A: mxp matris.

% B: pxn matris.

% Utvariabel:

% C: A*B (via skal¨arprodukt) [m,p] = size(A);

[p,n] = size(B);

C = zeros(m,n);

for j = 1:n

% Ber¨aknar j:te kolumnen i C.

for i=1:m

C(i,j) = A(i,:)*B(:,j);

end end

Men vi vet ocks˚a, att den j:te kolumnen i C ¨ar lika med matrisen A g˚anger den j:te kolumnen i matrisen B. Om vi utnyttjar den kolumnvisa matris–vektoroperationen (i ers¨atts med ’:’), s˚a f˚ar vi funktionen

(14)

function C = matmuls(A,B)

% Invariabler:

% A: mxp matris.

% B: pxn matris.

% Utvariabel:

% C: A*B (saxpy metoden) [m,p] = size(A);

[p,n] = size(B);

C = zeros(m,n);

for j = 1:n

for k = 1:p

C(:,j) = C(:,j) + A(:,k)*B(k,j);

end end

Denna version av matrismultiplikationen anv¨ander saxpy metoden.

(15)

Om den innersta slingan ers¨atts av en enda matrisvektorprodukt f˚ar programmet f¨oljande form:

function C = matmulv(A,B)

% Invariabler:

% A: mxp matris.

% B: pxn matris.

% Utvariabel:

% C: A*B (matrisvektorprodukt) [m,p] = size(A);

[p,n] = size(B);

C = zeros(m,n);

for j = 1:n

C(:,j) = A*B(:,j);

end

(16)

Man kan ytterligare tolka matrismultiplikationen som en summa av yttre produkter. Yttre produkten av en kolumnvektor u med m element och radvektor v med n element definieras genom

uv⁰ = 2 66 4

u₁ u₂ ...

u_m 3 77 5

v₁ v₂ · · · v_n

= 2 66 4

u₁v₁ u₁v₂ · · · u₁v_n u₂v₁ u₂v₂ · · · u₂v_n

... ... . . . ...

u_mv₁ u_mv₂ · · · u_mv_n 3 77 5

Detta kan uppfattas som en vanlig multiplikation av en m × 1 matris och en 1 × n–matris, t.ex.

2 41

2 3

3

5 4 5 6 = 2

4 4 5 6

8 10 12

12 15 18

3 5

Matrismultiplikationsproblemet kan d˚a formuleras som

C = AB = [A(:, 1)|A(:, 2)| · · · |A(:, p)]

2 66 4

B(1, :) B(2, :)

...

B(p, :) 3 77 5 =

Xp k=1

A(:, k)B(k, :).

(17)

Ett exempel p˚a en s˚adan framst¨allning av matrisprodukten ¨ar 2

41 2

3 4

5 6

3 5

7 8

9 10

= 2

41 · 7 + 2 · 9 1 · 8 + 2 · 10 3 · 7 + 4 · 9 3 · 8 + 4 · 10 5 · 7 + 6 · 9 5 · 8 + 6 · 10

3 5

= 2

41 · 7 1 · 8 3 · 7 3 · 8 5 · 7 5 · 8

3 5 +

2

42 · 9 2 · 10 4 · 9 4 · 10 6 · 9 6 · 10

3 5 =

2 41

3 5

3

5 [7 8] + 2 42

4 6

3

5 [9 10]

(18)

Vi kan allts˚a ocks˚a konstruera en annan version av matrismultiplikationsprogrammet, som utnyttjar yttre produkten:

function C = matmulo(A,B)

% Invariabler:

% A: mxp matris.

% B: pxn matris.

% Utvariabel:

% C: A*B (som yttre produkt) [m,p] = size(A);

[p,n] = size(B);

C = zeros(m,n);

for k = 1:p

% Bildar yttre produkten.

C = C + A(:,k)*B(k,:);

end

Denna metod är dock den l˚angsammaste av de fyra metoderna för att multiplicera matriser, vilket man visa genom att skriva ett testprogram, som jämför beräkningstiderna.

(19)

1.4. Linj¨ ara ekvationssystem

Lösningen av ekvationssystem är ett av de vanligaste problemen man stöter p˚a i fysiken. Ofta gäller det ekvationer med m˚anga obekanta, s˚a det är inte konstigt att de första datorerna just användes för lösning av ekvationssystem. Symboliskt kan ett ekvationssystem framställas i matrisform:

Ax = b,

d¨ar A ¨ar en kvadratisk (koefficientmatris) av ordningen n, b en given kolumnvektor med n rader och x en kolumnvektor med n obekanta.

Man kan lösa ekvationssystem p˚a m˚anga olika sätt, t. ex. med Cramers regel ¹, där varje obekant är kvoten av tv˚a determinanter. Denna metod är praktisk bara för ett litet antal obekanta. Om vi t.ex.

skulle till¨ampa metoden p˚a ett system med 30 ekvationer, s˚a skulle vi vara tvungna att ber¨akna 31 st.

30-radiga determinanter. Om man skulle upplösa determinanterna direkt, skulle det behövas 31 · 30! · 29 multiplikationer, och ett motvarande antal additioner. Även p˚a en snabb dator som kan göra 10⁹ flops, skulle det räcka ca 10¹⁹ ˚ar! Man skulle ocks˚a kunna tänka sig att lösa ekvationssystemet genom matrisinversion:

x = A⁻¹y, men det l¨onar sig inte, bl.a. p˚a grund av det finns snabbare metoder.

1efter Gabriel Cramer, schweizisk matematiker (1704-1752), som beskrev metoden i Introduction `a l’analyse des lignes courbes alg´ebraique

(20)

En stor matris kan ta upp mycket utrymme i datorns minne. M˚anga matriser, som används i fysiken, har ett stort antal element lika med noll, och den kallas d˚a gles. Om matrisen är gles, lönar det sig endast att lagra de matriselement som är olika noll. I s˚adana fall m˚aste man använda n˚agon metod för att ange matriselementens läge i den ursprungliga matrisen. Om matrisen är varken gles eller singulär (dvs dess determinant är olika 0), s˚a kan ekvationssystemet lösas enligt en metod som upptäcktes av C.F. Gauss. Vi skall dock börja med att studera triangulära ekvationssystem, som kan lösas p˚a ett enklare sätt.

Antag, att vi har triangul¨art ekvationssystem med tre obekanta:

2

4a₁₁ 0 0

a₂₁ a₂₂ 0 a₃₁ a₃₂ a₃₃

3 5

2 4x₁

x₂ x₃

3 5 =

2 4b₁

b₂ b₃

3 5

Dessa ekvationer kan lösas s˚a, att man först räknar ut x1 ur den första ekvationen, substituerar lösningen i den andra ekvationen, som sedan kan lösas i avseende p˚a x₂, och till sist substituerar x₁ och x₂ i den tredje ekvationen, och beräknar x3:

x₁ = b₁/a₁₁

x₂ = (b₂ − a21x₁)/a₂₂

x₃ = (b₃ − a₃₁x₁ − a₃₂x₂)/a₃₃.

(21)

Denna typ av algoritm kallas fram˚atsubstitution. F¨or att den skall fungera, m˚aste det(A) = a₁₁a₂₂a₃₃ vara olika noll.

Vi skall nu skriva ett MATLAB-program för att lösa systemet Lx = b, där L är en undre triangulär matris. Den i:te ekvationen i detta system kan d˚a skrivas

`_i1x₁ + . . . + `_iix_i = b_i, och l¨osningen kan skrivas

x_i = 0

@b_i − Xi−1

j=1

`_ijx_j 1

A /`_ii.

P˚a basen av denna ekvation kan vi skriva ett programfragment f¨or att ber¨akna x:

for i = 1:n

x(i) = b(i);

for j=1:i-1

x(i) = x(i) - L(i,j)*x(j);

end

x(i) = x(i)/L(i,i);

end