Equações Lineares em MatLab

Funções do Matlab

Resolução de Equações Lineares

‰ Como determinar a solução de um sistema de equações lineares AX = B ?

1º. Introduza a matriz A

» A=[4 1 2; 0 1 0; 8 4 5]
A=
4 1 2
0 1 0
8 4 5

2º. Introduza o vector dos temos independentes b

← não esquecer que é transposto
» b=[1 0 2] '

b=
1
0
2

3º. Resolva o sistema de equações lineares Ax = b (o Matlab implementa o método de

eliminação de Gauss para resolver este sistema)

» x=A\b
x=

0.2500
0
0

4º. Verifique que x é a solução do sistema Ax = b

» A*x
ans =
1
0
2

‰ função LU: factoriza uma matriz

(Esta função encontram-se no Matrix Functions, numerical linear algebra toolbox)

[L,U]=LU(A) – factoriza a matriz A=L*U, onde L – uma matriz triangular

inferior e U – uma matriz triangular superior.

Como determinar a solução de um sistema de equações lineares Ax = B usando o

método de factorização LU?

» A=[4 1 2; 0 1 0; 8 4 5]

A=
4 1 2
0 1 0
8 4 5

» b=[1 0 2]' % introduza o vector b transposto

1 Gladys Castillo Jordan

 b=
1
0
2

» [L,U]=lu(A) % factorize a matriz A usando o comando lu

L=
0.5000 -1.0000 1.0000
0 1.0000 0
1.0000 0 0
U=

8.0000 4.0000 5.0000

0 1.0000 0
0 0 -0.5000

» L*U % verifique que A=L*U

ans =

4 1 2
0 1 0
8 4 5

» y=L\b % resolva Ly=b

y=
2
0
0

» x=U\y % resolva Ax=y

x=
0.2500
0
0
obtém-se a mesma solução do sistema

Utilizando o método de factorização LU com matriz de permutação:

A = [1 2 3; 2 4 7; 3 5 3];

A=
1 2 3
2 4 7
3 5 3

Factorizando a matriz A na forma LU sem matriz de permutação

» [L,U]=lu(A)

L=
0.3333 0.5000 1.0000
0.6667 1.0000 0
1.0000 0 0
U=
3.0000 5.0000 3.0000
0 0.6667 5.0000
0 0 -0.5000

Verificando se L*U = A
ans =
3 5 3
2 4 7
1 2 3

2 Gladys Castillo Jordan

As linhas 1 e 3 estão trocadas

Factorizando a matriz A na forma LU com matriz de permutação

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

L=
1.0000 0 0
0.6667 1.0000 0
0.3333 0.5000 1.0000

U=
3.0000 5.0000 3.0000
0 0.6667 5.0000
0 0 -0.5000

P=
0 0 1
0 1 0
1 0 0

Obtemos a matriz de permutação que corresponde à permutação das linhas 1 e 3.

Resolver o sistema de equaçoes Ax=b usando a factorizaçao LU com a matriz de permutação P

1º: Aplicar a permutação P ao vector b

» b=[1, -1, 0.5]'

b=
1.0000
-1.0000
0.5000

» pb=P*b

pb =
0.5000
-1.0000
1.0000

2º: Resolver o sistema de equações Ly=pb

» y=L\pb

y=

0.5000
-1.3333
1.5000
3º: Resolver o sistema de equações Ux=y

» x=U\y

x=

-31.0000
20.5000
-3.0000

3 Gladys Castillo Jordan

‰ função NORM: determina a norma de uma matriz ou um vector

