teoria, Apuntes de Cálculo

¡Descarga teoria y más Apuntes en PDF de Cálculo solo en Docsity!

Càlcul numèric i matricial

Pep Mulet

Índex

Capítol 1

Preliminars

1.1 Notació

Els escalars es representaran habitualment per lletres gregues minúscules: α, β, γ,.. .. Interpretarem els vectors com a vectors columna quan no s’especifique el contrari i els representarem habitualment per lletres romanes minúscules: x, y, a, b,.. .. Les matrius les representarem per lletres romanes majúscules: A, B, X,

.. .. L’element i d’un vector v es representa per vi, v(i) o ai si s’ha especificat prèviament que v = (ai). L’element (i, j) (fila i, columna j) d’una matriu A es representa per Aij o A(i, j) o aij , si hem especificat prèviament que A = (aij ). La notació que utilitzarem per referir-nos a les files, les columnes i les submatrius d’una matriu és de tipus matlab , per exemple - A(i, :) ≡ fila i-èsima de A. - A(:, i) ≡ columna i-èsima de A. - en general, si I i J són subconjunts dels índex de files i columnes, respectivament, aleshores A(I, J) denota la submatriu de A formada a partir de les files de I i les columnes de J.

Així mateix, si v 1 ,... , vn ∈ Rm, aleshores representarem la matriu, les columnes de la qual són v 1 ,... , vn, per [v 1... vn]. De la mateixa manera, si u 1 ,... , um ∈ Rn^ són vectors fila , aleshores repre- sentarem la matriu, les files de la qual són u 1 ,... , um, per

  

u 1 .. . um

1.2 Errors numèrics i nombres en punt flotant

Quan aproximem una quantitat x per una altra xˆ, els errors que es cometen en aquesta aproximació són:

• error absolut: x − xˆ
• error relatiu: x−x ˆx(suposant x 6 = 0).

Els nombres en punt flotant (en una base donada) venen caracteritzats essencialments com aquells que tenen mantissa i exponent de longituds finites característiques. Els processadors actuals implementen dos tipus de nombres en punt flotant, que s’anomenen precissió simple i precissió doble. Aquests nombres en base dos venen caracteritzats per:

1.3. OPERACIONS AMB VECTORS I MATRIUS 4

• precissió simple : mantissa de 24 bits, rang d’exponents (en base 2) entre -127 i 128, codificació en 32 bits.
• precissió doble : mantissa de 53 bits, rang d’exponents (en base 2) entre -1023 i 1024, codificació en 64 bits.

Com que n’hi ha una quantitat finita de nombres en punt flotant (una volta fixada la longitud de la mantissa i de l’exponent) n’hi ha un mínim m i un màxim M dels nombres en punt flotant positius. Qualsevol real positiu inferior a m produeix una situació que s’anomena underflow i qualsevol real positiu superior a M produeix una situació de overflow. Un nombre real en l’interval [m, M ] es pot aproximar òptimament per un nombre en punt flotant, per un procés anomenat arrodoniment. Per exemple si, per simplicitat, treballem en base 10 i amb una mantissa de 4 dígits, el nombre 12.345 té mantissa 0.12345 i exponent 2. Com que la mantissa té longitud major que 4, no es tracta d’un nombre en punt flotant com els que hem establert. El seu arrodoniment f l(x) es calcula retenint les quatre primeres xifres significatives de la mantissa, augmentant una unitat l’última si aquesta és superior a 4:

x = 0. 12345 × 102 ⇒ fl(x) = 0. 1235 × 102 = 12350.

Aquesta estratègia d’aproximació té la propietat d’assegurar la fita superior de l’error relatiu: ∣ ∣ ∣ ∣

x − fl(x) x

∣ ≤^

base^1 −longitud mantissa

Les operacions aritmètiques exactes amb nombres en punt flotant no donen obligatòriament dins del conjunt de nombres en punt flotant amb el qual es treballa. Les operacions aritmètiques en punt flotant corresponents a aquestes es solen modelar com l’arrodoniment del resultat de l’operació exacta, amb la qual cosa s’obté una fitació superior de l’error relatiu de l’operació en punt flotant com a aproximació de la seua operació corresponent.

