Fazendo um ajuste não linear em dados experimentais - FORTRAN 90
Publicado por Iago Lira (última atualização em 19/04/2017)
[ Hits: 2.046 ]
Homepage: https://notabug.org/iagolira/
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.
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 !---------------------------------------------------------------------
Verifica se o link caiu e manda aviso por email - MIkrotik v5.*
MoonScript - Agenda telefônica semifuncional em 101 linhas
Software via GPO no Logon de Usuário - SAMBA e AD
Nenhum comentário foi encontrado.
Compartilhando a tela do Computador no Celular via Deskreen
Como Configurar um Túnel SSH Reverso para Acessar Sua Máquina Local a Partir de uma Máquina Remota
Configuração para desligamento automatizado de Computadores em um Ambiente Comercial
Efeito "livro" em arquivos PDF
Como resolver o erro no CUPS: Unable to get list of printer drivers
Flatpak: remover runtimes não usados e pacotes
Mudar o gerenciador de login (GDM para SDDM e vice-versa) - parte 2
Acessar ftp por nome ao invés de ip [RESOLVIDO] (5)
Estou com sede em aprender sobre o nosso querido Linux. (2)
[Python] Automação de scan de vulnerabilidades
[Python] Script para analise de superficie de ataque
[Shell Script] Novo script para redimensionar, rotacionar, converter e espelhar arquivos de imagem
[Shell Script] Iniciador de DOOM (DSDA-DOOM, Doom Retro ou Woof!)
[Shell Script] Script para adicionar bordas às imagens de uma pasta