0

试图在 fortran77 中实现单参数最小二乘最小化。这是代码；它可以编译并且似乎可以工作，除了......它陷入了h1 = 1.8E-2和3.5E-2的值之间的无限循环。现在看看，但很可能，我自己不会有太多运气来解决这个问题。欢迎大家帮忙！

      PROGRAM assignment

      ! A program designed to fit experiemental data, using the method
      ! of least squares to minimise the associated chi-squared and
      ! obtain the four control parameters A,B,h1 and h2.
      !*****************************************************************

      IMPLICIT NONE
      INTEGER i
      DOUBLE PRECISION t(17),Ct(17),eCt(17)
      DOUBLE PRECISION h1loop1,h1loop2,deltah,Cs
      DOUBLE PRECISION chisqa,chisqb,dchisq

      OPEN(21, FILE='data.txt', FORM='FORMATTED', STATUS='OLD')
      DO i=1,17
          READ(21,*)t(i),Ct(i),eCt(i)
      END DO
      CLOSE(21)

      !Read in data.txt as three one dimensional arrays.
      !*****************************************************************
      !OPEN(21, FILE='outtest.txt', FORM='FORMATTED', STATUS='NEW')
      !DO i=1,17
      !   WRITE(21,*)t(i),Ct(i),eCt(i)
      !END DO
      !CLOSE(21)
      !
      !Just to check input file is being read correctly.
      !*****************************************************************
      !**********************Minimising Lamda1 (h1)*********************
      deltah= 0.0001
      h1loop2= 0.001
      h1loop1= 0.0 !Use initial value of 0 to calculate start-point chisq
      DO 10
           chisqa= 0.0
          DO 20 i= 1, 17
              Cs= exp(-h1loop1*t(i))
              chisqa= chisqa + ((Ct(i) - Cs)/eCt(i))**2
 20       END DO

           chisqb= 0.0
          DO 30 i= 1, 17
              h1loop2= h1loop2 + deltah
              Cs= exp(-h1loop2*t(i))
              chisqb= chisqb + ((Ct(i) - Cs)/eCt(i))**2
 30       END DO

         !Print the two calculated chisq values to screen.
          WRITE(6,*) 'Chi-squared a=',chisqa,'for Lamda1=',h1loop1
          WRITE(6,*) 'Chi-squared b=',chisqb,'for Lamda1=',h1loop2

          dchisq= chisqa - chisqb
          IF (dchisq.GT.0.0) THEN
              h1loop1= h1loop2
          ELSE
              deltah= deltah - ((deltah*2)/100)
          END IF

          IF (chisqb.LE.6618.681) EXIT
 10   END DO

      WRITE(6,*) 'Chi-squared is', chisqb,' for Lamda1 = ', h1loop2
      END PROGRAM assignment

编辑：再看一遍我决定我不知道是什么搞砸了。应该从中得到 6618.681 的卡方，但它只是停留在 6921.866 和 6920.031 之间。帮助！

4

1 回答 1

1 
do i=1

没有开始循环，对于循环，您还需要指定上限：

do i=1,ub

这就是为什么你会收到关于 doi 没有类型的错误消息，在固定格式中，空格是微不足道的......

编辑：如果你想有一个无限循环，只需完全跳过“i =”声明。当达到某个条件时，您可以使用 exit 语句离开循环：

do
  if (min_reached) EXIT
end do

Edit2：我不知道你为什么坚持 F77 固定格式。这是你的自由格式的程序，有一些地方的修复，看起来很奇怪，没有深入研究细节：

PROGRAM assignment

  ! A program designed to fit experiemental data, using the method
  ! of least squares to minimise the associated chi-squared and
  ! obtain the four control parameters A,B,h1 and h2.
  !*****************************************************************

  IMPLICIT NONE

  integer, parameter :: rk = selected_real_kind(15)
  integer, parameter :: nd = 17

  integer :: i,t0

  real(kind=rk) :: t(nd),t2(nd),Ct(nd),eCt(nd),Ctdiff(nd),c(nd)
  real(kind=rk) :: Aa,Ab,Ba,Bb,h1a,h1b,h2a,h2b,chisqa,chisqb,dchisq
  real(kind=rk) :: deltah,Cs(nd)

  OPEN(21, FILE='data.txt', FORM='FORMATTED', STATUS='OLD')
  DO i=1,nd
     READ(21,*) t(i),Ct(i),eCt(i)
  END DO
  CLOSE(21)

  !Read in data.txt as three one dimensional arrays.
  !*****************************************************************
  !OPEN(21, FILE='outtest.txt', FORM='FORMATTED', STATUS='NEW')
  !DO i=1,17
  !   WRITE(21,*)t(i),Ct(i),eCt(i)
  !END DO
  !CLOSE(21)
  !
  !Just to check input file is being read correctly.
  !*****************************************************************
  !****************************Parameters***************************

  Aa= 0
  Ba= 0
  h1a= 0
  h2a= 0

  !**********************Minimising Lamda1 (h1)*********************
  deltah= 0.001_rk
  h1b= deltah

  minloop: DO

    chisqa= 0
    DO i= 1,nd
       Cs(i)= exp(-h1a*t(i))!*Aa !+ Ba*exp(-h2a*t(i))
       Ctdiff(i)= Ct(i) - Cs(i)
       c(i)= Ctdiff(i)**2/eCt(i)**2
       chisqa= chisqa + c(i)
       h1a= h1a + deltah
    END DO

    ! Use initial h1 value of 0 to calculate start-point chisq.

    chisqb= 0
    DO i= 1,nd
       h1b= h1b + deltah
       Cs(i)= exp(-h1b*t(i))!*Ab !+ Bb*exp(-h2b*t(i))
       Ctdiff(i)= Ct(i) - Cs(i)
       c(i)= Ctdiff(i)**2/eCt(i)**2
       chisqb= chisqb + c(i)
    END DO

    ! First-step h1 used to find competing chisq for comparison.
    WRITE(6,*) 'Chi-squared a=', chisqa,'for Lamda1=',h1a
    WRITE(6,*) 'Chi-squared b=', chisqb,'for Lamda1=',h1b
    ! Prints the two calculated chisq values to screen.

    dchisq= chisqa - chisqb
    IF (dchisq.GT.0) THEN
       h1a= h1b
    ELSE IF (dchisq.LE.0) THEN
       deltah= (-deltah*2)/10
    END IF

    IF (chisqb.LE.6000) EXIT minloop

  END DO minloop

  WRITE(6,*) 'Chi-squared is', chisqb,'for Lamda1=',h1b

END PROGRAM assignment
于 2014-01-26T20:01:45.987 回答