Para as matrizes:
NORM(A,1)= max(sum(abs((X)))-
determina ||A||1 - o máximo das somas por coluna dos valores absolutos
dos elementos de A
NORM(A,inf) =max(sum(abs((X')))
determina ||A||∞ - o máximo das somas por linha dos valores absolutos
dos elementos de A
Para os vectores:
NORM(X,1) = max(sum(abs((X))). Determina ||X||1 - norma absoluta
NORM(X,inf) = max(abs(X)). Determina ||X||∞ - norma de máximo
NORM(X,-inf) = min(abs(X)). Determina ||X||-∞ - norma de mínimo

For exemplo:

» A=[1 2 4; 4 -1 1; 2 5 2]

A=

1 2 4
4 -1 1
2 5 2

» norm(A,1) % o máximo das somas por coluna dos valores absolutos dos
elementos de A
ans =

» norm(A,inf) % o máximo das somas por linha dos valores absolutos dos
elementos de A
ans =

» b=[11 8 3]'

b=

11
8
3

» norm(b,1) % a soma dos valores absolutos dos elementos de b

ans =

22

» norm(b,inf) % o máximo em valor absoluto dos elementos de b

ans =

11

» norm(b,-inf) % o mínimo em valor absoluto dos elementos de b

ans =

‰ função COND: determina o número de condição de uma matriz

cond(A)= ||A-1|| ||A||

COND(A,P) – devolve o número de condição de X na P-norma, P=1,2, inf:

4 Gladys Castillo Jordan

 For exemplo:

» A=[1 2 4; 4 -1 1; 2 5 2]
A=
1 2 4
4 -1 1
2 5 2

» cond(A,inf)

ans =

4.17391304347826

Exercícios em Matlab

1. Escreva em Matlab a função X=residuos (A,b,X0,eps,maxit) que implementa

o método dos resíduos para melhorar uma aproximação da solução do sistema
de equações lineares Ax=b dada uma solução inicial X0, e tal que ||b - Axk || <
eps. Alem disso um número máximo de iterações, maxit, também é dado, para
terminar o processo iterativo quando não se verificar ||b-Axk || < eps.

function [X,error, k, Z]=residuos(A,b,X0,eps,maxit)

%-------------------------------------------------------------------
% Esta função implementa o método dos residuos para melhorar uma
% aproximação da solução de um sistema de equações lineares Ax=b
% Para executar chamar:
% X = residuo(A,b,X0,eps,maxit)
% [X,error,k, Z] = residuo(A,b,X0,eps,maxit)
% Parâmetros de entrada
% A - a matriz do sistema
% b - o vector com os termos independentes
% X0 - a solução que se pretende melhorar
% eps - a tolerancia para o criterio de paragem ||b-A*X||<eps
% maxit - numero máximo de iterações
% Parâmetros de saída
% X - a aproximação da solução
% error - a norma do residuo
% k - o número de iterações
% Z - Matrix com todas as soluções (uma por linha)
%-------------------------------------------------------------------
echo on;
format long;
X =X0;
Z =X0'; % inicializo Z (armazena a historia)
for k=1:maxit;
r = b - A*X;
if (norm(r,inf)<eps)
k=k-1;
break, end;
e=A\r;
X=X+e;
Z = [Z;X']; % armazeno a historia
end;
error=(norm(r,inf));
end

5 Gladys Castillo Jordan

2. Escreva em Matlab a função M=majorrx(A,b, delta) que determina um
majorante do erro relativo (em norma infinita) da solução x do sistema de
equações lineares Ax=b quando o vector b dos termos independentes é
ligeiramente perturbado segundo um delta dado.
function M=majorrx(A,b,delta);
%--------------------------------------------------------------
% Esta função determina um majorante do erro relativo
%( em norma infinito)da solução x do sistema Ax=b
% quando o vector b dos termos independentes é ligeiramente
% perturbado segundo um delta definido: ||b-b'||<= delta
% O majorante do erro relativo em norma pode ser calculado pela
% fórmula: rx=||delta_x||/||x||<= cond(A)* (||delta_b||/||b||)
% Para executar chamar:
% M = majorrx(A,b,delta)
% Parâmetros de entrada
% A - a matriz A do sistema
% b - o vector com os termos independentes
% delta - o majorante para a perturbação do vector b
% Parâmetros de saída
% M - o majorante do erro relativo da solução para o
% sistema perturbado
%--------------------------------------------------------------
echo on;
format long;

cA=cond(A,inf) % número de condição da matriz (norma inf)

nb=norm(b,inf) % norma inf de b
M=cA / nb * delta % o majorante

6 Gladys Castillo Jordan

3. Escreva em Matlab a função jacobi(A,b, X0, eps, max) que implementa o
método iterativo de Jacobi para determinar a solução aproximada de Ax=b,
dada uma solução inicial X0 tal que ||xk - xk-1 || < eps. Alem disso um número
máximo de iterações, maxit, também é dado, para terminar o processo
iterativo quando não se verificar ||xk - xk-1 || < eps.

function [X,deltaX,k,Z] = jacobi(A,b,X0,eps,maxit)

%-----------------------------------------------------------------------
% Implementa o método iterativo de Jacobi para determinar uma solução
% aproximada de Ax=b:
% X(k)(i)= (b(i)-A(i,1)*X(k-1)(1)-...-A(i,i-1)*X(k-1)(i-1)-
% A(i,i+1)*X(k-1)(i+1)-...-A(n,n)*X(k-1)(n))/ Aii
% Executar
% [X,deltaX,k] = jacobi(A,B,P,eps,maxit)
% [X,deltaX,k,Z] = jacobi(A,B,P,eps,maxit)
% Entrada
% A a matriz A do sistema
% b o vector dos termos independentes
% X0 a solução inicial
% eps tolerância de convergência para critério de paragem
% maxit número máximo de iterações para critério de paragem
% Devolve
% X o vector com a aproximação da solução
% error a norma do vector abs(X(k)-X(k-1))
% k o número de iterações
% Z matriz com todas as soluções (uma por linha)
%-----------------------------------------------------------------------
n = length(b);
Xant = X0; % armazena a solução anterior X(k-1)
X=X0; % inicializo X
Z = X0'; % inicializo Z
for k=1:maxit, % iterar até max veces
for i = 1:n, % para cada equação
X(i) = (b(i) - A(i,[1:i-1,i+1:n])*Xant([1:i-1,i+1:n]) ) / A(i,i);
end
Z = [Z;X']; % armazeno a historia
deltaX = norm(X-Xant,inf);
if (deltaX <eps) break, end
Xant = X;
End

7 Gladys Castillo Jordan

4. Escreva em Matlab a função gseidel(A,b, X0, eps, max) que implementa o
método iterativo de Gauss-Seidel para determinar a solução aproximada de
Ax=b, dada uma solução inicial X0, e tal que ||xk - xk-1 || < eps. Alem disso um
número máximo de iterações, maxit, também é dado, para terminar o processo
iterativo quando não se verificar ||xk - xk-1 || < eps.

function [X,deltaX, k, Z] = gseidel(A,b,X0,eps,maxit)

%--------------------------------------------------------------------------
% Implementa o método iterativo de Gauss-Seidel para determinar uma solução
aproximada de Ax=b
% X(k)(i)= (b(i)-A(i,1)*X(k)(1)-...-A(i,i-1)*X(k)(i-1)
% -A(j,j+1)*X(k-1)(j+1)-...-A(n,n)*X(k-1)(n))/ Aii
% Executar
% [X, deltaX, k] = gseidel(A,B,P,delta,maxit)
% [X, deltaX, k, Z] = dseidel(A,B,P,delta,maxit)
% Entrada
% A a matriz A do sistema
% b o vector dos termos independentes
% X0 a solução inicial
% eps tolerancia de convergencia do criterio de paragem
% maxit número máximo de iterações
% Devolve
% X o vector com a aproximação da solução
% deltaX a norma do vector abs(X(k)-X(k-1))
% k o número de iterações
% Z Matrix com todas as soluções (uma por linha)
%--------------------------------------------------------------------------
n = length(b);
Xant = X0; % armazena a solução anterior X(k-1)
X = X0; % inicializo X
Z = X0'; % inicializo Z
for k=1:maxit, % iterar até max veces
for i = 1:n, % para cada equação
if i==1
X(i) =(b(1) - A(1,2:n)*Xant(2:n))/ A(i,i);
elseif i==n
X(i) = (b(n) - A(n,1:n-1)*X(1:n-1) )/ A(i,i);
else
X(i) = (b(j)-A(i,1:i-1)*X(1:i-1)-A(i,i+1:n)*Xant(i+1:n))/ A(i,i);
end
end
Z = [Z;X']; % armazeno a historia
deltaX = norm(X-Xant,inf);
if (deltaX < eps) break, end
Xant = X;
end

8 Gladys Castillo Jordan

5. Partindo da aproximação inicial X0=[5/4,-1/4]' aproxime a solução do sistema

2x1 + 0.5 x2= 2.5

x1 + 3 x2= 0.5

aplicando o método de Jacobi e o método de Gauss-Seidel com o número de iterações

suficientes por forma a que ||xk - xk-1 || < 0.5 x 10-1

% -------------------------------------------------------------------------
% Este programa partindo da aproximação inicial X0=[5/4,-1/4]
% aproxima a solução do sistema
% 2x1 + 0.5 x2= 2.5
% x1 + 3 x2= 0.5
% aplicando o método de Jacobi e o método de Gauss-Seidel por forma a que
% ||x(k) – x(k-1) || < 0.5 x 10-1
% -------------------------------------------------------------------------
A=[2 0.5; 1 3] % introduzir a matriz A
b=[2.5 0.5]' % introduzir o vector b
X0=[5/4,-1/4]' % introduzir a solução inicial
[X,dX,k,Z]=jacobi(A,b,X0,0.05,2);
m=length(Z);
k=1:m;
disp('');
disp('---------------------Método de Jacobi --------------------');
disp(' Iteração x(1) x(2) ');
disp('--------------------------------------------------------------');
disp([k',Z])
disp('--------------------------------------------------------------');
disp(' ');
disp(['A solução: ','x1= ',num2str(Z(m,1),8),' x2= ',num2str(Z(m,2),8)]);
[X,dX,k,Z]=gseidel(A,b,X0,0.05,2);
m=length(Z);
k=1:m;
disp(' ');
disp('-------------------Método de Gauss-Seidel --------------------');
disp(' Iteração x(1) x(2) ');
disp('--------------------------------------------------------------');
disp([k',Z])
disp('--------------------------------------------------------------');
disp(' ');
disp(['A solução: ','x1= ',num2str(Z(m,1),8),' x2= ',num2str(Z(m,2),8)]);

» ex5
A=

2.00000000000000 0.50000000000000
1.00000000000000 3.00000000000000
b=

2.50000000000000
0.50000000000000
X0 =

1.25000000000000
-0.25000000000000

---------------------Método de Jacobi --------------------

Iteração x(1) x(2)
--------------------------------------------------------------
1.00000000000000 1.25000000000000 -0.25000000000000
2.00000000000000 1.31250000000000 -0.25000000000000
3.00000000000000 1.31250000000000 -0.27083333333333

--------------------------------------------------------------

A solução: x1= 1.3125 x2= -0.27083333

9 Gladys Castillo Jordan

-------------------Método de Gauss-Seidel --------------------
Iteração x(1) x(2)
--------------------------------------------------------------
1.00000000000000 1.25000000000000 -0.25000000000000
2.00000000000000 1.31250000000000 -0.27083333333333
3.00000000000000 1.31770833333333 -0.27256944444444

--------------------------------------------------------------

A solução: x1= 1.3177083 x2= -0.27256944

10 Gladys Castillo Jordan