Análise Numérica - PSist2 - Resolução Sistemas Lineares com Matlab

Nota: O Matlab tem o operador \ que permite com uma sintaxe simples resolver sistemas de

equações lineares de qualquer tipo.

Se quisermos resolver o sistema Ax = b basta fazer x=A\b. Se o sistema for determinado o método

que é utilizado é idêntico à factorização LU dada. Vamos, no entanto, ver outras funções que nos

fornecem a própria decomposição.

Faça o download, com copy-past, dos seguintes �cheiros (com funções Matlab) do Moodle (direc-

tório sistemas de equações lineares) para a sua área de trabalho: jacobi.m; seidel.m

1. Pretendemos resolver o sistema de equações Ax = b:

A =

1 1 1 12 2 1 11 2 1 23 1 1 1

, b =

4666

.

1.1. Execute os seguintes comandos Matlab:

A=[1 1 1 1;2 2 1 1;1 2 1 2;3 1 1 1]

b=[4;6;6;6]

xx=A\b

Resolução:

xx = (1.0000 · · · , 1.0000 · · · , 1.0000 · · · , 1.0000 · · · )Nota: Como a resposta aparece com c.d. a solução não é inteira. Se calcularmos a

diferença entre o valor obtido e a solução exacta (1, 1, 1, 1), usando o comando Matlab

dif=xx-[1;1;1;1] vamos obter o valor

dif==

1.0e-015 *

0

-0.1110

-0.2220

0.2220

que mostra que a diferença para os valores exactos ocorre na 16 a c.d..

1.2. Execute o comando Matlab: [L,U,P]=lu(A)

Este comando efectua a decomposição LU de Doolittle da matriz A com pivotagem

parcial sendo L, U e P tais que LU = PA. Indique uma solução aproximada, x, dosistema Ax = b usando esta decomposição. Escreva todos os comandos Matlab que

executou.

Resolução:

[L,U,P]=lu(A);

y=L\(P*b)

x=U\y

Resultados:

L = U =

1.0000 0 0 0 3.0000 1.0000 1.0000 1.0000

0.3333 1.0000 0 0 0 1.6667 0.6667 1.6667

0.3333 0.4000 1.0000 0 0 0 0.4000 0

0.6667 0.8000 -0.5000 1.0000 0 0 0 -1.0000

P =

0 0 0 1

0 0 1 0

1 0 0 0

0 1 0 0

y =

6.0000

4.0000

0.4000

-1.0000

x =

1.0000

1.0000

1.0000

1.0000

1.3. Porque é que LU 6= A?

Resolução:

Porque lu usa a decomposição LU com escolha parcial de pivot. Logo obtém a decom-

posição de uma matriz Aperm constituída pelas linhas de A permutadas.

Nota: Pela análise de P , a 1a linha de Aperm é a 4a linha de A, a 2a linha de Aperm é a

3a linha de A, a 3a linha de Aperm é a 1a linha de A e a 4a linha de Aperm é a 2a linha

de A. De facto:

>> D=L*U >> Aperm=P*A

D = Aperm =

3 1 1 1 3 1 1 1

1 2 1 2 1 2 1 2

1 1 1 1 1 1 1 1

2 2 1 1 2 2 1 1

1.4. Calcule o número de condição da matriz A. Estime a precisão da solução x calculada.

Resolução:

cond(A) = 16.9847 ≈ 101. Como o algoritmo é estável, espera-se que o problema que

na realidade resolvemos aproxime o problema que pretendíamos resolver com cerca de

16 a.s., que é a precisão do Matlab. Como a condição é da ordem dos 101 espera-se

perder apenas 1 a.s. e, por isso, temos con�ança em 15 a.s. do valor da solução.

(Estamos a usar a fórmula: e.m.r.(resultado) ≤ cond(A) × e.m.r.(dados))

1.5. Para este sistema, qual foi o método que obteve melhores resultados?

Resolução:

Os resultados são os mesmos. De facto, se determinar a diferença entre x e xx obtém-se:

