Metodos_Nadja

Prévia do material em texto

UNIVERSIDADE DE PERNAMUCO
ESCOLA POLITÉCNICA DE PERNAMBUCO
DEPARTAMENTO DO CICLO BÁSICO
MÉTODOS NUMÉRICOS DIRETOS PARA RESOLUÇÃO DE 
SISTEMAS LINEARES
Recife, 24 de outubro de 2008.
UNIVERSIDADE DE PERNAMUCO
ESCOLA POLITÉCNICA DE PERNAMBUCO
DEPARTAMENTO DO CICLO BÁSICO 
MÉTODOS NUMÉRICOS DIRETOS PARA RESOLUÇÃO DE 
SISTEMAS LINEARES
Trabalho acadêmico apresentado à Universidade de Pernambuco – Escola Politécnica de Pernambuco pel
a aluna Nadja Juliana da Silva Lima
, solicitado pel
o
 Professor 
Jornandes Dias
, como requisito parcial na disciplina de 
Métodos Computacionais I
.
Recife, 24 de outubro de 2008.
INTRODUÇÃO
A resolução de sistemas lineares é desafio que surge em diversas áreas do conhecimento. Este trabalho vem com o propósito de apresentar métodos numéricos para a resolução de sistemas lineares com igual número de equações e incógnitas, através da resolução de três questões.
Os métodos aqui utilizados são os classificados como diretos, que são aqueles que fornecem a solução exata do sistema linear após um número finito de operações.
Dentre esses métodos, serão abordados o Método de Eliminação de Gauss, a Fatoração LU e a Fatoração de Cholesky, métodos eficientes. Dessa forma, cada método será utilizado na resolução de uma questão, com a abordagem de toda a sua metodologia.
OBJETIVOS
Este trabalho vem com o objetivo de determinar as soluções das variáveis decisórias de um dado sistema linear com número de equações iguais ao de incógnitas através de um método direto pré-estabelecido, desenvolver um algoritmo para a resolução de sistemas pelo método direto utilizado, calcular os resíduos e mostrar o comportamento das variáveis do sistema ao longo do processo de iteração com a utilização de gráficos e tabelas.
Resolver o sistema de EAL’s utilizando o método de eliminação de Gauss.
1,234 f1(k) + 5,331f2(k) + 0,793 f3(k) - 0,027f4(k) = 6,456
1,333f1(k) + 2,346f2(k ) - 4,234f3(k) - 0,973f4(k) = 2,502 Sf(0) = (1,02 1,12 2,10 1,41)t
0,231f1(k) -2,021 f2(k ) + 6,001f3(k) – 1,125f4(k) = 1,101
3,237f1(k) + 0,027f2(k )- 0,234f3(k) – 5,003f4(k) = -2,115
O processo começa com a construção da matriz aumentada [A,b]:
[A(0),b(0)]= 	a1,1(0) a1,2(0) a1,3(0) a1,4(0) b1(0)		L1
a2,1(0) a2,2(0) a2,3(0) a2,4(0) b2(0) 		L2
a3,1(0) a3,2(0) a3,3(0) a3,4(0) b3(0) 		L3
a4,1(0) a4,2(0) a4,3(0) a4,4(0) b4(0) 		L4
1,234 5,331 0,793 - 0,027 6,456
1,333 2,346 -4,234 -0,973 2,502 
0,231 2,021 6,001 1,125 1,101
3,237 0,027 0,234 –5,003 -2,115
Esse método consiste na aplicação de uma seqüência de operações elementares sobre as linhas da matriz aumentada com o objetivo de transformá-la em um sistema triangular superior.
Passo 1: Inicialmente, vamos eliminar o coeficiente x1 nas linhas 2, 3 e 4, isto é, eliminar os elementos a2,1(0), a3,1(0) e a4,1(0). Para isso multiplicamos as linhas da matriz aumentada pelos fatores
m2,1(1) = = =1,0802 m3,1(1)= = =0,1872 m4,1(1)= = =2,6232
e definimos as linhas:
L1(1) 	 L1(0), L2(1)	 L2(0) - m2,1(1) L1(0) , L3(1)	 L3(0) – m3,1(1) L1(0),
L4(1) 	 L4(0) - m4,1(1) L1(0)
Calculo dos elementos com a1,1(0) 0 para as linhas L2(1), L3(1) e L4(1):
Linha L2(1)
a2,1(1) = a2,1(0) - m2,1(1) a1,1(0) = 1,333 – 1,0802 x 1,234 = 1,333 – 1,3329 = 0
a2,2(1) = a2,2(0) - m2,1(1) a1,2(0) = 2,346 – 1,0802 x 5,331 = -3,4126
a2,3(1) = a2,3(0) - m2,1(1) a1,3(0) = -4,234 – 1,0802 x 0,793 = -5,0906
a2,4(1) = a2,4(0) - m2,1(1) a1,4(0) = -0,973 – 1,0802 x (- 0,027) = -0,9438
b2(1) = b2(0) - m2,1(1) b1(0) = 2,502 – 1,0802 x 6,456 = -4,4718
Linha L3(1)
a3,1(1) = a3,1(0) – m3,1(1) a1,1(0) = 0,231 – 0,1872 x 1,234 = 0
a3,2(1) = a3,2(0) – m3,1(1) a1,2(0) = -2,021 – 0,1872 x 5,331 = -3,0189
a3,3(1) = a3,3(0) – m3,1(1) a1,3(0) = 6,001 – 0,1872 x 0,793 = 5,8526
a3,4(1) = a3,4(0) – m3,1(1) a1,4(0) = -1,125 – 0,1872 x (-0,027) = -1,1199
b3(1) = b3(0) – m3,1(1) b1(0) = 1,101 – 0,1872 x 6,456 = -0,1076
Linha L4(1)
a4,1(1) = a4,1(0) – m4,1(1) a1,1(0) = 3,237 – 2,6232 x 1,234 = 0
a4,2(1) = a4,2(0) – m4,1(1) a1,2(0) = 0,027 – 2,6232 x 5,331 = -13,9572
a4,3(1) = a4,3(0) – m4,1(1) a1,3(0) = -0,234 – 2,6232 x 0,793 = -2,3142
a4,4(1) = a4,4(0) – m4,1(1) a1,4(0) = 5,003 – 2,6232 x (-0,027) = 5,0738
b4(1) = b4(0) – m4,1(1) b1(0) = -2,115 – 2,6232 x 6,456 = -19,0504
Assim, a matriz fica: 
[A(1),b(1)]= 	1,234 5,331 0,793 -0,027 6,456	
0 -3,4126 -5,0906 -0,9438 -4,4718
0 	 - 3,0189 5,8526 -1,1199 -0,1076
0 -13,9572 -2,3142 5,0738 -19,0504
Passo 2: Definição dos multiplicadores e das linhas para o passo k=2.
m3,2(2) = = =0,8847 m4,2(2)= = = 4,0915
L1(2) 	 L1(0), L2(2)	 L2(1) , L3(2)	 L3(1) – m3,2(2) L2(1), L4(2) 	 L4(1) - m4,2(2) L2(1)
Linha L3(2)
a3,2(2) = a3,2(1) – m3,2(2) a2,2(1) = -3,0189 – 0,8847 x (-3,4126) = 0
a3,3(2) = a3,3(1) – m3,2(2) a2,3(1) = 5,8526 – 0,8847 x (-5,0906) = 10,3563
a3,4(2) = a3,4(1) – m3,2(2) a2,4(0) = -1,1199 – 0,8847 x (-0,9438) = -0,2849
b3(2) = b3(1) – m3,2(2) b2(1) = -0,1076 – 0,8847 x (-4,4718) = 3,8486
Linha L4(2)
a4,2(2) = a4,2(1) – m4,2(2) a2,2(1) = -13,9572 – 4,0915 x (-3,4126) = 0
a4,3(2) = a4,3(1) – m4,2(2) a2,3(1) = -2,3142 – 4,0915 x (-5,0906) = 18,5140
a4,4(2) = a4,4(1) – m4,2(2) a2,4(1) = 5,0738 – 4,0915 x (-0,9438) = 8,9354
b4(2) = b4(1) – m4,2(2) b2(1) = -19,0504 – 4,0915 x (-4,4718) = -1,2076
A matriz expandida para esse passo é:
[A(2),b(2)]= 	1,234 5,331 0,793 -0,027 6,456	
0 -3,4126 -5,0906 -0,9438 -4,4718
0 10,3563 -0,2849 3,8486
0 0 18,5140 8,9354 -1,2076
Passo 3: Definição dos multiplicadores e linhas para o passo k=3.
m4,3(3)= = = 1,7877
L1(3) 	 L1(0), L2(3)	 L2(1) , L3(3)	 L3(2) , L4(3) 	 L4(2) - m4,3(3) L3(2)
Linha L4(3)
a4,3(3) = a4,3(2) – m4,3(3) a3,3(2) = 18,5140 – 1,7877x 10,358 = 0
a4,4(3) = a4,4(2) – m4,3(3) a3,4(2) = 8,9354 - 1,7877 x (-0,2849) = 9,4447
b4(3) = b4(2) – m432(3) b3(2) = -1,2076 – 1,7877 x 3,8486 = -8,0877
[A(3),b(3)]= 	1,234 5,331 0,793 -0,027 6,456	
0 -3,4126 -5,0906 -0,9438 -4,4718
0 	 0 10,3563 -0,2849 3,8486
0 0 0 9,4447 -8,0877
A partir dessa matriz obtém-se o sistema de EAL’s triangular superior:
a1,1(0)f1(k) + a1,2(0) f2(k) + a1,3(0) f3(k) + a1,4(0 f4(k)) = b1(0)			
 a2,2(1) f2(k) + a2,3(1) f3(k) + a2,4(1) f4(k) = b2(1)
 a3,3(2) f3(k) + a3,4(2) f4(k) = b3(2) 		
