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.

• Download FitDate.f90
• Enviar nova versão

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
!---------------------------------------------------------------------

Zfehwallpaper - wallpaper no Openbox

Raízes - Fortran

Welcome!

Método da Bissecção em SCILAB

impressora brother

Nenhum comentário foi encontrado.