>> xx-x

ans =

0

0

0

0

o que mostra que as operações executadas quer pelo operador \ quer pela função lu

seguida dos 2 sistemas triangulares são as mesmas. O operador \ usa o método de Gauss

com escolha parcial de pivot e a função lu usa a decomposição de Doolitle com escolha

parcial de pivot. Logo o operador \ pode ser usado quando pretendemos resolver um

só sistema com a matriz A. A função lu é mais vantajosa quando pretendemos resolver

vários sistemas com a mesma matriz.

1.6. Calcule as normas 1,2 e ∞ da matriz do sistema. Porque obteve resultados diferentes?

Resolução:

norm(A,1)= 7, norm(A)(ou norm(A,2))=5.7016, norm(A,Inf)= 6.

Uma função pode ser designada por norma desde que respeite as respectivas condições

(v. pp. 100 da sebenta). Há uma in�nidade de possibilidades, sendo as normas indica-

das (v pp.101) as mais usadas para vectores e matrizes fornecendo valores diferentes.

Geralmente é o valor da componente maior que mais contribui para o valor da norma.

Este facto, faz com que por vezes as normas só digam o que acontece às componentes

maiores em valor absoluto.

2. Para resolver um dado problema de Engenharia é necessário determinar x tal que Ax = b.Sabese que

A =

1 0 0 0

2.9444 −6.2222 2.9444 00 2.9444 −6.2222 2.94440 0 0 1

e que b = (0.0185, 0.1111, 0.2222, 0.1481)T .

2.1. Faça

A=[1 0 0 0;2.9444 -6.2222 2.9444 0; 0 2.9444 -6.2222 2.9444 ;0 0 0 1];

b=[0.0185;0.1111;0.2222;0.1481];

xx=A\b

c=cond(A)

Resolução:

Qual a solução que obteve? x= 0.018500000000000, 0.009230809998256, 0.038739487152272,

0.148100000000000

Atendendo à condição de A = 13.1008 temos con�ança em 15 a.s. desta solução. Por-

que cond(A) ≈ 101, esperamos perder um a.s. dos cerca de 16 a.s. com que o Matlab

trabalha.

Nota: A solução é apresentada com 15 c.d. para cada componente e não com 15 a.s.

pois, tal como referimos anteriormente, quando se trabalha com normas só �camos

praticamente com a a informação da componente maior. Como o erro relativo é dado

pela norma do erro absoluto a dividir pela norma do vector x, apesar de não sabermos

qual a componente que tem maior erro, sabemos que a componente maior de x é a

última. Assim, na pior das hipóteses, é nessa componente que temos o maior erro

absoluto que nos permite ter con�ança em 15 a.s.. Como a ordem de grandeza dessa

componente é das décimas, os 15 a.s. são 15 c.d. e por essa razão só podemos garantir

15 c.d. para as restantes componentes.

2.2. Faça

format short

[xj,k,rerr,C,Y]=jacobi(A,b,[0;0;0;0],1.e-3,20,1);

para resolver Ax=b pelo método de Jacobi, usando para ponto x0 o vector nulo e usando

para critério de paragem‖xi−x(i−1)‖∞‖xi‖∞ < 1.e− 3 ou, no máximo, 20 de iterações.

O programa fornece a solução em xj, obtida em k iterações. Fornece também rerr=‖xk−x(k−1)‖∞‖xk‖∞ , C a matriz da fórmula de recorrência de Jacobi e a tabela Y com todos

os iterados. Na 1a coluna x0 e na última linha∥∥xi − x(i−1)

∥∥∞.

Resolução:E obteve-se o quadro:

iter # 0 1 2 3 4 5 6 7

x1 = 0.000000 0.018500 0.018500 0.018500 0.018500 0.018500 0.018500 0.018500

x2 = 0.000000 -0.017855 -0.026000 0.003165 0.001342 0.007873 0.007464 0.008927

x3 = 0.000000 -0.035711 0.025922 0.022068 0.035869 0.035006 0.038097 0.037904

