Resolução de Sistemas de Equações Lineares com Gauss e Cholesky, Trabalhos de Métodos Computacionais

Baixe Resolução de Sistemas de Equações Lineares com Gauss e Cholesky e outras Trabalhos em PDF para Métodos Computacionais, somente na Docsity! UNIVERSIDADE DE PERNAMBUCO 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 (turma BE) Recife, 24 de outubro de 2008. 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 obje�vo 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 mul�plicamos 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,0277 6,456 0 -3,4126 -5,0906 -0,9438 -4,4718 0 -3,0189 5,8526 -1,1199 -0,1076 0 -13,957 -2,3142 5,0738 -19,0504 Passo 2: Definição dos mul�plicadores 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 0 10,3563 -0,2849 3,8486 0 0 18,5140 8,9354 -1,2076 Passo 3: Definição dos mul�plicadores 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 par�r 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 par�r 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) 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 prin�('\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 prin�('\n'); solucaoinicial=input('Digite a solução inicial da seguinte forma: [x1 x2 ... xn]:'); prin�('\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 prin�('\nMatriz expandida para a etapa %1.0f : \n',etapa); sistemanovo=sistemaacumulado; 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: i) Ly=b ii) Uf=y E obtemos as seguintes equações itera�vas: i) 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) ii) 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 end end prin�('\nDigite os termos independentes da seguinte forma: [b1 b2 ... bn] \n'); termosindependentes=input('Termos independentes:'); prin�('\nDigite a solução inicial da seguinte forma: [y1 y2 ... yn] \n'); solucaoinicialyi=input('S(y) inicial:'); prin�('\nDigite a solução inicial da seguinte forma: [x1 x2 ... xn] \n'); solucaoinicialxi=input('S(x) inicial:'); prin�('\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 //ro�na 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 //ro�na 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 //ro�na de varredura de linhas mult=sistemanovoU(i,etapa)/sistemanovoU(etapa,etapa); sistemaacumuladoL(i,etapa)=mult; //ro�na de escrita de mul�plicadores na matriz L for j=1:neq //ro�nade 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 prin�('\nMatriz triangular superior para a etapa %1.0f : \n',etapa); sistemanovoU=sistemaacumuladoU; disp(sistemanovoU) end prin�('\nMatriz triangular inferior para a etapa %1.0f : \n',etapa); sistemanovoL=sistemaacumuladoL; disp(sistemanovoL) novasolucaoinicialyi=solucaoinicialyi; novasolucaoinicialxi=solucaoinicialxi; for cit=1:iteramax //ro�na 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)*nov asolucaoinicialxi(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)*novas olucaoinicialyi(1,varreLclcy); end acumultermoLyb=novasolucaoinicialyi(1,varretermoylinha)+(termind- acumultermoLyb); //aqui também foi adotado sinal + II. + Equações Itera�vas: i) i) Tomando como referência o sinal (+); Para K=0 K=1 K=2 K=3 K=4 K=5 K=6 //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); prin�("\n\nMatriz de Cholesky:\n"); //imprime a matriz de cholesky for(i = 1; i <= n; i++) { for(j = 1; j <= n; j++) { prin�(" 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++) { prin�("\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]; } } prin�("\n\nResultado do sistema Ax=b: "); for(i=0; i<n; i++) { prin�("\n x(%i): %f", i+1, x[i]);