1.3 Operacions amb vectors i matrius

Producte com a combinació lineal

La multiplicació d’una matriu per un vector s’efectua habitualment mitjançant:

(A ∗ x)(i) = A(i, :) ∗ x =

∑^ n

k=

A(i, k) ∗ x(k), i = 1,... , m,

on A és m × n i x és n × 1. La següent igualtat sol resultar molt útil:

A ∗ x = x(1) ∗ A(:, 1) + x(2) ∗ A(:, 2) + · · · + x(n) ∗ A(:, n),

és a dir,

1.4. MATRIUS PER BLOCS 6

1.4 Matrius per blocs

Podem considerar una matriu una matriu de matrius , efectuant una partició de les files i columnes de la matriu en diversos grups. Per ilustrar aquest fet considerarem particions de dos grups, encara que, en general, el nombre potser arbitrari. Per exemple, la matriu A anterior es pot particionar de la següent manera:    

[

A 11 A 12

A 21 A 22

]

A 11 =

 A 12 =

 A 21 =

[

]

A 22 =

[

]

Operacions per blocs

Considerem la matriu m × n

A =

[

A 11 A 12

A 21 A 22

]

on Aij és una matriu mi × nj , on m 1 + m 2 = m, n 1 + n 2 = n. La transposició de A es pot calcular com:

AT^ =

[

AT 11 AT 21

AT 12 AT 22

]

Si α és un escalar, aleshores es verifica:

α ∗ A =

[

α ∗ A 11 α ∗ A 12 α ∗ A 21 α ∗ A 22

]

Si B és una altra matriu m × n amb una partició de dimensions idèntiques a la de A

B =

[

B 11 B 12

B 21 B 22

]

on Bij és una matriu mi × nj , on m 1 + m 2 = m, n 1 + n 2 = n, es verifica

A + B =

[

A 11 + B 11 A 12 + B 12

A 21 + B 21 A 22 + B 22

]

Si B és una matriu n × l amb una partició de dimensions compatibles pel producte per A, és a dir,

B =

[

B 11 B 12

B 21 B 22

]

on Bij és una matriu ni × lj , on n 1 + n 2 =, l 1 + l 2 = l, es verifica

A ∗ B =

[

A 11 ∗ B 11 + A 12 ∗ B 21 A 11 ∗ B 12 + A 12 ∗ B 22

A 21 ∗ B 11 + A 22 ∗ B 21 A 21 ∗ B 12 + A 22 ∗ B 22

]

Resumint:

Les operacions amb matrius per blocs s’efectuen com si els blocs foren escalars, amb la precau- ció que cal assegurar que les operacions resultants entre els blocs es puguen efectuar i que el producte matricial no és commutatiu.

1.5. MATRIUS ORTOGONALS 7

1.5 Matrius ortogonals

Definició 1.5.1. Una matriu A es diu ortogonal si AT^ A = AAT^ = I o, equivalentment, si A és invertible i AT^ = A−^1.

Proposició 1.5.2. Siga A una matriu n × n. Les següents condicions són equivalents:

1. A és ortogonal.
2. Les columnes de A formen una base ortonormal de Rn.
3. Les files de A formen una base ortonormal de Rn.

Prova. L’element (i, j) de AT^ A és el producte escalar de la fila i de AT^ i la columna j de A, és a dir, de la columna i i la columna j de A. Si la matriu és ortogonal, aquest producte dóna 0 si i 6 = j i 1 si i = j, per tant les columnes de A són una base ortonormal de Rn^ (perquè té n elements). El recíproc es prova de manera similar i l’equivalència d’1 i 3, utilitzant AAT^.

1.6 Valors i vectors propis

Definició 1.6.1. Un vector propi d’una matriu A, associat al valor propi λ, és un vector x 6 = 0 tal que

Ax = λx.

Un escalar λ és valor propi de la matriu A n × n si i solament si és arrel del polinomi característic

