Fazendo um ajuste não linear em dados experimentais - FORTRAN 90

Publicado por Iago Lira (última atualização em 19/04/2017)

[ Hits: 572 ]

Homepage: NO

Download FitDate.f90




Olá pessoal, como sou amante do Fortran, resolvi criar um programa que faz um ajuste não linear em dados experimentais.

- Eu sei, já existe programas para tais! A grande utilidade é  quando usa-se muitos parâmetros a serem determinados, o que é vantajoso em relação aos demais.

No programa existe a função PLOT, onde nesta dá-se a entrada da função a fazer o ajuste. Exemplo:

FUNCTION PLOT(X,A,Qp)
IMPLICIT NONE
INTEGER, PARAMETER :: qpl = selected_real_kind(15, 16)
INTEGER :: Qp
real(KIND=qpl) :: X, A(Qp), PLOT

PLOT=1.0_qpl+A(1)*exp(-A(2)*X)*cos(A(3)*X)-&
                &A(4)*exp(-A(5)*X)*cos(A(6)*X)-&
                &A(7)*exp(-A(8)*X)*cos(A(9)*X)+&
                &A(10)*exp(A(11)*X)*cos(A(12)*X)-&
                &A(13)*exp(-A(14)*X)*cos(A(15)*X)
      
END FUNCTION PLOT

Percebe-se que a função PLOT têm 16 parâmetros a serem determinados, então percebe-se que é fácil entrar com os valores.

COMPILANDO

No terminal digite:

gfortran Fit.Date.f90 -o FitDate.x -O3

O "-O3" é opcional, pois é um parâmetro de otimização.

EXECUTANDO

Ainda no terminal, digite:

./FitDate.x

Então, aparecerá uma tela pedidos os arquivo que contém os dados a serem analisados, a quantidade de parâmetros e o erro que você quer cometer. Quanto menor o erro, mais demorado.

No programa existe uma variável chamada de "tol" (PARAMETER(tol=0.000000000001)), esta é a precisão do cálculo, então ajuste para suas necessidades.

  



Esconder código-fonte

PROGRAM FIT
IMPLICIT NONE
INTEGER, PARAMETER :: qpl = selected_real_kind(15, 16)
real(KIND=qpl), ALLOCATABLE :: X(:), Y(:), A(:), F(:), B(:)
real(KIND=qpl) :: tol, ai, af, PLOT, ERRO
REAL :: r
PARAMETER(tol=0.000000000001) 
!ai, af -> variacao inicial e final das constantes
INTEGER :: I, nlines, COUNTL, Qp, K, P, PLOTC
CHARACTER*20 :: filename
! LOGICAL :: FLAG

  !ABRE O ARQUIVO COM DADOS dexpERIMENTAIS
  WRITE(*,'(A)', ADVANCE="NO") "DIGITE O NOME DO ARQUIVO: "
  READ*, filename
  !CHAMA A FUNCAO COUNTL
  nlines=COUNTL(filename)
  
  !PERGUNTANDO A QUANTIDADE DE PARAMETROS
  write(*,'(A)', ADVANCE="NO") "DIGITE A QUANTIDADE DE PARAMETROS: "
  READ*, Qp
  
  !TAMANHO DO ERRO
  write(*,'(A)', ADVANCE="NO") "DIGITE O ERRO EM PORCENTAGEM (0<ERRO<100): "
  READ*, ERRO
  
  ALLOCATE(X(nlines), Y(nlines), A(Qp), F(nlines), B(Qp))
  
  OPEN(UNIT=1, FILE=TRIM(ADJUSTL(filename)), STATUS="OLD", ACTION="READ")
   DO I=1,nlines
    READ(1,*) X(I), Y(I)
   ENDDO
  CLOSE(UNIT=1)
   
   
  OPEN(UNIT=2, FILE="parametros.dat")
  A=1
  ai=0; af=1
  
  CALL SYSTEM('clear')
  PRINT*, 'Aguarde um pouco...'
  PRINT*, '---------------------------------------------------------------------'
  
  P=0
  CALL init_random_seed()
  DO I=1,nlines
   DO

    !============CHAMA FUNCAO==============
    F(I)=PLOT(X(I),A,Qp)
    !================================    
        IF (F(I)>Y(I)) af=af-tol
        IF (F(I)<Y(I)) af=af+tol
        
        IF (P.eq.1) THEN
             WRITE(2,*) A
             WRITE(*,*) I, "Experimento=",Y(I), "Fit=",F(I)
             EXIT
        ELSE
!           DO
           DO K=1,Qp
             CALL RANDOM_NUMBER(r)
             !aleatorio entre (x,y)
             A(K)=(r*(af-ai))+ai
           ENDDO
               P=PLOTC(X,Y,A,Qp,nlines,ERRO)