a4,4(3) f4(k) = b4(3) 
Resolvendo o sistema:
f4(k) = 
f3(k) = 
f2(k) = 
f1(k) = 
Integração computacional
 A partir de um estado inicial (EI) após um período de tempo finito obtém-se um estado final previsível (EFP). Sendo assim:
EIn(k) = fn(k)
EFPn(k+1) = fn(k+1)
f1(k+1) = f1(k) + 
f2(k+1) = f2(k) + 
f3(k+1) = f3(k) + 
f4(k+1) = f4(k) + 
Resolvendo
Para k=0:
f1(0+1) = 0,0946
f2(0+1) = -1,0922
f3(0+1) =2,5104
f4(0+1) = 0,5837
Sf(1) = ( 0,0946 -1,0922 2,5104 0,5837)t
Para k=1:
f1(1+1) = 8,4442
f2(1+1) = -3,6880
f3(1+1) = 2,8981
f4(1+1) = -0,2726
Sf(2) = ( 8,4442 -3,6880 2,8981 -0,2726)t
Para k=2:
f1(2+1) = 27,7401
f2(2+1) = -6,6254
f3(2+1)= 3,2622
f4(2+1) = -1,1289
Sf(3) = (27,7401 -6,6254 3,2622 -1,1289)t
Para k=3:
f1(3+1) = 59,4732
f2(3+1) = -9,8691
f3(3+1) = 3,6028
f4(2+1) = -1,9763
Sf(4) = (59,4732 -9,8691 3,6028 -1,9763)t
Para k=4:
f1(4+1) = 105,7912
f2(4+1) = -13,3865
f3(4+1) = 4,0288
f4(4+1) = -2,8326
Sf(5) = (105,7912 -13,3865 4,0288 -2,8326)t
Para k=5:
f1(5+1) = 166,2030
f2(5+1) = -17,3025
f3(5+1) = 4,3225
f4(5+1) = -3,6889
Sf(6) = (1,66,2030 -17,3025 4,3225 -3,6889)t
Para k=6:
f1(6+1) = 243,3248
f2(6+1) =-21,4198
f3(6+1) = 4,6226
f4(6+1) =-4,5452
Sf(7) = (243,3248 -21,4198 4,6226 -4,5452)t
Para k=7:
f1(7+1) = 338,0157
f2(7+1) = -25,7479
f3(7+1) = 4,892
f4(7+1) = -5,4015
Sf(8) = (338,0157 -25,7479 4,8692 -5,4015)t
Para k=8:
f1(8+1) = 451,2337
f2(8+1) = -15,6802
f3(8+1) = 5,0922
f4(8+1) = -4,5452
Sf(9) = (451,2337 -15,6802 5,0922 -4,5452)t
Para k=9:
f1(9+1) = 520,8094 
f2(9+1) = -20,7089
f3(9+1) = 5,3388
f4(9+1) = -3,6889
Sf(10) = (520,8094 -20,7089 5,3388 -3,6889)t
Tabela dos resultados de f
	k\fn
	f1
	f2
	f3
	f4
	0
	1,02
	1,12
	2,1
	1,41
	1
	0,0946
	-1,0922
	2,5104
	0,5837
	2
	8,4442
	-3,688
	2,8981
	-0,2726
	3
	27,7401
	-6,6254
	3,2622
	-1,1289
	4
	59,4732
	-9,8691
	3,6028
	-1,9763
	5
	105,7912
	-13,3865
	4,0288
	-2,8326
	6
	166,203
	-17,3025
	4,3225
	-3,6889
	7
	243,3248
	-21,4198
	4,6226
	-4,5452
	8
	338,0157
	-25,7479
	4,8692
	-5,4015
	9
	451,2337
	-15,6802
	5,0922
	-4,5452
	10
	520,8094
	-20,7089
	5,3388
	-3,6889
