Trab02.Labo.Raiz

Prévia do material em texto

CEFET-RJ 
Cálculo Numérico (Laboratório) Profª Natalia Silveira 
 
T2 – Raízes Reais 
 
Os códigos fontes apresentados a seguir calculam as raízes da função 
)x4sen()x(sen)x(f 2  para   ,x através do Método da Bisseção (bissecao.for) e 
Método de Newton-Rapson (newton.for). Em ambos os programas, o domínio é dividido em 32 
partes. 
No programa bissecao.for é verificada a existência de raiz em cada subdivisão, 
através do critério 0)b(f).a(f  e em newton.for encontram-se as raízes estimando-se 
valores iniciais 0x . 
A partir dos modelos abaixo, determine as raízes de uma função polinomial de grau 
qualquer pelos Métodos da Bisseção e de Newton-Raphson e compare os resultados. 
O grau e os coeficientes do polinômio serão lidos em um arquivo de entrada e as 
raízes gravadas em arquivo de saída. 
A defesa deste trabalho será realizada em ..../..../......., e constará de um Teste 
aplicado ao grupo formado por dois alunos. 
 
Arquivo Fonte: bissecao.for 
 
C ******** Método da Bisseção ******** 
 implicit double precision(a-h,o-z) 
c Arquivo de saída: contém as raízes em cada intervalo 
 OPEN (2,FILE='bis.2',STATUS='UNKNOWN') 
 tole_x =1.d-15 
 zero=1.d-15 
 pi = 4.d0*atan(1.d0) 
 xa=-pi 
 xb=pi 
 deltax=xb-xa 
 ndiv=32 
 h=deltax/ndiv 
 write(2,5) 
 do i=1,ndiv 
 diferenca=1.d0 
 xb = xa+h 
 fa=fx(xa) 
 fb=fx(xb) 
 produto=fa*fb 
 xa_aux=xa 
 xb_aux=xb 
 if(produto.le.zero) then ! Tem raiz neste intervalo 
 do while(diferenca.gt.tole_x) 
 xmed=(xa_aux+xb_aux)/2.d0 
 produto=fx(xa_aux)*fx(xmed) 
 if (produto.le.zero) then ! A raiz está na 1ª metade 
 xb_aux=xmed 
 else ! A raiz está na 2ª metade 
 xa_aux=xmed 
 end if 
 diferenca=abs(xa_aux-xb_aux) 
 end do 
 write(2,10)i,xa,fa,xb,fb,xmed 
 else ! Não tem raiz neste intervalo 
 write(2,15)i,xa,fa,xb,fb 
 end if 
 xa=xb 
 enddo 
 5 format(2x,'i',9x,'a',12x,'f(a)',11x,'b',12x,'f(b)',9x,'raiz') 
 10 format(i3,1x,5f14.8) 
 15 format(i3,1x,4f14.8,10x,'-') 
 STOP 
 end 
****************************************************************** 
 function fx(x) 
 implicit double precision(a-h,o-z) 
 fx=sin(x)**2+sin(4.d0*x) 
 return 
 end 
Arquivo de saída de dados: bis.2 
 
 i a f(a) b f(b) raiz 
 1 -3.14159265 0.00000000 -2.94524311 0.74516701 -3.14159265 
 2 -2.94524311 0.74516701 -2.74889357 1.14644661 - 
 3 -2.74889357 1.14644661 -2.55254403 1.01576507 - 
 4 -2.55254403 1.01576507 -2.35619449 0.50000000 - 
 5 -2.35619449 0.50000000 -2.15984495 -0.01576507 -2.16794290 
 6 -2.15984495 -0.01576507 -1.96349541 -0.14644661 - 
 7 -1.96349541 -0.14644661 -1.76714587 0.25483299 -1.86146340 
 8 -1.76714587 0.25483299 -1.57079633 1.00000000 - 
 9 -1.57079633 1.00000000 -1.37444679 1.66904655 - 
 10 -1.37444679 1.66904655 -1.17809725 1.85355339 - 
 11 -1.17809725 1.85355339 -0.98174770 1.39844850 - 
 12 -0.98174770 1.39844850 -0.78539816 0.50000000 - 
 13 -0.78539816 0.50000000 -0.58904862 -0.39844850 -0.68298271 
 14 -0.58904862 -0.39844850 -0.39269908 -0.85355339 - 
 15 -0.39269908 -0.85355339 -0.19634954 -0.66904655 - 
 16 -0.19634954 -0.66904655 0.00000000 0.00000000 -0.00000001 
 17 0.00000000 0.00000000 0.19634954 0.74516701 - 
 18 0.19634954 0.74516701 0.39269908 1.14644661 - 
 19 0.39269908 1.14644661 0.58904862 1.01576507 - 
 20 0.58904862 1.01576507 0.78539816 0.50000000 - 
 21 0.78539816 0.50000000 0.98174770 -0.01576507 0.97364975 
 22 0.98174770 -0.01576507 1.17809725 -0.14644661 - 
 23 1.17809725 -0.14644661 1.37444679 0.25483299 1.28012925 
 24 1.37444679 0.25483299 1.57079633 1.00000000 - 
 25 1.57079633 1.00000000 1.76714587 1.66904655 - 
 26 1.76714587 1.66904655 1.96349541 1.85355339 - 
 27 1.96349541 1.85355339 2.15984495 1.39844850 - 
 28 2.15984495 1.39844850 2.35619449 0.50000000 - 
 29 2.35619449 0.50000000 2.55254403 -0.39844850 2.45860995 
 30 2.55254403 -0.39844850 2.74889357 -0.85355339 - 
 31 2.74889357 -0.85355339 2.94524311 -0.66904655 - 
 32 2.94524311 -0.66904655 3.14159265 0.00000000 3.14159265 
 
Arquivo Fonte: newton.for 
 
c Este programa calcula as ra¡zes de uma funçao pelo método de Newton-Raphson 
 implicit double precision(a-h,o-z) 
 open(6,file='new.2',status='unknown') 
 tole_x=1.d-15 
 ndiv=32 ! Numero de divisoes (ndiv) 
 pi=4.d0*atan(1.d0) 
 xa = -pi ! Intervalo xa e xb 
 xb = pi 
 deltax=xb-xa 
 h=deltaX/dfloat(ndiv) 
 write(6,5) ! Imprimir cabeçalho 
 do i=1,ndiv 
 diferenca=1.d0 
 x0=xa+dfloat(i-1)*h 
 xinicial=x0 
 do while(diferenca>tole_x) 
 fa=fx(x0) 
 dfa=dfx(x0) 
 x1=x0-fa/dfa 
 diferenca=abs(x1-x0) 
 x0=x1 
 enddo 
 write(6,10)i,xinicial,x0,fa 
 enddo 
5 format('divisao',5x,'x_inicial',13x,'raiz',20x,'f(raiz)') 
10 format(i5,f16.8,5x,f16.8,5x,e25.15) 
 stop 
 end 
c ****************************************************************** 
 function fx(x) 
 implicit double precision(a-h,o-z) 
 fx=sin(x)**2+sin(4.d0*x) 
 return 
 end 
c ****************************************************************** 
 function dfx(x) 
 implicit double precision(a-h,o-z) 
 dfx=2.d0*sin(x)*cos(x)+4.d0*cos(4.d0*x) 
 return 
 end 
 
Arquivo de saída de dados: newton.2 
 
divisao x_inicial raiz f(raiz) 
 1 -3.14159265 -3.14159265 0.489842541528951E-15 
 2 -2.94524311 -3.14159265 0.489842541528951E-15 
 3 -2.74889357 -17.56942666 0.206215253206743E-14 
 4 -2.55254403 -3.82457535 -0.347215745738483E-16 
 5 -2.35619449 -2.16794289 -0.741594285980085E-16 
 6 -2.15984495 -2.16794289 -0.741594285980085E-16 
 7 -1.96349541 -1.86146339 0.714489231667947E-16 
 8 -1.76714587 -1.86146339 0.714489231667947E-16 
 9 -1.57079633 -1.86146339 0.714489231667947E-16 
 10 -1.37444679 -2.16794289 -0.741594285980085E-16 
 11 -1.17809725 1.28012926 -0.123002736468480E-15 
 12 -0.98174770 -0.68298270 0.726415455565288E-17 
 13 -0.78539816 -0.68298270 0.726415455565288E-17 
 14 -0.58904862 -0.68298270 0.726415455565288E-17 
 15 -0.39269908 -1.86146339 0.714489231667947E-16 
 16 -0.19634954 0.00000000 0.189326617253043E-27 
 17 0.00000000 0.00000000 0.000000000000000E+00 
 18 0.19634954 0.00000000 0.315544362088405E-29 
 19 0.39269908 -14.42783401 0.164565626750712E-14 
 20 0.58904862 -0.68298270 0.726415455565288E-17 
 21 0.78539816 0.97364976 -0.407660016854550E-16 
 22 0.98174770 0.97364976 -0.407660016854550E-16 
 23 1.17809725 1.28012926-0.123002736468480E-15 
 24 1.37444679 1.28012926 -0.123002736468480E-15 
 25 1.57079633 1.28012926 0.451570204840213E-15 
 26 1.76714587 0.97364976 -0.407660016854550E-16 
 27 1.96349541 4.42172192 0.351281503885303E-16 
 28 2.15984495 2.45860995 -0.932007292522852E-15 
 29 2.35619449 2.45860995 -0.932007292522852E-15 
 30 2.55254403 2.45860995 0.114147515224705E-14 
 31 2.74889357 1.28012926 0.451570204840213E-15 
 32 2.94524311 3.14159265 -0.489842541528951E-15 
 
Gráfico da função )x4sen()x(sen)x(f 2  obtido pelo Excel: 
 
-1.0
-0.5
0.0
0.5
1.0
1.5
2.0
-3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5
f(x)
 
Gráfico 1 –   ,x . 
 
 
Algoritmo – Função fx (Função Polinomial): 
 Função fx(x,a,ngrau) 
 fx=0 
 Faça i=1,ngrau+1 
 fx = fx + )1i(x).i(a  
 Fim faça 
 Fim da Função fx 
 
Algoritmo – Função dfx (Derivada da Função Polinomial): 
 Função dfx(x,a,ngrau) 
 fx=0 
 Faça i=1,ngrau+1 
 dfx = dfx + )2i(x.)i(a.)1i(  
 Fim faça 
 Fim da Função dfx