det(A − tI) = (−1)n(tn^ − traça(A)tn−^1 + · · · + (−1)n^ det(A)).

Com que tot polinomi real té arrels en C, en deduïm que tota matriu real té valors propis, no necessària- ment reals. Per exemple, la matriu (^) [ 0 1 − 1 0

]

té com a únics valors propis les arrels i, −i del seu polinomi característic ∣ ∣ ∣ ∣

−t 1 − 1 −t

∣ =^ t

Definició 1.6.2. Una matriu A es diu diagonalitzable si existeix una matriu invertible P tal que P −^1 AP = D, amb D diagonal o, equivalentment, si AP = P D. (1.2)

Com que una matriu és invertible si i solament si les seues columnes formen una base, la igualtat (1.2) implica que una matriu A és diagonalitzable si i solament si existeix una base formada per vectors propis: les columnes de P són aquesta base de vectors propis i l’element i-èsim de la diagonal de D és el valor propi de la i-èsima columna de P.

Definició 1.6.3. El radi espectral d’una matriu és el màxim dels valors absoluts del valors propis de la matriu, és a dir ρ(A) = max λv.p.A

|λ|. (1.3)

Com veurem més endavant, el radi espectral està relacionat amb la grandària de la matriu. Aquesta relació és molt clara per a matrius simètriques, però no en general, ja que es pot donar que el radi espectral d’una matriu siga 0 (quins valors propis tindrà la matriu ...?) sense que la matriu ho siga, per exemple: [ 0 a 0 0

]

té el radi espectral 0, independentment de a, el qual pot ser arbitràriament gran. Adoneu-se’n que la matriu d’aquest exemple no és simètrica.

1.8. NORMES DE MATRIUS 9

3. ‖A + B‖ ≤ ‖A‖ + ‖B‖ ( desigualtat triangular ).

Com que una matriu m × n es pot identificar amb un vector de mn components, a partir d’una norma vectorial es pot construir una norma matricial; però aquesta construcció no sol donar normes útils. Una excepció és:

Exemple 1.8.2. La norma de Frobenius d’una matriu A es defineix com a:

‖A‖F =