!            IF (P == 1) EXIT   
!           ENDDO
        END IF
   ENDDO
  ENDDO
    PRINT*, '---------------------------------------------------------------------'
    PRINT*, 'Terminado!'
    PRINT*, '---------------------------------------------------------------------'
  CLOSE(UNIT=2)
  
  OPEN(UNIT=1, FILE="fit.dat")
    CALL CALCPAR(Qp,nlines, B)
    print*, B
    DO I=1, nlines
!     XX=Grid*X
      WRITE(1,*) X(I), PLOT(X(I),B,Qp)
    ENDDO
  CLOSE(UNIT=1)
  PRINT*, '---------------------------------------------------------------------'
END PROGRAM FIT
!---------------------------------------------------------------------

!---------------------------------------------------------------------
INTEGER FUNCTION PLOTC(X,Y,A,Qp,nlines,ERRO)
IMPLICIT NONE
INTEGER :: I,Qp, nlines
INTEGER :: P1(nlines)
INTEGER, PARAMETER :: qpl = selected_real_kind(15, 16)
real(KIND=qpl) :: X(nlines), Y(nlines), A(Qp), G(nlines), PLOT
real(KIND=qpl) :: Et, Eb, ERRO, PP(nlines)

P1=0
 DO I=1,nlines
  G(I)=PLOT(X(I),A,Qp)
  PP(I)=ABS(((Y(I)-G(I))/Y(I)))*100 !Erro percentual
  IF(PP(I)>=0.and.PP(I)<=ERRO) P1(I)=1
 ENDDO

  IF(SUM(P1) .eq. NINT(ABS((nlines-((nlines*ERRO)/100))))) THEN 
   PLOTC=1
  ELSE
   PLOTC=0
  ENDIF
END FUNCTION PLOTC
!---------------------------------------------------------------------


!---------------------------------------------------------------------
!ESCREVA A FUNÇÃO AQUI!
FUNCTION PLOT(X,A,Qp)
IMPLICIT NONE
INTEGER, PARAMETER :: qpl = selected_real_kind(15, 16)
INTEGER :: Qp
real(KIND=qpl) :: X, A(Qp), PLOT
 
 PLOT=1.0_qpl+A(1)*exp(-A(2)*X)*cos(A(3)*X)-&
                &A(4)*exp(-A(5)*X)*cos(A(6)*X)-&
                &A(7)*exp(-A(8)*X)*cos(A(9)*X)+&
                &A(10)*exp(A(11)*X)*cos(A(12)*X)-&
                &A(13)*exp(-A(14)*X)*cos(A(15)*X)
      
END FUNCTION PLOT
!---------------------------------------------------------------------

!---------------------------------------------------------------------
!CONTAR NUMERO DE LINHAS NO ARQUIVO
INTEGER FUNCTION COUNTL(filename)
 CHARACTER*20 :: filename

 OPEN(UNIT=1, FILE=TRIM(ADJUSTL(filename)), STATUS="OLD", ACTION="READ")

    COUNTL = 0 
    DO 
        READ(1,*, END=10) 
        COUNTL = COUNTL + 1 
    END DO 
    10 CLOSE (1) 

END FUNCTION COUNTL
!---------------------------------------------------------------------

!---------------------------------------------------------------------
SUBROUTINE init_random_seed()
INTEGER :: i, n, clock
INTEGER, DIMENSION(:), ALLOCATABLE :: seed

CALL RANDOM_SEED(size = n)
ALLOCATE(seed(n))

CALL SYSTEM_CLOCK(COUNT=clock)

seed = clock + 37 * (/ (i - 1, i = 1, n) /)
CALL RANDOM_SEED(PUT = seed)
DEALLOCATE(seed)
END SUBROUTINE 
!---------------------------------------------------------------------

!---------------------------------------------------------------------
SUBROUTINE CALCPAR(Qp,nlines, B)
IMPLICIT NONE
INTEGER, PARAMETER :: qpl = selected_real_kind(15, 16)
INTEGER :: Qp, nlines
real(KIND=qpl) :: A(Qp,nlines), B(Qp)

 OPEN(UNIT=10, FILE="parametros.dat")
  read(10,*) A
 CLOSE(UNIT=10)
 
   B=A(:,nlines-1)
 
END SUBROUTINE CALCPAR
!---------------------------------------------------------------------

Scripts recomendados

Método de Gauss-Seidel em SCILAB

Fortran - Integrais

Crivo de Eratóstenes Simples em XBase (Clipper)

SSH NAO NAVEGA

impressora brother


  

Comentários

Nenhum comentário foi encontrado.


Contribuir com comentário




Patrocínio

Site hospedado pelo provedor HostGator.
Linux banner
Linux banner
Linux banner

Artigos

Dicas

Viva o Android

Tópicos

Top 10 do mês

Scripts