Fazendo um ajuste não linear em dados experimentais - FORTRAN 90
Publicado por Iago Lira (última atualização em 19/04/2017)
[ Hits: 1.943 ]
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 !---------------------------------------------------------------------
trocar permissão, dono e grupo de arquivos ou diretórios
Tranposta da matriz em Haskell
PJEOffice - Baixa automaticamente última versão do CNJ (Conselho Nacional de Justi&cce
Nenhum comentário foi encontrado.
Melhorando o tempo de boot do Fedora e outras distribuições
Como instalar as extensões Dash To Dock e Hide Top Bar no Gnome 45/46
E a guerra contra bots continua
Tradução do artigo do filósofo Gottfried Wilhelm Leibniz sobre o sistema binário
Conheça o firewall OpenGFW, uma implementação do (Great Firewall of China).
Instalando o FreeOffice no LMDE 6
Anki: Remover Tags de Estilo HTML de Todas as Cartas
Colocando uma opção de redimensionamento de imagem no menu de contexto do KDE
Não consigo acessar os modos de desempenho (1)
Criar um script para testar pen drive (5)
Problema com alias usando locate (4)
[Shell Script] Script para desinstalar pacotes desnecessários no OpenSuse
[Shell Script] Script para criar certificados de forma automatizada no OpenVpn
[Shell Script] Conversor de vídeo com opção de legenda
[C/C++] BRT - Bulk Renaming Tool
[Shell Script] Criação de Usuarios , Grupo e instalação do servidor de arquivos samba