267 pág.
Pré-visualização | Página 33 de 35
de dados: do 10 i=1,ni. 244 Cap´ıtulo 10 - Algoritmos implementados em FORTRAN 90 Tabela 10.1: Algumas func¸o˜es matema´ticas do FORTRAN 90 Func¸a˜o Definic¸a˜o func¸a˜o definic¸a˜o Exponenciais sqrt(r) √ R exp(r) eR alog(r) Log2(R) alog10 Log10(R) tanh(r) tanh(R) atanh tanh−1(R) sinh(r) senh(R) cosh(r) cosh(R) Trigonome´tricas cos(r) cos(R) acos(r) cos−1(R) sin(r) sen(R) asin(r) sin−1(R) tan(r) tan(R) atan(r) tan−1(R) Complexas abs(r) |(a± bi)| iabs |bi| aimag(r) bi conjg (a± bi) Nume´ricas ifix(r) R(R) →֒ Z(R) float(r) Z(R) →֒ R(R) dble(r) →֒ dupla precisa˜o (R) cmplx(r) →֒ para C(R) amax1(r) maior valor de (R) amin1(r) menor valor de (R) END - Indica fim do programa. ENDIF - Encerra o comando de condic¸a˜o. ELSE - Avaliac¸a˜o da expressa˜o lo´gica quando existe um bloco a ser executado quando a expressa˜o e´ falsa. EXTERNAL - Identifica um nome definido pelo usua´rio como uma subrotina ou func¸a˜o externa. FORMAT - Indica o formato: format(2e16.7) - 2 colunas com elementos reais. FUNCTION - Declarac¸a˜o de uma func¸a˜o externa ao programa. Fund. de Ca´lculo Nume´rico para Engenheiros 245 GOTO - Indica transfereˆncia: goto a - vai para a posic¸a˜o ”a”. IF - THEN Comando de condic¸a˜o, usado para avaliac¸a˜o de expressa˜o lo´gica. Este utiliza-se das seguintes expresso˜es lo´gicas: if(a.lt.b) lt − menor que (a < b) if(a.le.b) le − menor ou igual (a <= b) if(a.eq.b) eq − igual (a == b) if(a.gt.b) gt − maior que (a > b) if(a.ge.b) ge − maior ou igual (a >= b) INTEGER - Converte um argumento de entrada para um inteiro por trunca- mento. OPEN - Estabelece uma conexa˜o entre uma unidade e um arquivo externo e determina as propriedades da conexa˜o. PARAMETER - Atributo que determina o conteu´do da varia´vel declarada, im- pedindo sua modificac¸a˜o durante a execuc¸a˜o. PAUSE - Pa´ra o programa temporariamente. READ - Comando de leitura de dados: read(a,b) - a = identificador; b = formato. REAL - Declara o tipo intr´ınsico dos dados. STOP - Indica que o programa deve parar. WRITE - Comando de impressa˜o de dados: write(a,b) - a = identificador; b = formato. ! - Representa comenta´rio (todo caracter colocado apo´s esse sinal na˜o sera´ lido). Na pro´xima sec¸a˜o apresenta-se algumas implementac¸o˜es utilizadas para a resoluc¸a˜o dos exemplos resolvidos nos cap´ıtulos anteriores. 246 Cap´ıtulo 10 - Algoritmos implementados em FORTRAN 90 10.3 Exemplos de implementac¸o˜es em FORTRAN 90 10.3.1 Me´todo de Newton-Raphson Esta implementac¸a˜o calcula a ra´ız de uma func¸a˜o qualquer. Esta e´ calculada com toleraˆncia (tol = 10−4), utilizando ni iterac¸o˜es. Como exemplo determina-se a raiz da func¸a˜o f(x) = 2x− cos(x). ! ************************************************ ! ****** Me´todo de Newton-Raphson ****** ! ****** Calcula ra´ızes de func¸o˜es ****** ! ************************************************ ! *** Paraˆmetros de entrada: *** ! ni : Nu´mero de iterac¸o˜es; ! tol : Erro ma´ximo permitido; ! x(0) : Estimativa inicial; ! *** Paraˆmetros de Sa´ıda: *** ! x(ni) : ra´ız da func¸a˜o. program newton dimension x(100) tol = 0.0001 ni = 100 x(0) = acos(-1)/8 do i = 1,ni x(i) = x(i-1) - ((2.0*x(i-1) - cos(x(i-1)))/(2.0 + sin(x(i-1)))) write(*,*)i,x(i) if(i.eq.ni) write(*,*) ’Na˜o convergiu.’ Fund. de Ca´lculo Nume´rico para Engenheiros 247 dxr = abs((x(i)-x(i-1))/x(i-1)) if(dxr.lt.tol) goto 40 enddo 40 continue stop end 10.3.2 Gauss-Seidel Este programa calcula um sistema Ax = B de qualquer ordem, neste caso resolve-se 1, 0 −2, 0 −7, 04, 0 2, 9 −4, 6 0, 0 3, 0 7, 5 x1 x2 x3 = −2, 0 8, 1 0, 0 Avalia-se as ra´ızes do sistema com precisa˜o de 4 d´ıgitos significativos exatos. ! ********************************************************** ! ****** Me´todo Iterativo de Gauss-Seidel ****** ! ****** Resolve um sistema Ax=B de ordem 3× 3 ****** ! ********************************************************** ! *** Paraˆmetros de entrada: *** ! a : Matriz dos coeficientes; ! b : Matriz dos termos indepedentes; ! n : Ordem da matriz; ! nk : Nu´mero ma´xino de iterac¸o˜es; ! x : Vetor soluc¸a˜o; ! tol : Ordem de precisa˜o. ! *** Paraˆmetros de Sa´ıda: *** 248 Cap´ıtulo 10 - Algoritmos implementados em FORTRAN 90 ! x : Vetor soluc¸a˜o; ! k : Nu´mero de iterac¸o˜es; ! eps : Erro. Program GaussSeidel integer i,j,n,k,nk real :: a(3,3),x(5),b(3),v(5) real :: r,soma,norma1,norma2,eps,tol tol = 0.0001 n = 3 nk = 1000 a(1,1) = 1. a(1,2) = -2. a(1,3) = -7. a(2,1) = 4. a(2,2) = 2.9 a(2,3) = -4.6 a(3,1) = 0. a(3,2) = 3. a(3,3) = 7.5 b(1) = -2. b(2) = 8.1 b(3) = 0. do i = 1,n r = 1./a(i,i) do j = 1,n if((i.lt.j).or.(i.gt.j)) then Fund. de Ca´lculo Nume´rico para Engenheiros 249 a(i,j) = a(i,j)*r endif enddo b(i) = b(i) * r x(i) = b(i) enddo do k = 1,nk do i = 1,n soma = 0. do j = 1,n if((i.lt.j).or.(i.gt.j)) then soma = soma + a(i,j) * x(j) endif enddo v(i) = x(i) x(i) = b(i) - soma enddo norma1 = 0. norma2 = 0. do i = 1,n if(abs(x(i) - v(i)).gt.norma1) then norma1 = abs(x(i) - v(i)) 250 Cap´ıtulo 10 - Algoritmos implementados em FORTRAN 90 endif if(abs(x(i)).gt.norma2) then norma2 = abs(x(i)) endif enddo eps = norma1/norma2 if((eps.lt.tol).or.(k.ge.nk)) then goto 30 endif enddo ! *** fim iteracao *** 30 continue write(*,*) ’ k ’, ’ x(1) ’, ’ x(2) ’, ’ x(3) ’ , ’ Erro ’ write(*,*) k,x(1),x(2),x(3),eps end program 10.3.3 Regressa˜o linear !******************************************************** !*** PROGRAMA PARA AJUSTE DE CURVA ******* !*** PELO METODO DOS MINIMOS QUADRADOS ****** !******************************************************** ! *** Paraˆmetros de entrada: *** ! x : Valores da abcissa; ! y : Valores da ordenada; Fund. de Ca´lculo Nume´rico para Engenheiros 251 ! a : Coeficientes; ! sx : Somato´rios dos valores em x; ! nj : Nu´mero de elementos; ! *** Paraˆmetros de Sa´ıda: *** ! reta : ax + b; dimension x(6),y(6) nj = 6 open(12, file=’ajuste.dat’, status=’unknown’) ! *** valores tabelados em y *** y(1) = -0.01 y(2) = 0.51 y(3) = 0.82 y(4) = 0.88 y(5) = 0.81 y(6) = 0.49 ! *** valores tabelados em x *** x(1) = -2. x(2) = -1. x(3) = 0. x(4) = 1. x(5) = 2. x(6) = 3. ! **** Inicio regressa˜o linear **** sx = 0. sx2 = 0. sy = 0. sxy = 0. sx2y = 0. 252 Cap´ıtulo 10 - Algoritmos implementados em FORTRAN 90 sx3 = 0. sx4 = 0. do j = 1,nj sx = sx + x(j) sx2 = sx2 + x(j)**2 sy = sy + y(j) sxy = sxy + x(j)*y(j) sx2y = sx2y + x(j)*x(j)*y(j) sx3 = sx3 + x(j)*x(j)*x(j) sx4 = sx4 + x(j)*x(j)*x(j)*x(j) enddo a0 = (sy*sx2 - sx*sxy)/(nj*sx2 - sx*sx) a1 = (nj*sxy - sx*sy)/(nj*sx2 - sx*sx) write(12,*)’ Linear = ’,a1,’x + ’,a0 end 10.3.4 Interpolac¸a˜o de Lagrange ! **************************************************************** ! ****** Me´todo de Interpolac¸a˜o de Lagrange ****** ! ****** Interpola 1 ou mais valores em uma func¸a˜o tabelada ****** ! **************************************************************** ! *** Paraˆmetros de entrada: *** ! tab : Matriz que conte´m os pontos conhecidos de uma func¸a˜o; ! n : Nu´mero de pontos da tabela; Fund. de Ca´lculo Nume´rico para Engenheiros 253 ! npi : Nu´mero de pontos a ser interpolado; ! x : Vetor que conte´m as abscissas dos pontos interpolados; ! *** Paraˆmetros de Sa´ıda: *** ! y : Vetor que conte´m as ordenadas dos pontos interpolados. program lagrange integer i,j,k,n,npi real tab(10,2),x(10),y(10),sum n = 8 npi = 3 open(1,file=’tabela.dat’,status=’old’)