Gráfico dos resultados de f
Cálculo dos resíduos:
|
|
|
|
Tabela do calculo dos resíduos
	k\fn
	f1
	f2
	f3
	f4
	1
	0,9254
	2,2122
	0,4604
	0,8263
	2
	8,3496
	2,5958
	0,3877
	0,8563
	3
	19,2959
	2,9374
	0,3641
	0,8563
	4
	31,7331
	3,2437
	0,3406
	1,7037
	5
	46,318
	3,5174
	0,426
	0,8563
	6
	60,4118
	3,916
	0,2937
	0,8563
	7
	77,1218
	4,1173
	0,3001
	0,8563
	8
	94,6909
	4,3281
	0,2466
	0,8563
	9
	113,218
	10,0677
	0,223
	0,8563
	10
	69,5757
	5,0287
	0,2466
	0,8563
 
Gráfico dos resíduos
Algoritmo (Gauss):
//resolução de sistema de equações pelo método de eliminação de Gauss.
clear
clc
disp('Resolução de sistema de equações pelo método de eliminação de Gauss',);
neq=input('Digite o número de Equações Algébricas Lineares:');
for i=1:neq
 printf('\nDigite os coeficientes da equação %1.0f da seguinte forma: [a1 a2 ... an bn] \n',i);
 Equacao=input('Equação:');
 if i>=2
 sistema=[sistema;Equacao];
 else
 sistema=Equacao;
 end