( (^) ∑m

i=

∑^ n

j=

|aij |^2

Amb la identificació natural de la matriu

A =

a 11... a 1 n

............... am 1... amn

amb el vector a = (a 11 ,... , a 1 n, a 21 ,... , a 2 n,... , am 1 ,... , amn),

tenim que ‖A‖F = ‖a‖ 2.

Una diferència fonamental entre vectors i matrius és que aquestes últimes admeten l’operació producte, doncs resultaria interessant poder relacionar el producte matricial amb les normes. Idealment desitjariem que es verificara ‖AB‖ = ‖A‖‖B‖,

però aquest requeriment és massa fort.

Definició 1.8.3. Una norma matricial ‖ · ‖ satisfà la propietat multiplicativa si

‖AB‖ ≤ ‖A‖‖B‖.

Cal tenir una mica de compte amb les normes involucrades en la desigualtat anterior, ja que, encara que les denotem de la mateixa manera, rigorosament poden ser diferents si les matrius no són quadrades. Tot seguit veurem en quin sentit es verifica aquesta propietat.

Proposició 1.8.4. La norma de Frobenius satisfà la propietat multiplicativa.

Prova. Cal veure que ‖AB‖F ≤ ‖A‖F ‖B‖F per a matrius qualssevol A, B (amb grandàries compatibles amb el producte matricial):

‖AB‖^2 F =

i

j

(AB)(i, j)^2 =

i

j

k

A(i, k)B(k, j)

i

j

k

A(i, k)^2

k

B(k, j)^2

i

k

A(i, k)^2

j

k

B(k, j)^2 = ‖A‖^2 F ‖B‖^2 F ,

on hem utilitzant la desigualtat de Cauchy-Schwartz per a la fita: ∑

k

A(i, k)B(k, j) = A(i, :)B(:, j)

≤ ‖A(i, :)‖ 2 ‖B(:, j)‖ 2 =

k

A(i, k)^2

k

B(k, j)^2

Resulta molt útil identificar A amb l’operador x 7 → Ax, aleshores es pot considerar la grandària de la matriu A com el major factor d’amplificació de la norma al passar de x a Ax. Més precisament:

1.9. TRANSFORMACIONS ELEMENTALS 10

Definició 1.8.5. Siga ‖ · ‖α una norma vectorial, on α = 1, 2 , ∞, p. Definim la norma matricial associada com a:

‖A‖α = sup x 6 =

‖Ax‖α ‖x‖α

Noteu que si A no és quadrada, A ∈ Rm×n, aleshores Ax ∈ Rm^ i x ∈ Rn, amb la qual cosa ‖Ax‖α fa referència a la norma ‖ · ‖α en Rm^ i ‖x‖α a la norma ‖ · ‖α en Rn.

Nota 1.8.6. Una matriu ortogonal A satisfà que ‖A‖ 2 = 1: si x ∈ Rn, x 6 = 0,

‖Ax‖^22 = (Ax)T^ (Ax) = xT^ AT^ Ax = ‖x‖^22 ,

d’on ‖Ax‖ 2 = ‖x‖ 2 , la qual cosa implica que

‖A‖ 2 = max x 6 =

‖Ax‖ 2 ‖x‖ 2

Proposició 1.8.7. Siga ‖ · ‖ una norma vectorial i A ∈ Rn×n, aleshores existeix x∗^ ∈ Rn^ amb ‖x∗‖ = 1 tal que ‖Ax∗‖ = ‖A‖, és a dir,

sup x 6 =

‖Ax‖ ‖x‖

= ‖A‖ = max x 6 =

‖Ax‖ ‖x‖

= max ‖x‖=

‖Ax‖.

Proposició 1.8.8. (Propietat multiplicativa) Siga A una matriu m × n i representem per ‖ · ‖ tant una norma vectorial com la seua norma matricial associada, aleshores:

• ‖Ax‖ ≤ ‖A‖‖x‖, ∀x ∈ Rn.
• ‖AB‖ ≤ ‖A‖‖B‖, ∀B matriu n × p.

Proposició 1.8.9. Siga A una matriu m × n. Aleshores

‖A‖ 1 = max 1 ≤j≤n

∑m

i=

|aij | = max(‖A(:, 1)‖ 1 ,... , ‖A(:, n)‖ 1 ) (1.4)

‖A‖∞ = max 1 ≤i≤m

∑n

j=

|aij | = max(‖A(1, :)‖ 1 ,... , ‖A(m, :)‖ 1 ). (1.5)

La norma 2, ‖ · ‖ 2 també admet una caracterització que permet el seu càlcul. Aquesta caracterització utilitza els valors propis de la matriu AT^ A.

Proposició 1.8.10. La norma 2 d’una matriu real m × n A coincideix amb l’arrel quadrada del radi espectral de AT^ A, és a dir,

‖A‖ 2 =

ρ(AT^ A).

Si la matriu és simètrica, aleshores ‖A‖ 2 = ρ(A).

1.9 Transformacions elementals

Les transformacions elementals efectuades a les equacions d’un sistema lineal deixen invariant la seua solu- ció. Mitjançant la matriu ampliada, aquestes transformacions es poden considerar en un contexte matricial. Les transformacions elementals (per files) són de tres tipus:

1. Intercanvi de dues files.
2. Producte d’una fila per un nombre no nul.

1.10. TRANSFORMACIONS ELEMENTALS ORTOGONALS 12

Aquest problema concret es pot solucionar amb una tècnica que s’anomena pivotatge parcial (veure més endavant el capítol 2), que consisteix en permutar files de manera que el pivot siga el màxim possible, amb la qual cosa els multiplicadors són menors que 1. En el cas anterior convertiríem el sistema en

x + y = 2 10 −^8 x + y = 1

En aquest cas la transformació elemental de tipus III que produeix un 0 en la posició (2, 1) és

[ 1 0 − 10 −^8

]

que té norma ≈ 1 , doncs és una matriu “inofensiva” des del punt de vista numèric.

1.10 Transformacions elementals ortogonals

Per evitar el problema de la possible amplificació d’errors de les transformacions elementals de tipus III trac- tarem de trobar transformacions ortogonals que facen una tasca semblant a aquelles, és a dir, que produisquen zeros en posicions determinades. Per tractar-se de transformacions ortogonals, la seua norma (euclídia) és 1, per tant, el mínim possible. No hi haurà, aleshores, problemes d’amplificació d’errors. Com tot no poden ser avantatges, les transformacions ortogonals que anem a veure tenen un cost computacional més alt que les transformacions anteriors.

Transformacions de Givens

L’objectiu de les transformacions de Givens és efectuar un gir d’un angle a determinar per a que un vector donat es transforme en un altre vector amb una de les seues components 0. Per simplificar, suposem que tenim un vector (a, b) i que volem calcular un gir que el transforme en (c, 0)

Quina matriu te aquesta transformació? En general, les matrius dels girs són del tipus [ cos α − sin α sin α cos α

]

on α és l’angle de gir. La matriu que busquem és

1 c

[

a b −b a

]

c =

a^2 + b^2.

Comproveu com a exercici que és ortogonal i que fa l’efecte dessitjat.

1.10. TRANSFORMACIONS ELEMENTALS ORTOGONALS 13

Vegem com actuaríem en el cas general. Suposem que tenim que transformar  

a b c

t 0 0

mitjançant transformacions de Givens. El que cal fer és produir un primer zero on és b: [ a/r b/r −b/r a/r

] [

a b

]

[

r 0

]

(r =

a^2 + b^2 ) ⇒  

a/r b/r 0 −b/r a/r 0 0 0 1

a b c

r 0 c

El següent pas produeix un zero on és c, de la manera següent: [ r/t c/t −c/t r/t

] [

r c

]

[

t 0

]

(t =

r^2 + c^2 ) ⇒  

r/t 0 c/t 0 1 0 −c/t 0 r/t

r 0 c

t 0 0

Com exercici , compareu el nombre d’operacions d’una transformació de Givens amb el d’una transfor- mació elemental de tipus III.

Transformaciones de Householder

L’objectiu d’aquestes transformacions és calcular una simetria per convertir un vector donat en altre que es trobe sobre l’eix de les x, és a dir, que tinga totes les components zero excepte la primera. Vegem gràficament en la figura 1.1 com calcular aquesta simetria. Les relacions que es dedueixen de la

Figura 1.1: Transformació de Householder

gràfica són les següents: x = a + b x′^ = a − b b = (x, v)v

on (x, v) és el producte escalar de x amb un vector v ortogonal a l’hiperplà de simetria. Deduïm l’equació de la simetria, coneixent el vector v ortogonal a l’hiperplà de simetria:

x′^ = x − 2 b = x − 2(x, v)v = (I − 2 vvT^ )x, (1.7)

1.11. PROBLEMES 15

satisfà que ‖x‖ 1 = 1 i ‖Ax‖ 1 = ‖A‖ 1 = 12, ja que

Ax =

‖Ax‖ 1 = 12.

2. ‖A‖∞ = major suma de valors absoluts de files = max{ 7 , 10 , 8 , 4 } = 10. El vector

x =

satisfà que ‖x‖∞ = 1 i ‖Ax‖∞ = ‖A‖∞ = 10, ja que

Ax =

‖Ax‖∞ = 10.

Problema 2. Dibuixeu en R^2 les boles unitat B(0, 1) = {x : ‖x‖ ≤ 1 } per a les normes ‖ · ‖ 1 , ‖ · ‖∞, ‖ · ‖ 2.

Solució.

1. ‖ · ‖ 1

B 1 (0, 1) = {(x 1 , x 2 ) : |x 1 | + |x 2 | ≤ 1 } = {(x 1 , x 2 ) : x 1 + x 2 ≤ 1 , x 1 − x 2 ≤ 1 , −x 1 + x 2 ≤ 1 , −x 1 − x 2 ≤ 1 }.

B 2 (0, 1) = {(x 1 , x 2 ) :

x^21 + x^22 ≤ 1 } = {(x 1 , x 2 ) : x^21 + x^22 ≤ 1 }.

1.11. PROBLEMES 16

B∞(0, 1) = {(x 1 , x 2 ) : max(|x 1 |, |x 2 |) ≤ 1 } = {(x 1 , x 2 ) : x 1 ≤ 1 , −x 1 ≤ 1 , x 2 ≤ 1 , −x 2 ≤ 1 }.

Problema 3. Proveu que xT^ x = ‖x‖^22 , x ∈ Rn.

Problema 4. Proveu que ‖A‖F = traça(AT^ A) (^12) .

Problema 5. Siguen P ∈ Rm×m, Q ∈ Rn×n^ matrius ortogonals i A ∈ Rm×n. Proveu que

1. ‖ P AQ ‖ 2 =‖ A ‖ 2.
2. ‖ P AQ ‖F =‖ A ‖F.

Solució. 1. Com que Q és invertible i ‖P y‖ 2 = ‖y‖ 2 i ‖Qx‖ 2 = ‖x‖ 2 per a y ∈ Rm^ i x ∈ Rn, tenim que {x ∈ Rn^ : ‖x‖ 2 = 1} = {Qx ∈ Rn^ : ‖x‖ 2 = 1}, d’on es dedueix:

‖P AQ‖ 2 = max ‖x‖ 2 =

‖P AQx‖ 2 = max ‖Qx‖ 2 =

‖A(Qx)‖ 2 = max ‖z‖ 2 =

‖Az‖ 2 = ‖A‖ 2.

2. Pel problema 4 i tenint en compte que P T^ = P −^1 , QT^ = Q−^1 ,

‖P AQ‖^2 F = traça

(P AQ)T^ P AQ

= traça

QT^ AT^ (P T^ P )AQ

traça

QT^ AT^ AQ

= traça

Q−^1 AT^ AQ

= traça

AT^ A

= ‖A‖^2 F ,

ja que la traça és invariant per semblança.

Problema 6. Demostreu que una matriu diagonal és definida positiva si i solament si tots els elements de la diagonal són positius.

Problema 7. Si A és una matriu n × n, demostreu que són equivalents:

1. Existeix una matriu regular (invertible) P i una matriu D diagonal tal que P −^1 AP = D.
2. Existeix una base v 1 ,... , vn de Rn^ tal que vi és vector propi de A.

Prova. Per demostrar l’equivalència, podeu seguir els passos següents:

1. Utilitzeu l’equivalència P −^1 AP = D si i sols si AP = DP.
2. Deduïu igualtats vectorials, igualant les columnes de la igualtat matricial AP = DP.

Problema 8. Demostreu que el radi espectral d’una matriu quadrada està fitat superiorment per la norma de la matriu, per a normes que satisfacen la propietat multiplicativa.

Problema 9. Si A = [a 1 · · · an] ∈ Rm×n^ ; ai ∈ Rm^ i B = [b 1 · · · bn] ∈ Rp×n^ ; bi ∈ Rp. Proveu que

A · BT^ =

∑^ n

k=

ak · bTk.

2.1. RESOLUCIÓ DE SISTEMES TRIANGULARS 18

El camp d’aplicació d’aquest resultat és aquell en que ‖∆A‖, ‖∆b‖ ≈ 0 , per la qual cosa podem suposar r ≈ 0. En aquest cas, podem assegurar que el factor d’ampliació de l’error relatiu (el denominador es pot interpretar com l’error relatiu en les dades i el numerador és l’error relatiu en la solució)

‖∆x‖ ‖x‖ ‖∆A‖ ‖A‖ +^

‖∆b‖ ‖b‖

≤ κ(A)

La manera d’interpretar i aplicar aquest resultat és, per exemple, la següent: si

‖∆A‖ ‖A‖

≤ 10 −^15

‖∆b‖ ‖b‖

≤ 10 −^15

r ≈ 0 κ(A) = 10d

aleshores ‖∆x‖ ‖x‖

≤ 2 · 10 −(15−d),

és a dir, d’una exactitud en les dades de 15 xifres significatives en les dades passem a una exactitud en la solució de 15-d xifres significatives, és a dir, podem tenir una pèrdua de d xifres significatives a l’hora d’obtenir la solució. Diguem podem perquè potser que en alguns casos la pèrdua siga menor, ja que sols tenim una afitació superior. Adoneu-se’n també que κ(A) ≥ 1 i que el que sí que assegura aquest resultat és que si una matriu és ben condicionada, és a dir, κ(A) és petit, la pèrdua de xifres significatives és mínima.

2.1 Resolució de sistemes triangulars

Resoldre sistemes lineals (compatibles i determinats) no és una tasca trivial, encara que hi ha casos en els quals si que ho és, com, per exemple, el que té una matriu diagonal. No tan trivial, però quasi, resulta la resolució de sistemes amb matriu triangular superior amb l’algoritme de substitució cap endarrere :

Input: A triangular superior n × n, b n × 1 for i = n,... , 1 for j = i + 1,... , n b(i) = b(i) − A(i, j)b(j) end b(i) = b(i)/A(i, i) end Output: b

Algoritme 1: Algoritme de substitució cap endarrere. Després de l’execució d’aquest algoritme, la solució queda en la variable b.

De manera dual, la resolució de sistemes amb matriu triangular inferior es pot efectuar amb l’algoritme de substitució cap endavant Com es dedueix, el principi en el qual es basen aquests algoritmes és la resolució de l’equació més senzilla de totes (amb una única incògnita) primer.

2.2. ELIMINACIÓ DE GAUSS 19

Input: A triangular inferior n × n, b n × 1 for i = 1,... , n for j = 1,... , i − 1 b(i) = b(i) − A(i, j)b(j) end b(i) = b(i)/A(i, i) end Output: b

Algoritme 2: Algoritme de substitució. Després de l’execució d’aquest algoritme, la solució queda en la variable b.

Input: A matriu n × n, b vector n × 1 for k = 1,... , n − 1 for i = k + 1,... , n m = A(i, k)/A(k, k) for j = k + 1,... , n A(i, j) = A(i, j) − m A(k, j) end b(i) = b(i) − mb(k) end end

Algoritme 3: Algoritme de l’eliminació de Gauss.

2.2 Eliminació de Gauss

La manera d’obtenir sistemes amb la mateixa solució que un donat Ax = b és efectuar transformacions ele- mentals a les equacions o, equivalentment, pre-multiplicant per una matriu invertible M ambdos membres de l’equació per obtenir el sistema equivalent Ax˜ = ˜b, on A˜ = M A, ˜b = M b. L’algoritme 3 de l’eliminació de Gauss essencialment calcula una matriu triangular inferior M , amb uns en la diagonal, tal que M A = U és triangular superior. El seu cost computacional (el nombre d’operacions) és 23 n^3 + O(n^2 ). A l’eixida d’aquest algoritme, en la part triangular superior de A obtenim la part triangular superior de U = M A i en b obtenim M b, amb la qual cosa podem aplicar als nous A i b l’algoritme de la substitució cap endarrere. Noteu que guardem en la part estrictament triangular inferior de A els multiplicadors utilitzats en el procés d’eliminació. Es pot comprovar que la matriu triangular inferior amb uns en la diagonal formada a partir d’aquestes entrades és precisament L = M −^1 , una matriu triangular inferior amb uns en la diagonal. Resulta aleshores que LU = M −^1 M A = A, la qual cosa és la interpretació algebraica de l’eliminació de Gauss. De fet, qualsevol factorització A = LU amb L triangular inferior i U triangular inferior permet una resolució ràpida de sistemes que tinguen A com matriu de coeficients, mitjançant la solució de sistemes amb matrius L (utilitzant substitució cap endavant) i U (utilitzant substitució cap endarrere):

Ax = b ⇔ (LU )x = b ⇔ L (U x) ︸ ︷︷ ︸ y

= b

Ly = b, U x = y

Exemple 2.2.1. Considerem el sistema que té per matriu ampliada la següent:

 