Calculo Numérico com Fortran – Parte 2

Continuando com a série de artigos sobre Fortran, hoje demonstrarei um algoritmo para calculo da soma de Riemman, que é um metodo para calculo da área abaixo da curva sobre um gráfico. Ela é usada para definir a idéia da Integral Definida.

Sugiro que leiam http://www.uff.br/webmat/Calc1_LivroOnLine/Cap21_Calc1.html

Para implementar o método a primeira etapa e definir o intervalo fechado [a,b] que no programa chamamos de limite_inferior e limite_superior, no caso do código abaixoo limite ficou definido enre [-10,10]

 limite_inferior = -10
 limite_superior = 10

Depois precisamos calcular a distância entre a e b, nesse ponto tive o primeiro erro, pois para calcular a distancia não basta somar o valor de a e b, se somarmos -10 e 10 por exemplo a distancia seria 0. Portanto antes de mais nada devemos calcular o modulo.

 !faz o modulo dos limite
 if (limite_inferior.lt.0) then
 modulo_inferior = -1 * limite_inferior
 else
 modulo_inferior = limite_inferior
 endif

 if (limite_superior.lt.0) then
 modulo_superior = -1 * limite_superior
 else
 modulo_superior = limite_superior;
 endif

A precisão desse método vária de acordo com o numero de divisões para cada intervalo x. Sendo que a solução fica mais precisa quanto mais X+n-X tende a 0. Nesse caso n tende a infinito sendo essa a definição da integral definida. Embaixo segue a implementação dessa idéia, quanto maior for n, mais precisa é a solução.

! numero de interavlos entre o limite inferior e o superior
 n=100
 precisao = (modulo_inferior+modulo_superior)/n

Depois definimos x0 como sendo o limite_inferior e x1 o resultado da soma de x0 com a precisão. A partir dai iniciamos um loop, onde o critério de parada é x1 ser igual ou por questões de arredondamento maior que o limite superior.
Dentro do loop, calculamos y_esquerdo e y_direito o valor de ambos e multiplicado pelo valor da distancia entre x0 e x1 (precisão) ao final do loop temos valores para soma_esquerda e soma_direita que indicam o intervalo em que a solução esta contida.

 DO

 !soma esquerda
 if (x0+precisao.lt.limite_superior) then
 f_esquerdo = f(x0)
 soma_esquerda = soma_esquerda + (f_esquerdo*precisao)
 endif

 !soma direita
 if (x1+precisao.lt.limite_superior) then
 f_direito = f(x1)
 soma_direita = soma_direita + (f_direito*precisao)
 endif

 print *,soma_direita,",",soma_esquerda

 x0 = x1
 x1 = x1 + precisao

 ! Se o x1 alcançou o limite superior pare
 if (x1.GE.limite_superior) then
 exit
 endif

 END DO

Definimos a lei que determina a função no trecho de código abaixo:

 FUNCTION F(X)
 REAL*4, INTENT (IN) :: X
 REAL*4 Y/0/
 !definimos a funçao que desejamos avaliar.
 Y = x**4-10*x**3+130*x
 F=Y
 RETURN
 END FUNCTION F

Segue abaixo o código completo:

PROGRAM RIEMMAN

REAL*4 soma_direita/0/
 REAL*4 soma_esquerda/0/
 REAL*4 n/0/
 REAL*4 precisao/0/
 REAL*4 limite_inferior/0/
 REAL*4 limite_superior/0/
 REAL*4 f_direito/0/
 REAL*4 f_esquerdo/0/
 REAL*4 modulo_inferior/0/
 REAL*4 modulo_superior/0/
 REAL*4 x0/0/
 REAL*4 x1/0/

 PRINT *,"Entre um Byte e Outro - Aprendendo Fortran e Calculo Numerico"
 PRINT *,"Anthony M Collucci - Algoritmo para Calculo da soma de Riemman"
 PRINT *,"=============================================================="

 limite_inferior = -10
 limite_superior = 10

 !faz o modulo dos limite
 if (limite_inferior.lt.0) then
 modulo_inferior = -1 * limite_inferior
 else
 modulo_inferior = limite_inferior
 endif

 if (limite_superior.lt.0) then
 modulo_superior = -1 * limite_superior
 else
 modulo_superior = limite_superior;
 endif

 ! Calcula a distancia entre o limite superior e limite inferior e divide pelo numero de intervalos
! numero de interavlos entre o limite inferior e o superior
 n=100
 precisao = (modulo_inferior+modulo_superior)/n
 ! define o x inicial e o x final
 x0 = limite_inferior
 x1 = limite_inferior + precisao

 print *,"soma_direita",",","soma_esquerda"
 DO

 !soma esquerda
 if (x0+precisao.lt.limite_superior) then
 f_esquerdo = f(x0)
 soma_esquerda = soma_esquerda + (f_esquerdo*precisao)
 endif

 !soma direita
 if (x1+precisao.lt.limite_superior) then
 f_direito = f(x1)
 soma_direita = soma_direita + (f_direito*precisao)
 endif

 print *,soma_direita,",",soma_esquerda

 x0 = x1
 x1 = x1 + precisao

 ! Se o x1 alcançou o limite superior pare
 if (x1.GE.limite_superior) then
 exit
 endif

 END DO
 !Debug
 print *,"O valor da área esta entre",soma_esquerda,"e",soma_direita

 STOP

 CONTAINS

 FUNCTION F(X)
 REAL*4, INTENT (IN) :: X
 REAL*4 Y/0/
 !definimos a funçao que desejamos avaliar.
 Y = x**4-10*x**3+130*x
 F=Y
 RETURN
 END FUNCTION F

END PROGRAM RIEMMAN

Vamos avaliar o programa, primeiro acessamos o cmd, e digitamos gnuplot na console digitamos plot x**4-10*x**3+130*x e será exibido o gráfico abaixo:

Sempre gosto de dar uma olhada no gráfico antes de começar, agora imagine que no intervalo de -10 a 10 exibido acima dividimos ele em 100 partes e traçamos gráficos de barra para y(x0) e y(x1) e somamos as áreas desses intervalos. Se fizessemos isso manualmente teríamos um resultado semelhante à saída abaixo:

Se fosse necessário mais precisão bastaria aumentar n, veja abaixo o valor para n=10000

Perceba que o meu erro já e muito baixo, a diferença entre a estimativa para cima e estimativa por baixo, e de apenas 375 milésimos, no meu caso decidi parar por aqui, se você executar o programa verá que ele se torna mais lento quanto maior for N pois o numero de cálculos cresce exorbitantemente.

Deixe uma resposta

O seu endereço de e-mail não será publicado. Campos obrigatórios são marcados com *

Please type the characters of this captcha image in the input box

Por favor, digite os caracteres desta imagem na caixa de entrada