end
printf('\n');
solucaoinicial=input('Digite a solução inicial da seguinte forma: [x1 x2 ... xn]:');
printf('\n');
iteramax=input('Digite o número de iterações:');
disp('Término da alimentação de dados.')
disp('Aguarde...')
tamanhovar=length(sistema);
conf=tamanhovar/neq;
if conf==neq+1
 disp(sistema,'Matriz expandida:');
 sistemanovo=sistema;
 for etapa=1:neq-1
 pivo=sistemanovo(etapa,etapa);
 for b=1:etapa
 if b==1
 sistemaacumulado=sistemanovo(1,:);
 else
 sistemaacumulado=[sistemaacumulado;sistemanovo(b,:)];
 end
 end
 for i=etapa+1:neq
 mult=sistemanovo(i,etapa)/sistemanovo(etapa,etapa);
 for j=1:neq+1
 if j<etapa+1
 novalinha=[0];
 else
 novalinha=[novalinha,(sistemanovo(i,j)-(mult*(sistemanovo(etapa,j))))];
 end
 end
 tamanhonovalinha=length(novalinha);
 for a=tamanhonovalinha:neq+1
 if tamanhonovalinha<neq+1
 novalinha=[0,novalinha];
 tamanhonovalinha=length(novalinha);
 else
 end
 end
 sistemaacumulado=[sistemaacumulado;novalinha];
 end
 printf('\nMatriz expandida para a etapa %1.0f : \n',etapa);
 sistemanovo=sistemaacumulado;
 disp(sistemanovo)
 end
gsolucaoinicial=solucaoinicial;
for cit=1:iteramax
 for termox=1:neq
 solucinicvar=solucaoinicial(1,termox);
 pivonovo=sistemanovo(termox,termox);
 termoindep=sistemanovo(termox,conf);
 acumulatermo=0;
 for termvar=termox+1:neq
 acumulatermo=acumulatermo+(sistemanovo(termox,termvar)*solucaoinicial(1,termvar));
 end
 acumulatermo=solucinicvar+(1/pivonovo)*(termoindep-acumulatermo);
 if termox==1
 solucinicvar1=[acumulatermo];
 else
 solucinicvar1=[solucinicvar1,acumulatermo];
 end
 end
 if cit==1
 iteracaoacumul=[solucaoinicial;solucinicvar1];
 else
 iteracaoacumul=[iteracaoacumul;solucinicvar1];
 end
solucaoinicial=solucinicvar1;
end
disp('Resultados das iterações:')
disp(iteracaoacumul)
eixox=1:iteramax+1;
for nsol=1:neq
 plot(eixox,iteracaoacumul(:,nsol))
end
else
 disp('Matriz não quadrada')
