Esta é uma pré-visualização de arquivo. Entre para ver o arquivo original
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