Fazendo um ajuste não linear em dados experimentais - FORTRAN 90
Publicado por Iago Lira (última atualização em 19/04/2017)
[ Hits: 2.238 ]
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 !---------------------------------------------------------------------
Programação para sistemas embarcados em Assembly
Zfehwallpaper - wallpaper no Openbox
MoonScript - Agenda telefônica semifuncional em 101 linhas
Nenhum coment�rio foi encontrado.
Conciliando o uso da ZRAM e SWAP em disco na sua máquina
Servidor de Backup com Ubuntu Server 24.04 LTS, RAID e Duplicati (Dell PowerEdge T420)
Visualizar câmeras IP ONVIF no Linux sem necessidade de instalar aplicativos
Realizar overclock no Miyoo Mini (plus ou normal)
Otimização de memória para máquinas modestas
Unbuntu não atualiza o firmware (1)
linux mint reconhece microfone de lapela como fone de ouvido sem micro... (0)
Dúvidas sobre a originalidade de conteúdos online (10)
Erro de interface de Rede no Virt Manager dentro Debian 13 KDE (12)