Universidade Federal de Minas Gerais
Departamento de Ciência da Computação
Análise Numérica
Trabalho Prático 1
Decomposição LU e Integração Numérica
Professor: Frederico Ferreira Campos Filho
Alunos: Bruno Xavier da Silva ([email protected])
Diego de Moura Duarte ([email protected])
Belo Horizonte
22 de Abril de 2008
Decomposição LU:
O objetivo da primeira parte do trabalho era a implementação, na forma de subprogramas, dos
algoritmos de decomposição LU do livro Algoritmos Numéricos 2ªedição.
O programa lê um arquivo de nome matriz.dat onde a primeira linha é a ordem N da matriz, as
linhas subsequentes devem conter as linhas da matriz de ordem N em seguida as linhas devem conter as
linhas do vetor independente que também possui comprimento N. A matriz a ser decomposta deve ter
no máximo ordem 100 se for necessário uma decomposição de uma matriz de ordem superior é
necessário apenas modificar o formato de leitura do arquivo matriz.dat o que é bastante simples.
Conforme pedido na especificação o programa após ler os parâmetros de entrada realiza a
decomposição PA=LU, resolve o sistema Ly=Pb, resolve o sistema Ux=y, calcula o determinante da
matriz A, acha a norma-2 do vetor resíduo que é dada por r=||b-Ax|| 2 e grava os resultados no arquivo
matriz.saida.
Utilizando o programa para resolver o sistema que foi passado na especificação, as
configurações do arquivo de entrada e a saída gerada para tal sistema são mostradas logo abaixo. Temos
primeiro a ordem da matriz, depois a matriz e abaixo deste, o vetor independente. O NC(número de
chamada do primeiro aluno da dupla) é igual a quatro.
Matriz.dat
5
0. 1. 3. -2. 4.
-2. 4. 5. -1. 4.
5. 1. 1. -7. -2.
8. -2. 9. 1. 2.
7. -5. 3. 1. 0.
-1.
0.4
2.
0.8
-5.
Matriz.saida
Decomposicao LU
Matriz=
0.0000
-2.0000
5.0000
8.0000
7.0000
1.0000 3.0000 -2.0000 4.0000
4.0000 5.0000 -1.0000 4.0000
1.0000 1.0000 -7.0000 -2.0000
-2.0000 9.0000 1.0000 2.0000
-5.0000 3.0000 1.0000 0.0000
Vetor Independente=
-1.0000
0.4000
2.0000
0.8000
-5.0000
Matriz apos decomposicao LU
8.0000 -2.0000 9.0000 1.0000
-0.2500 3.5000 7.2500 -0.7500
0.6250 0.6429 -9.2857 -7.1429
0.0000 0.2857 -0.1000 -2.5000
0.8750 -0.9286 -0.2000 0.8000
Solucao=
13.7141
14.7859
-11.5397
7.5500
8.4833
Determinante=
312.0000
Residuo=
0.000000000000000000000000
0.000000000000001443289932
0.000000000000003552713679
2.0000
4.5000
-6.1429
2.1000
-0.4800
0.000000000000020650148258
0.000000000000018651746814
Norma Segunda=
0.000000000000028089520156
Na impressão dos valores foram utilizados quatro casas decimais exceto para o vetor resíduo e
para a norma segunda com o objetivo de mostra que apesar de eles serem muito pequenos ele não são
iguais a zero o que seria a conclusão obtida devido ao erro de arredondamento do computador.
Código em Fortran da Decomposição LU
program main
integer i, j, N
real*8 Det
real*8, dimension(:,:), allocatable :: Matriz
real*8, dimension(:), allocatable :: Pivot
real*8, dimension(:), allocatable :: Solucao
real*8, dimension(:), allocatable :: Vetor_Ind
real*8, dimension(:), allocatable :: y
real*8, dimension(:), allocatable :: Residuo
real*8, dimension(:,:), allocatable :: A
real*8 Norma
open(25, file='Matriz.dat', status='old')
open(10, file='Matriz.saida', status='old')
write(10, *) "\t\tDecomposicao LU\n"
read(25, *) N
allocate(Matriz(N, N))
allocate(Pivot(N))
allocate(Solucao(N))
allocate(Vetor_Ind(N))
allocate(y(N))
allocate(Residuo(N))
allocate(A(N,N))
do i=1, n
read(25, *) (Matriz(i,j), j=1, n)
enddo
do i=1, n
read(25, 6) Vetor_Ind(i)
enddo
close(25, status='keep')
do i=1, n
do j=1, n
A(i,j) = Matriz(i,j)
enddo
enddo
write(10,*) "\nMatriz="
do i=1, n
write (10, 6) (Matriz(i,j), j=1, n)
enddo
write(10,*) "\nVetor Independente="
do i=1, n
write (10, 6) Vetor_Ind(i)
enddo
call Decomposicao_LU(N, Matriz, Det, Pivot)
call Subs_Suce_Pivotal(N, Matriz, Vetor_Ind, Pivot, y)
call Subs_Retroativas(N, Matriz, y, Solucao)
call CalResiduo(N, A, Solucao, Vetor_Ind, Residuo)
call Norma2(N, Residuo, Norma)
write(10,*) "\nMatriz apos decomposicao LU"
do i=1, n
write (10, 6) (Matriz(i,j), j=1, n)
enddo
write(10,*) "\nSolucao="
do i=1, n
write (10, 5) Solucao(i)
enddo
write(10,*) "\nDeterminante="
write(10, 5) Det
write(10,*) "\nResiduo="
do i=1, n
write (10, 7) Residuo(i)
enddo
write(10,*) "\nNorma Segunda="
write(10, 7) Norma
close(10, status='keep')
5 format(f10.4)
6 format(100f10.4)
7 format(f30.24)
end
C-----------Decomposição LU-----------C
subroutine Decomposicao_LU(N, A, Det, Pivot)
real*8 Det, Amax, Mult, t, m, r
real*8 A(N, N)
real*8 Pivot(N)
integer N, i, j, k, p
do i=1, n
Pivot(i)=i
enddo
Det=1
do j=1, n-1
p=j
Amax=abs(A(j,j))
do k=j+1, n
if (abs(A(k,j))>Amax) then
Amax=abs(A(k,j))
p=k
endif
enddo
if (p/=j) then
do k=1, n
t=A(j,k)
A(j,k)=A(p,k)
A(p,k)=t
enddo
m=Pivot(j)
Pivot(j)=Pivot(p)
Pivot(p)=m
Det=-Det
endif
Det=Det*A(j,j)
if(abs(A(j,j))/=0) then
r=1/A(j,j)
do i=j+1, n
Mult= A(i,j)*r
A(i,j)=Mult
do k=j+1, n
A(i,k)=A(i,k) - (Mult*A(j,k))
enddo
enddo
endif
enddo
Det=Det*A(n,n)
end
C-----------Substituições Sucessivas Pivotal-----------C
subroutine Subs_Suce_Pivotal(ordem, Matriz, Vetor_Ind, Pivot, y)
integer ordem
real*8 Matriz(ordem,ordem)
real*8 Vetor_Ind(ordem)
real*8 Pivot(ordem)
real*8 y(ordem)
real*8 Soma
integer k, i, j
k= Pivot(1)
y(1)=Vetor_Ind(k)
do i=2, ordem
Soma=0
do j=1, i-1
Soma= Soma + Matriz(i,j) * y(j)
enddo
k=Pivot(i)
y(i)= Vetor_Ind(k) - Soma
enddo
end
C-----------Substituições Retroativas-----------C
subroutine Subs_Retroativas(N, U, d, x)
integer N, i, j
real*8 U(N,N)
real*8 d(N)
real*8 x(N)
real*8 Soma
x(N)=d(N)/U(N,N)
i = N-1
do while(i>=1)
Soma=0
j=i+1
do while(j<=n)
Soma = Soma+U(i,j)*x(j)
j = j+1
enddo
x(i) = (d(i) - Soma)/U(i,i)
i = i-1
enddo
end
C-----------Calcula Resíduo-----------C
subroutine CalResiduo(N, Matriz, x, Vetor_Ind, Residuo)
integer N
real*8 Matriz(N,N)
real*8 x(N)
real*8 Vetor_Ind(N)
real*8 Residuo(N)
real*8 Produto(N)
call Mult(N, Matriz, x, Produto)
call Sub(N, Vetor_Ind, Produto, Residuo)
end
C-----------Multiplica Matriz por Vetor-----------C
subroutine Mult(N, Matriz, vetor, Produto)
integer N, i, j
real*8 temp
real*8 Matriz(N,N)
real*8 vetor(N)
real*8 Produto(N)
do i=1, N
temp = 0.
do j=1, N
temp = temp + Matriz(i,j)*vetor(j)
enddo
Produto(i) = temp
enddo
end
C-----------Subtrai Vetores-----------C
subroutine Sub(N, Vetor1, Vetor2, Subtracao)
integer N, i
real*8 Vetor1(N)
real*8 Vetor2(N)
real*8 Subtracao(N)
do i=1, N
Subtracao(i) = Vetor1(i) - Vetor2(i)
enddo
end
C-----------Calcula Norma-2-----------C
subroutine Norma2(N, Vetor, Norma)
integer N, i
real*8 Vetor(N)
real*8 temp
real*8 Norma
real*8 Somatorio
Somatorio = 0.
do i=1, N
Somatorio = Somatorio + (Vetor(i)*Vetor(i))
enddo
Norma = dsqrt(Somatorio)
end
Integração Numérica:
O objetivo da segunda parte do trabalho era a implementação, na forma de subprogramas, dos
algoritmos de integração numérica do livro Algoritmos Numéricos 2ªedição.
O programa lê um arquivo de nome integral.dat onde a primeira linha é o limite inferior de
integração, a segunda é o limite superior de integração, a terceira é a tolerância e a quarta é o número
máximo de iterações. A integral a ser calculada deve a ter a sua função de integração escrita no código
do programa enquanto os outros parâmetros devem ser modificados no arquivo integral.dat o que é
bastante simples.
Conforme pedido na especificação, o programa após ler os parâmetros de entrada ele chama a
subrotina para integração iterativa, imprime os resultados intermediários, o valor da integral, a menor
diferença relativa obtida e a condição de erro e grava os resultados no arquivo integral.saida.
Utilizando o programa para resolver a integral que foi passada na especificação, as
configurações do arquivo de entrada e a saída gerada para tal integral são mostradas logo abaixo. O
NC(número de chamada do primeiro aluno da dupla) é igual a quatro.
Integral.dat
1.
4.
1e-8
20.
Integral.saida
Integracao
Limite Inferior =
1.0000
Limite Superior =
4.0000
Tolerancia =
0.000000010000000000000000
Iteracoes =
20.0000
Iter
1
2
3
4
5
n
8
13
21
34
55
Integral
6.42596843753310
6.34109130065408
6.34391038565535
6.34398037927947
6.34398035142837
Dif.Relativa
1.338525702512337E-002
4.443765485161463E-004
1.103307701670910E-005
4.390161363289855E-009
Delta =
4.390161363289855E-009
Integral =
6.34398035142837
CondErro =
0
Código em Fortran da Integração Numérica
program main
integer n
integer CondErro
real*8 a, b, Integral, Toler, IterMax, Delta
open(100, file='integral.dat', status='old')
read(100, *) a
read(100, *) b
read(100, *) Toler
read(100, *) IterMax
close(100, status='keep')
open(200, file='integral.saida', status='old')
write(200,*) "\t\tIntegracao\n"
write(200,'(a, f10.4)')"Limite Inferior =\n",a
write(200,'(a, f10.4)')"\nLimite Superior =\n",b
write(200,'(a, f33.24)')"\nTolerancia =\n",Toler
write(200,'(a, f10.4)')"\nIteracoes =\n",IterMax
call GL_Iter(a,b,Toler,IterMax,Integral,Delta,CondErro)
close(200, status='keep')
5 format(f10.4)
end
C-----------Gauss-Legendre-Iterativo-----------C
subroutine GL_Iter(a,b,Toler,IterMax,Integral,Delta,CondErro)
integer CondErro
real*8 a, b, Toler, IterMax
real*8 Integral,Delta,Int2
integer Iter, n1, n2
Iter = 1
n1 = 5
n2 = 8
call Gauss_Legendre(a,b,n2,Int2,CondErro)
write(200, *) "\n\tIter\t
n\tIntegral\t\tDif.Relativa"
write(200, *) Iter, n2, Int2
c
Repita
10 Iter = Iter+1
n = n1 + n2
call Gauss_Legendre(a,b,n,Integral,CondErro)
if(Integral /= 0) then
Delta = abs((Integral-Int2)/Integral)
else
Delta = abs(Integral-Int2)
endif
write(200, *) Iter, n, Integral, Delta
if(Delta <= Toler .or. Iter == IterMax) then
goto 20
endif
Int2 = Integral
n1 = n2
n2 = n
goto 10
20 write(200,*) "\nDelta ="
write(200, *) Delta
write(200,*) "\nIntegral ="
write(200, *) Integral
write(200,*) "\nCondErro ="
write(200, *) CondErro
if(Delta<=Toler)then
CondErro = 0
else
CondErro = 1
endif
end
C-----------PesAbsGL-----------C
subroutine PesAbsGL(n, A, T, CondErro)
integer n, m, i, j
integer CondErro
real*8 A(n)
real*8 T(n)
real*8 pi
real*8 z, z1
real*8 p1, p2, p3, p4, pp
if( n < 1 ) then
CondErro = 1
return
endif
CondErro = 0
pi = 3.14159265358979323846
m = aint(0.5 * (n+1))
do i=1, m
z=cos(pi*(i-0.25)/(n+0.5))
c
REPITA
30
p1 = 1
p2 = 0
do j=1, n
p3 = p2
p2 = p1
p1 = ((2*j-1)*z*p2-(j-1)*p3)/j
enddo
pp = n*(z*p1 - p2)/(z*z - 1)
z1 = z
z = z1-p1/pp
if(abs(z-z1) >= 1e-15)then
goto 30
endif
T(m + 1 - i) = z
A(m + 1 - i) = 2/((1-z*z)*pp*pp)
enddo
end
C-----------Gauss-Legendre-----------C
subroutine Gauss_Legendre(a, b, n, Integral, CondErro)
integer CondErro, n, i, k
real*8 a, b, c, Integral, e1, e2, c1, c2, t, x
real*8, dimension(:), allocatable :: Avet
real*8, dimension(:), allocatable :: Tvet
real*8 temp
allocate(Avet(n))
allocate(Tvet(n))
Integral=0.
call PesAbsGL(n, Avet, Tvet, CondErro)
if (CondErro/=0) then
return
endif
e1=(b-a)/2
e2=(a+b)/2
if (mod(n,2)==0) then
c1=1
c2=0.5
else
c1=0
c2=1
endif
do i=1, n
k=dint(i-0.5 * (n+1)+ Sinal(i-0.5 * (n+c1)) * c2)
temp = k
t=Sinal(temp)*Tvet(abs(k))
x= e1*t + e2
y=F(x)
c=Avet(abs(k))
Integral= Integral + y*c
enddo
Integral=e1*Integral
end
C-----------Sinal-----------C
function Sinal(Valor)
real*8 Valor
if(Valor>0.) then
Sinal=1
endif
if(Valor<0.) then
Sinal=-1
endif
if(Valor == 0) then
Sinal=0
endif
end
C-----------F(x) (Função a ser integrada)-----------C
function F(x)
real*8 x, temp, NC
NC=4.
temp=sin(5*x)*sin(5*x)
temp= NC + temp
F=dsqrt(temp)
end