x4 = 0.000000 0.148100 0.148100 0.148100 0.148100 0.148100 0.148100 0.148100

||xN-xV|| 0.0 0.148100 0.061633 0.029165 0.013801 0.006531 0.003090 0.001462

iter # 8 9 10 11

x1 = 0.018500 0.018500 0.018500 0.018500

x2 = 0.008835 0.009163 0.009142 0.009216

x3 = 0.038596 0.038552 0.038707 0.038698

x4 = 0.148100 0.148100 0.148100 0.148100

||xN-xV|| 0.000692 0.000327 0.000155 0.000073

� O método é convergente. Parou ao �m de 11 iterações pelo critério de paragem

rerr =

∥∥x(11) − x(10)∥∥∥∥x(11)∥∥ ≤ 10−3 que garante que da 10a para a 11a iteração houve

uma estabilização, na maior das componentes, de cerca de 3 a.s..

Nota: Leiam os comentários da função jacobi para identi�carem o signi�cado dos

parâmetros.

� Existe uma relação entre o e.m.a. cometido e∥∥x(n) − x(n−1)

∥∥ dada por ∥∥x− x(n)∥∥ ≈

K∥∥x(n) − x(n−1)

∥∥. Para obter essa relação faça:

∗ De�na um "script"com as instruções

for i=2:k

v(i-1)=norm(Y(1:4,k+1)-Y(1:4,i),Inf);

end

KR=v./Y(5,2:k)

Corra o "script", qual é a informação guardada em v e KR?

Resolução:

v guarda a norma ∞ da diferença entre a última iteração e as restantes ite-

rações. Quando a última iteração tiver pelo menos mais 1 a.s. do que a

outra iteração a aproximar a solução pretendia x então v também nos dá ∼=∥∥x− x(n)∥∥∞. Como na última linha de Y estão os valores de

∥∥x(n) − x(n−1)∥∥∞,

KR vai fornecer a constante de proporcionalidade K da relação anterior. Como

nas primeiras iterações aquela relação pode ainda não ser válida e nas últimas

o valor de v não aproxima bem o erro∥∥x− x(n)

∥∥∞, convém olhar para os

valores de KR intermédios.

Nota: se olharem para a última linha da tabela (ou Y ) vêem que são precisas

um pouco mais do que 3 iterações para a solução ganhar 1 a.s.. Deste modo,

v só aproxima bem o erro até à 7a iteração.

i 1 2 3 4 5 6 7 8 9 10

KR 0.5024 0.5714 0.5702 0.5705 0.5652 0.5667 0.5430 0.5496 0.4436 0.4732

v 7.44e-2 3.52e-2 1.66e-2 7.87e-3 3.69e-3 1.75e-3 7.94e-4 3.80e-4 1.45e-4 7.33e-5

∗ O valor de K = 0.58 e∥∥x− x(n)

∥∥ ≈0.58 ∥∥x(n) − x(n−1)∥∥.

� Faça

Erro=K*Y(5,k+1)

O erro na última iteração é 4.2532e− 005 < 0.5× 10−4.

� Qual a solução que obteve (apresente-a de modo conveniente)?

x= ( 0.0185(0), 0.0092(2), 0.0387(0), 0.1481(0)).

� Veri�que se os resultados estão correctos

Faça

xx,xj,dif=xx-xj

dif=

1.0e-004 *

-0.0000

0.1525

0.4192

0

Conclusão: o valor estimado (0.43 × 10−4) para o erro, neste caso, além de ser

bastante próximo do erro cometido (0.42×10−4), é também um limite superior do

erro, mas isso não podemos garantir sempre. (Usou-se a norma ∞).

2.3. Faça

format short

[xs,k,rerr,Y]=seidel(A,b,[0;0;0;0],1.e-3,20,1);

para resolver Ax=b pelo método de Gauss-Seidel, usando para ponto x0 o vector nulo e

usando para critério de paragem‖xi−x(i−1)‖∞‖xi‖∞ < 1.e−3 ou, no máximo, 20 de iterações.

