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