amoeba.txt

用VB实现的函数的极致和最优化 包括黄金分割搜索法、不用导数的布伦特法 用导数的布伦特法、多元函数的山单纯形法、多元函数的变尺度法等

Sub AMOEBA(P(), Y(), MP, NP, NDIM, FTOL, ITER)
    NMAX = 20
    ALPHA = 1#
    BETA = 0.5
    GAMMA = 2#
    ITMAX = 500
    Dim PR(20), PRR(20), PBAR(20)
    MPTS = NDIM + 1
    ITER = 0
    Do
      ILO = 1
      If Y(1) > Y(2) Then
          IHI = 1
          INHI = 2
      Else
          IHI = 2
          INHI = 1
      End If
      For I = 1 To MPTS
          If Y(I) < Y(ILO) Then ILO = I
          If Y(I) > Y(IHI) Then
              INHI = IHI
              IHI = I
          ElseIf Y(I) > Y(INHI) Then
              If I <> IHI Then INHI = I
          End If
      Next I
      RTOL = 2# * Abs(Y(IHI) - Y(ILO)) / (Abs(Y(IHI)) + Abs(Y(ILO)))
      If RTOL < FTOL Then Erase PBAR, PRR, PR: Exit Sub
      If ITER = ITMAX Then
          Print " Amoeba exceeding maximum iterations."
          Exit Sub
      End If
      ITER = ITER + 1
      For J = 1 To NDIM
          PBAR(J) = 0#
      Next J
      For I = 1 To MPTS
          If I <> IHI Then
              For J = 1 To NDIM
                  PBAR(J) = PBAR(J) + P(I, J)
              Next J
          End If
      Next I
      For J = 1 To NDIM
          PBAR(J) = PBAR(J) / NDIM
          PR(J) = (1# + ALPHA) * PBAR(J) - ALPHA * P(IHI, J)
      Next J
      YPR = FAMOEB(PR)
      If YPR <= Y(ILO) Then
          For J = 1 To NDIM
              PRR(J) = GAMMA * PR(J) + (1# - GAMMA) * PBAR(J)
          Next J
          YPRR = FAMOEB(PRR)
          If YPRR < Y(ILO) Then
              For J = 1 To NDIM
                  P(IHI, J) = PRR(J)
              Next J
              Y(IHI) = YPRR
          Else
              For J = 1 To NDIM
                  P(IHI, J) = PR(J)
              Next J
              Y(IHI) = YPR
          End If
      ElseIf YPR >= Y(INHI) Then
          If YPR < Y(IHI) Then
              For J = 1 To NDIM
                  P(IHI, J) = PR(J)
              Next J
              Y(IHI) = YPR
          End If
          For J = 1 To NDIM
              PRR(J) = BETA * P(IHI, J) + (1# - BETA) * PBAR(J)
          Next J
          YPRR = FAMOEB(PRR)
          If YPRR < Y(IHI) Then
              For J = 1 To NDIM
                  P(IHI, J) = PRR(J)
              Next J
              Y(IHI) = YPRR
          Else
              For I = 1 To MPTS
                  If I <> ILO Then
                      For J = 1 To NDIM
                          PR(J) = 0.5 * (P(I, J) + P(ILO, J))
                          P(I, J) = PR(J)
                      Next J
                      Y(I) = FAMOEB(PR)
                  End If
              Next I
          End If
      Else
          For J = 1 To NDIM
              P(IHI, J) = PR(J)
          Next J
          Y(IHI) = YPR
      End If
    Loop
End Sub