end
Aplicar o método da fatoração LU, admitindo 4 casas decimais nos resultados.
A matriz expandida é igual a: 
[A,b](0)= ; Sf(0) = (0,02 0,37 1,03 2,05)t ; Sy(0) = (0 0 0 0)t
O objetivo do método é fatorar a matriz dos coeficientes A, em um produto entre duas matrizes. Uma delas L, é uma matriz triangular inferior, com todos os elementos da diagonal principal iguais a 1. A outra U, é uma matriz triangular superior com os elementos da diagonal principal diferentes de zero. A matriz pode ser escrita como: 
A=LU
Podemos obter f resolvendo Ly=b resolvendo para y, e Uf=y resolvendo para f.
Aplica-se o método de eliminação de Gauss:
m2,1(1) = = =1,0425 m3,1(1)= = = 12,9792 m4,1(1)= = = 10,3258
e definimos as linhas:
L1(1) 	 L1(0), L2(1)	 L2(0) - m2,1(1) L1(0) , L3(1)	 L3(0) – m3,1(1) L1(0),
L4(1) 	 L4(0) - m4,1(1) L1(0)
Calculo dos elementos com a1,1(0) 0 para as linhas L2(1), L3(1) e L4(1):
Linha L2(1)
a2,1(1) = a2,1(0) - m2,1(1) a1,1(0) = 0
a2,2(1) = a2,2(0) - m2,1(1) a1,2(0) = 1,6807
a2,3(1) = a2,3(0) - m2,1(1) a1,3(0) = -5,0473
a2,4(1) = a2,4(0) - m2,1(1) a1,4(0) = 31,9147
Linha L3(1)
a3,1(1) = a3,1(0) – m3,1(1) a1,1(0) = 0
a3,2(1) = a3,2(0) – m3,1(1) a1,2(0) = -11,2619
a3,3(1) = a3,3(0) – m3,1(1) a1,3(0) = -35,5767
a3,4(1) = a3,4(0) – m3,1(1) a1,4(0) = 157,7152
Linha L4(1)
a4,1(1) = a4,1(0) – m4,1(1) a1,1(0) = 0
a4,2(1) = a4,2(0) – m4,1(1) a1,2(0) = 6,2985
a4,3(1) = a4,3(0) – m4,1(1) a1,3(0) = -37,7547
a4,4(1) = a4,4(0) – m4,1(1) a1,4(0) = 133,4081
[A,b](1)= ;
m3,2(2) = =-6,7011 m4,2(2)= = 3,7478
L1(2) 	 L1(0), L2(2)	 L2(1) , L3(2)	 L3(1) – m3,2(2) L2(1), L4(2) 	 L4(1) - m4,2(2) L2(1)
Linha L3(2)
a3,2(2) = a3,2(1) – m3,2(2) a2,2(1) =0
a3,3(2) = a3,3(1) – m3,2(2) a2,3(1) = -69,3992
a3,4(2) = a3,4(1) – m3,2(2) a2,4(0) = 371,5788
Linha L4(2)
a4,2(2) = a4,2(1) – m4,2(2) a2,2(1) = 0
a4,3(2) = a4,3(1) – m4,2(2) a2,3(1) = -18,8384
a4,4(2) = a4,4(1) – m4,2(2) a2,4(1) = 13,7979
[A,b](2)= ;
m4,3(3)= = 0,2714
L1(3) 	 L1(0), L2(3)	 L2(1) , L3(3)	 L3(2) , L4(3) 	 L4(2) - m4,3(3) L3(2)
Linha L4(3)
a4,3(3) = a4,3(2) – m4,3(3) a3,3(2) = 0
a4,4(3) = a4,4(2) – m4,3(3) a3,4(2) = -87,0486
[A,b](3)= ;
Decompondo nas matrizes Le U obtemos:
L= ; U= ;
Resolvemos os sistemas formados pelas equações:
Ly=b
Uf=y
E obtemos as seguintes equações iterativas:
y1(k+1) = y1(k) + 0,31
y2(k+1) = y2(k) + [0,57 – 1,0425y1(k)]
y3(k+1) = y3(k) + [-1,79 – 12,9792y1(k) + 6,7011 y2(k)]
y4(k+1) = y4 (k) + [2,3322 – 10,3258y1(k) – 3,7475y2(k) -0,2714 y3(k)
f1(k+1) = f1(k) + 
f2(k+1) = f2(k) + 
f3(k+1) = f3(k) + 
f4(k+1) = f4(k) + 
Para k=0:
y1(0+1) = 0,31
y2(0+1) = 0,57
y3(0+1) = -1,79
y4(0+1) = 2,3322
f1(0+1) = 10,0435
f2(0+1) = -35,4641
f3(0+1) =12,0062
f4(0+1) = 2,05
Para k=1:
y1(1+1) = 0,62
y2(1+1) = 0,8168
y3(1+1) = -3,7840
y4(1+1) = -0,1869
f1(1+1) = 32,2131
f2(1+1) = -37,9965
f3(1+1) =23,0082
f4(1+1) = 2,0232
Para k=2:
y1(2+1) = 0,93
y2(2+1) = 0,7405
y3(2+1) = -8,1476
y4(2+1) = -6,2907
f1(2+1) = 40,2434
f2(2+1)= -6,8391
f3(2+1) =33,8954
f4(2+1) = 2,0253
Para k=3:
y1(3+1) = 1,24
y2(3+1) = 0,3410
y3(3+1) = -17,0461
y4(3+1) = -14,1252
f1(3+1) = 8,0723
f2(3+1) = 56,9343
f3(3+1) =44,8567
f4(3+1) = 2,0976
Para k=4:
y1(4+1) = 1,55
y2(4+1) = -17,7688
y3(4+1) = -32,6452
y4(4+1) = -21,2486
f1(4+1) = -89,5597
f2(4+1) = 152,0149
f3(4+1) = 56,1021
f4(4+1) = 2,2599
Para k=5:
y1(5+1) = 1,86
y2(5+1) = -18,8147
y3(5+1) = -54,5530
y4(5+1) = 40,5271
f1(5+1) = -277,0799
f2(5+1) = 267,0094
f3(5+1) =68,6725
f4(5+1) = 2,5040
Para k=6:
y1(6+1) = 2,17
y2(6+1) = -20,1838
y3(6+1) = -206,5635
y4(6+1) = 109,0881
f1(6+1) = -571,5025
f2(6+1) = 414,4965
f3(6+1) =85,8656
f4(6+1) = 2,3843
Tabela dos resultados de f
	k\fn
	f1
	f2
	f3
	f4
	0
	0,02
	0,37
	1,03
	2,05
	1
	10,0435
	-35,4641
	12,0062
	2,05
	2
	32,2131
	-37,9965
	23,0082
	2,0232
	3
	40,2434
	-6,8391
	33,8954
	2,0253
	4
	8,0723
	56,9343
	44,8567
	2,0976
	5
	-89,5597
	152,0149
	56,1021
	2,2599
	6
	-277,0799
	267,0094
	68,6725
	2,504
	7
	-571,5025
	414,4965
	85,8656
	2,3843
Gráfico dos resultados de f
Cálculo dos resíduos:
|
|
|
|
Tabela dos resíduos
	k\fn
	f1
	f2
	f3
	f4
	1
	10,0235
	35,8341
	10,9762
	0
	2
	22,1696
	2,5324
	11,002
	0,0268
	3
	8,0303
	31,1574
	10,8872
	0,0021
	4
	32,1711
	63,7734
	10,9613
	0,0723
	5
	97,632
	95,0806
	11,2454
	0,1623
	6
	187,5202
	114,9945
	12,5704
	0,2441
	7
	294,4221
	147,4871
	17,1931
	0,1197
Gráfico dos resíduos
Algoritmo (L.U):
//resolução de sistema de equações pelo método da Fatoração LU.
clear
clc
disp('Resolução de sistema de EALs pelo método da Fatoração LU',);
printf('\nDigite o número de Equações Algébricas Lineares \n');
neq=input('Número de equações:');
for i=1:neq
 printf('\nDigite os coeficientes da equação %1.0f da seguinte forma: [a1 a2 ... an] \n',i);
 Equacao=input('Equação:');
 if i>=2
 sistema=[sistema;Equacao];
 else
 sistema=Equacao;
 end
end
printf('\nDigite os termos independentes da seguinte forma: [b1 b2 ... bn] \n');
termosindependentes=input('Termos independentes:');
printf('\nDigite a solução inicial da seguinte forma: [y1 y2 ... yn] \n');
solucaoinicialyi=input('S(y) inicial:');
printf('\nDigite a solução inicial da seguinte forma: [x1 x2 ... xn] \n');
solucaoinicialxi=input('S(x) inicial:');
printf('\nDigite o número de iterações \n');
iteramax=input('Iterações:');
disp('Término da alimentação de dados.')
disp('Aguarde...')
tamanhovar=length(sistema);
conf=tamanhovar/neq;
if conf==neq
 disp(sistema,'Matriz de coeficientes:');
 sistemanovoU=sistema;
 for linhasL=1:neq //rotina que gera matriz diagonal
 for colunasL=1:neq
 if linhasL==colunasL
 sistemaacumuladoL(linhasL,colunasL)=1;
 else
 sistemaacumuladoL(linhasL,colunasL)=0;
 end
 end
 end
 for etapa=1:neq-1
 pivo=sistemanovoU(etapa,etapa);
 for b=1:etapa //rotina de composição do novo sist linha/linha
 if b==1
 sistemaacumuladoU=sistemanovoU(1,:);
 else
 sistemaacumuladoU=[sistemaacumuladoU;sistemanovoU(b,:)];
 end
 end
 for i=etapa+1:neq //rotina de varredura de linhas
 mult=sistemanovoU(i,etapa)/sistemanovoU(etapa,etapa); 
 sistemaacumuladoL(i,etapa)=mult; //rotina de escrita de multiplicadores na matriz L 
 for j=1:neq //rotinade varredura de colunas
 if j<etapa+1
 novalinha=[0];
 else
 novalinha=[novalinha,(sistemanovoU(i,j)-(mult*(sistemanovoU(etapa,j))))];
 end
 end
 tamanhonovalinha=length(novalinha);
 for a=tamanhonovalinha:neq
 if tamanhonovalinha<neq
 novalinha=[0,novalinha];
 tamanhonovalinha=length(novalinha);
 else
 end
 end
 sistemaacumuladoU=[sistemaacumuladoU;novalinha];
 end
 printf('\nMatriz triangular superior para a etapa %1.0f : \n',etapa);
 sistemanovoU=sistemaacumuladoU;
 disp(sistemanovoU) 
 end
 printf('\nMatriz triangular inferior para a etapa %1.0f : \n',etapa);
 sistemanovoL=sistemaacumuladoL;
 disp(sistemanovoL)
 novasolucaoinicialyi=solucaoinicialyi;
 novasolucaoinicialxi=solucaoinicialxi;
 for cit=1:iteramax
 //rotina para cálculo de xn
 for varretermoxlinha=1:neq
 pivonovox=sistemanovoU(varretermoxlinha,varretermoxlinha);
 acumultermoUxb=0;
 if varretermoxlinha<neq
 for varreUclcx=varretermoxlinha+1:neq //varredura de colunas
 acumultermoUxb=acumultermoUxb+sistemanovoU(varretermoxlinha,varreUclcx)*novasolucaoinicialxi(1,varreUclcx);
 end
 acumultermoUxb=novasolucaoinicialyi(1,varretermoxlinha)-acumultermoUxb;
 else
 acumultermoUxb=novasolucaoinicialyi(1,varretermoxlinha);
 end
 acumultermoUxb=novasolucaoinicialxi(1,varretermoxlinha)+(acumultermoUxb/pivonovox); //aqui também foi adotado sinal +
 if varretermoxlinha==1
 solucxnovo=acumultermoUxb;
 else
 solucxnovo=[solucxnovo acumultermoUxb];
 end
 end
 novasolucaoinicialxi=solucxnovo;
 solucaoinicialxi=[solucaoinicialxi;novasolucaoinicialxi];
 for varretermoylinha=1:neq
 if varretermoylinha==1
 solucynovo=novasolucaoinicialyi(1,1)+termosindependentes(1,1); //aqui foi adotado sinal +
 else 
 termind=termosindependentes(1,varretermoylinha);
 acumultermoLyb=0; 
 for varreLclcy=1:varretermoylinha-1 //varredura de colunas
 acumultermoLyb=acumultermoLyb+sistemanovoL(varretermoylinha,varreLclcy)*novasolucaoinicialyi(1,varreLclcy);
 end
 acumultermoLyb=novasolucaoinicialyi(1,varretermoylinha)+(termind-acumultermoLyb); //aqui também foi adotado sinal +
 solucynovo=[solucynovo acumultermoLyb];
 end
 end 
 novasolucaoinicialyi=solucynovo;
 solucaoinicialyi=[solucaoinicialyi;novasolucaoinicialyi];
 end
 disp('Resultados das iterações:')
 disp(solucaoinicialyi,'Iterações de y')
 disp(solucaoinicialxi,'Iterações de x')
 eixox=1:iteramax+1;
 for nsol=1:neq
 plot(eixox,solucaoinicialxi(:,nsol))
 end
else
 disp('Matriz não quadrada')
end
Resolver o sistema utilizando a técnica de Cholesky, admitindo 4 casas decimais
 ; 
; 
Como A = , logo, a matriz A é simétrica e positiva. Sendo assim, aplica-se a fatoração de Cholesky. (A=G.)
*
-Cálculo do fator G por colunas:
-Coluna 1
-Coluna 2
-Coluna 3
-Coluna 4
Desta forma, o fator G e o fator transposto são dados por:
; 
-Obtidos os elementos para a matriz G, a resolução do sistema de EAL’s prossegue da seguinte forma:
 =3,71
 =1,45
 =7,27 
=1,37
+ 
Equações Iterativas:
Tomando como referência o sinal (+);
Para K=0
K=1
K=2
K=3
K=4
K=5
K=6
Tabela dos resultados de f
	K\f
	f1
	f2
	f3
	f4
	0
	0,57
	0,91
	0,11
	0,22
	1
	-2, 0666
	11, 7355
	0, 1137
	0,22
	2
	-20, 2446
	103, 4382
	0, 1604
	0, 4633
	3
	-288, 2983
	-314, 0748
	1, 0980
	0, 9714
	4
	669, 1786
	-1752, 5101
	-2, 4544
	2, 1898
	5
	3158, 4234
	-5250, 5997
	-22, 1525
	1, 8596
	6
	18450, 6828
	-12994, 1607
	-75, 9456
	-8, 1345
	7
	56195, 0329
	-28935, 4724
	-188, 1743
	-44, 9067
Gráfico dos resultados de fCálculo dos resíduos:
|
|
|
|
Tabela dos Resíduos:
	K\f
	f1
	f2
	f3
	f4
	1
	2,6366
	10,8255
	0,0037
	0
	2
	18,178
	91,7027
	11,9781
	0,2433
	3
	268,0537
	417,513
	0,9376
	0,5081
	4
	957,4769
	1438,465
	3,5524
	1,2182
	5
	2489,245
	3498,09
	19,6981
	0,33
	6
	15292,26
	7743,561
	53,7931
	9,9941
	7
	37744,35
	15941,31
	112,2287
	36,7722
Gráfico dos Resíduos
Algoritmo (Cholesky):
 //implementação da faroração de Cholesky
 #include<stdio.h>
 #include<stdlib.h>
 #include<conio.h>
 #include<math.h>
 
 //alocação de memória estática
 #define MAX 50
 
 main()
 {
 //declaração de variáveis
 float a[MAX][MAX], b[MAX], L[MAX][MAX];
 long double x[MAX], y[MAX],s[MAX], final[MAX][MAX];
 float soma , soma1, soma2, soma3, soma4 ;
 int i = 0, j = 0, k = 0, n, it, coef;
 
 FILE *saida;
 saida=fopen("CHOLESKY.dat","w+");
 
 //tamanho da matriz de entrada
 n = 4;
 
 //Coeficientes da matrix
 a[0][0]=0.2345;
 a [0][1]= 0.6793;
 a [0][2]= 0.0;
 a [0][3]= 0.0;
 a [1][0]= 0.6793;
 a [1][1]= 1,961;
 a [1][2]= 1.7723;
 a [1][3]= 0.0;
 a [2][0]= 0.0;
 a [2][1]= -1.77;
 a [2][2]= 5.1961;
 a [2][3]= (-2.8284);
 a [3][0]= 0.0;
 a [3][1]= 0.0;
 a [3][2]= (-2.8284);
 a [3][3]= (-5.785);
 
 //primeiro elemento da matriz(a11) 
 L[0][0] = sqrt(a[0][0]);
 for(j = 1; j < n; j++)
 {
 L[j][0] = a[j][0]/L[0][0]; //os outros elementos da coluna 1
 }
 for(i = 1; i < n -1; i++)
 {
 soma = 0.0;
 for(k = 0; k <= i -1; k++)
 {
 soma += pow(L[i][k],2); 
 }
 L[i][i] = sqrt(a[i][i] - soma); //elementos da diagonal principal
 for(j = i + 1; j < n; j++)
 {
 soma1 = 0.0;
 for(k = 0; k <= i-1 ; k++)
 {
 soma1 += L[j][k]*L[i][k];
 }
 L[j][i] = (1/L[i][i])*(a[j][i] - soma1);
 } 
 }
 for(k = 0; k <= i-1 ; k++)//definindo os outros elementos
 {
 soma = 0.0;
 soma2 += pow (L[n-1][k], 2);
 }
 L[n-1][n-1] = pow(a[n-1][n-1] - soma2, 0.5);
 
 printf("\n\nMatriz de Cholesky:\n"); //imprime a matriz de cholesky
 for(i = 1; i <= n; i++)
 {
 for(j = 1; j <= n; j++)
 {
 printf(" G( %d, %d ) = %f\n", i, j, L[i][j]);
 }
 }
 
 //Termos idependentes
 b[0] = 3.71;
 b[1] = 1.45;
 b[2] = 7.27;
 b[3] = 1.37;
 k = 0;
 y[k] = b[k]/L[k][k];// calculando a matriz y
 for(i = k+1 ; i < n ; i++)
 {
 soma3 = 0.0;
 for(j = 0; j < i; j++)
 {
 soma3 += L[i][j]*y[j];
 y[i] = (b[i] - soma3)/L[i][i];
 }
 }
 
 for (i=0; i<n; i++)
 {
 printf("\n y(%i): %f", i+1, y[i]); //imprimindo a matriz y
 }
 x[n-1] = y[n-1]/L[n-1][n-1];//calculando x
 
 for(i=n-2;i>-1;i--)
 {
 soma4= 0.0;
 
 for(j=i+1;j<n;j++)
 {
 soma4 += L[j][i]*x[j];
 x[i] = (y[i] - soma4)/ L[i][i]; 
 } 
 }
 printf("\n\nResultado do sistema Ax=b: ");
 for(i=0; i<n; i++)
 {
 printf("\n x(%i): %f", i+1, x[i]);
 }
 
 //Solução inicial dada no problema
 s[0] = 0.57;
 s[1] = 0.91;
 s[2] = 0.11;
 s[3] = 0.22;
 
 //Quantidade de iterações
 it = 40;
 
 //Carregar no vetor final a solução de k = 0 e k =1
 
 for(i = 0;i<n; i++)
 {
 final[0][i] = s[i]; //para k = 0
 final[1][i] = x[i]; //para k = 1
 } 
 
 for(i = 2; i <it + 2; i++)
 {
 for( j = 0; j< n; j++)
 {
 final[i][j] = final[i - 2][j] + final[i - 1][j];
 // x(k + 1) = x(k-1) + x(k), como k vai ser a partir de 2 então
 //x(i=2) = x(2 - 2 ) + x (2-1) 
 }
 }
 for (i = 0; i < it +2; i++) //imprimindo todas as iterações
 {
 printf("\n\n");
 for( j = 0; j< n; j++)
 {
 fprintf(saida, "\t %.4e\t" ,final[i][j]);
 printf( " %.4e", final[i][j]);
 }
 fprintf(saida,"\n"); 
 printf("\n\n");
 } 
 
 
 system("pause"); 
}