O programa fornece a solução em xj, obtida em k iterações. Fornece também rerr=‖xk−x(k−1)‖∞‖xk‖∞ e a tabela Y com todos os iterados. Na 1a coluna x0 e na última linha∥∥xi − x(i−1)

∥∥∞.

Resolução:E obteve-se o quadro:

iter # 0 1 2 3 4 5 6 7

x1 = 0.000000 0.018500 0.018500 0.018500 0.018500 0.018500 0.018500 0.018500

x2 = 0.000000 -0.009101 -0.028038 0.000885 0.007362 0.008812 0.009137 0.009210

x3 = 0.000000 -0.040018 0.021104 0.034790 0.037855 0.038541 0.038695 0.038730

x4 = 0.000000 0.148100 0.148100 0.148100 0.148100 0.148100 0.148100 0.148100

||xN-xV|| 0.0 0.148100 0.061121 0.028923 0.006477 0.001450 0.000325 0.000073

O método é convergente. Parou ao �m de 7 iterações pelo critério de paragem

rerr =

∥∥x(7) − x(6)∥∥∥∥x(7)∥∥ ≤ 10−3 que garante que da 6a para a 7a iteração houve uma

estabilização, na maior das componentes, de cerca de 3 a.s..

Nota: Leiam os comentários da função seidel para identi�carem o signi�cado dos

parâmetros.

* Existe uma relação entre o e.m.a. cometido e∥∥x(n) − x(n−1)

∥∥ dada por∥∥x− x(n)∥∥ ≈ K

∥∥x(n) − x(n−1)∥∥. Para obter essa relação faça:

* De�na um "script"com as instruções

for i=2:k

v(i-1)=norm(Y(1:4,k+1)-Y(1:4,i),Inf);

end

KR=v./Y(5,2:k)Corra o "script", qual é a informação guardada em v e KR?

Resolução:

v guarda a norma ∞ da diferença entre a última iteração e as restantes iterações.

Quando a última iteração tiver pelo menos mais 1 a.s. do que a outra iteração a

aproximar a solução pretendia x então v também nos dá ∼=∥∥x− x(n)

∥∥∞. Como na

última linha de Y estão os valores de∥∥x(n) − x(n−1)

∥∥∞, KR vai fornecer a constante

de proporcionalidade K da relação anterior. Como nas primeiras iterações aquela

relação pode ainda não ser válida e nas últimas o valor de v não aproxima bem o erro∥∥x− x(n)∥∥∞, convém olhar para os valores de KR intermédios.

Nota: se olharem para a última linha da tabela (ou Y ) vêem que, nas últimas iterações,

são precisas um pouco mais do que 3 iterações para a solução ganhar 2 a.s.. Deste

modo, v só aproxima bem o erro até à 5a iteração. Nas primeiras iterações parece que

o método ainda não estabilizou.i 1 2 3 4 5 6

KR 0.5317 0.6094 0.2878 0.2853 0.2741 0.2239

v 7.87e-002 3.72e-002 8.32e-003 1.85e-003 3.97e-004 7.27e-005

- O valor de K = 0.3 e∥∥x− x(n)

∥∥ ≈0.3 ∥∥x(n) − x(n−1)∥∥.

- Faça

Erro=K*Y(5,k+1)

O erro na última iteração é 2.1817e− 005 < 0.5× 10−4.

- Qual a solução que obteve (apresente-a de modo conveniente)?

x= ( 0.0185(0), 0.0092(1), 0.0387(3), 0.1481(0))

- Veri�que se os resultados estão correctos

Faça

xx,xs,dif=xx-xs

dif = 1.0e-004 * (-0.0000, 0.2098, 0.0993, 0)

Conclusão: o valor estimado (0.22×10−4) para o erro, neste caso, além de ser bastante

próximo do erro cometido (0.21 × 10−4), é também um limite superior do erro, mas

isso não podemos garantir sempre. (Usou-se a norma ∞).