Attribute VB_Name = "ModSteepestAlgorithm"
 '????????????????????????????????
Option Explicit
 Public Sub SteepEstMethod(canvas As PictureBox, AAA As String, nEq As Integer, stEqs() As String, VvarsSt() As String, VvarsI() As Double, NTry As Integer, Tol As Double, XAns() As Double)


'C***********************************************************************
'C                                                                      *
'C                  STEEPEST DESCENT ALGORITHM 10.3                     *
'C                                                                      *
'C***********************************************************************
'C
'C
'C
'C     TO APPROXIMATE A SOLUTION P TO THE MINIMIZATION PROBLEM
'C                 G(P) = MIN ( G(X) : X IN R(n) )
'C     GIVEN AN INITIAL APPROXIMATION X:
'C
'C     INPUT NUMBER N OF VARIABLES; INITIAL APPROXIMATION X;
'C           TOLERANCE TOL; MAXIMUM NUMBER OF ITERATIONS N.
'C
'C     OUTPUT APPROXIMATE SOLUTION X OR A MESSAGE OF FAILURE.

'C
      Dim FLAG As Boolean
      Dim Xx() As Double, Z() As Double, C() As Double, G() As Double, A() As Double
      ReDim Xx(1 To nEq), Z(1 To nEq), C(1 To nEq), G(1 To nEq), A(1 To nEq)
      'Dim FLAG1 As Boolean
      Dim OK As Boolean, K As Integer, i As Integer, j As Integer
      Dim Z0 As Double, X0 As Double, G0 As Double, G1 As Double, A0 As Double
      Dim Zero As Double
      Dim H1 As Double, H2 As Double, H3 As Double
      
'C     THE FUNCTION F USED HERE CORRESPONDS TO G USED IN THE ALGORITHM
      'pi = 4 * Atan(1#)
      Zero = 1E-20
      'OPEN(UNIT=5,FILE='CON',ACCESS='SEQUENTIAL')
      'OPEN(UNIT=6,FILE='CON',ACCESS='SEQUENTIAL')
      'WRITE(6,*) 'This is the Steepest Descent Method'
      'WRITE(6,*) 'NOTE THAT THE FUNCTIONS ARE VERY COMPLICATED.
      'WRITE(6,*) 'Have the functions F1,...,Fn and F (G) been'
      'WRITE(6,*) 'defined at the beginning of the program and.'
      'WRITE(6,*) 'has the gradient been implemented at line 118.'
      'WRITE(6,*) 'Enter Y or N '
      'WRITE(6,*) ' '
      'READ(5,*)  AAA
      If LCase(AAA) = "y" Then
         OK = False
101      If OK Then GoTo 11
            'WRITE(6,*) 'Input the number n of equations.'
            'WRITE(6,*) ' '
            'READ(5,*) N
            If nEq >= 2 Then
               OK = True
            Else
               MsgBox ("Neq must be greater than 1.")
            End If
         GoTo 101
11       OK = False
14       If (OK) Then GoTo 15
            'WRITE(6,*) 'Input tolerance '
            'WRITE(6,*) ' '
            'READ(5,*) TOL
            If Tol <= 0# Then
               MsgBox ("Tolerance must be positive ")
            Else
               OK = True
            End If
         GoTo 14
15       OK = False
12       If OK Then GoTo 13
            'WRITE(6,*) 'Input the maximum number of iterations.'
            'WRITE(6,*) ' '
            'READ(5,*) M
            If (NTry <= 0) Then
              MsgBox ("Niteration Must be positive integer ")
            Else
              OK = True
            End If
         GoTo 12
13
'13       For I = 1 To NTry
            'WRITE(6,*) 'Input initial approximation X(',I,')'
            'WRITE(6,*) ' '
            'READ(5,*) X(I)
'16       CONTINUE
      Else
         MsgBox ("The program will end so that the functions can be created ")
         OK = False
      End If
      If Not (OK) Then GoTo 400
      'WRITE(6,*) 'Select output destination: '
      'WRITE(6,*) '1. Screen '
      'WRITE(6,*) '2. Text file '
     ' WRITE(6,*) 'Enter 1 or 2 '
      'WRITE(6,*) ' '
      'READ(5,*) FLAG1
      
      canvas.Print "STEEPEST DESCENT METHOD FOR NONLINEAR SYSTEMS"
      canvas.Print "Iter. No.   Vector    G(vector)"
'C     STEP 1
      K = 0
      For i = 1 To nEq 'DO 16 I=1,NEq
            'WRITE(6,*) 'Input initial approximation X(',I,')'
            'WRITE(6,*) ' '
            'READ(5,*) X(I)
            Xx(i) = VvarsI(i)
           ' canvas.Print "i=" & i & "x(" & i & ")=" & Xx(i)
16       Next i ' CONTINUE
     'canvas.Print (" K= " & k & ";X(1)= " & X(1) & ";X(2)= " & X(2) & ";X(3)= " & X(3))
      
      K = 1
'C     STEP 2
100   If (K > NTry) Then GoTo 200
       'MsgBox ("kstep = " & K)
'C        STEP 3
      '////////////////////////////////
      'Dim Fxyz1 As Double, Fxyz2 As Double, Fxyz3 As Double, fxyz As Double
      'Fxyz1 = 3 * Xx(1) - Cos(Xx(2) * Xx(3)) - 0.5 '在點(x1,x2,x3)處方程式1值
      'Fxyz2 = Xx(1) * Xx(1) - 81 * (Xx(2) + 0.1) ^ 2 + Sin(Xx(3)) + 1.06
      'Fxyz3 = Exp(-Xx(1) * Xx(2)) + 20 * Xx(3) + (10 * PI - 3) / 3
      'F(X1,X2,X3)=F1(X1,X2,X3)**2+F2(X1,X2,X3)**2+F3(X1,X2,X3)**2 '在點(x1,x2,x3)處方程式1,方程式2,方程式3平方和
      'fxyz = Fxyz1 ^ 2 + Fxyz2 ^ 2 + Fxyz3 ^ 2
      'canvas.Print "fxyz1~3  A=" & Fxyz1 & ";" & Fxyz2 & ";" & Fxyz3 & ";f^2 A=" & fxyz
      '//////////////////////////////


         'C(1) = F1(X(1), X(2), X(3))    '在點(x1,x2,x3)處方程式1值
         'C(2) = F2(X(1), X(2), X(3))    '在點(x1,x2,x3)處方程式2值
         'C(3) = F3(X(1), X(2), X(3))    '在點(x1,x2,x3)處方程式3值
          'F(X1,X2,X3)=F1(X1,X2,X3)**2+F2(X1,X2,X3)**2+F3(X1,X2,X3)**2 '在點(x1,x2,x3)處方程式1,方程式2,方程式3平方和
         Dim sumF As Double
         sumF = 0
         For i = 1 To nEq
         C(i) = FunctionValues(glVbScriptA(i - 1), stEqs(i), VvarsSt, Xx)
         sumF = sumF + C(i) * C(i)
         Next i
         'canvas.Print "fxyz1~3 B =" & C(1) & ";" & C(2) & ";" & C(3) & ";f^2 B=" & sumF
'C        THE VECTOR C HOLDS F(X) FROM THE EXAMPLE AND G1 HOLDS G(X)
        ' G1 = F(X(1), X(2), X(3))
        
         
         G(1) = sumF 'G(1) = G1
         'Dim Q As Double, P As Double, Z1 As Double, Z2 As Double, Z3 As Double
         'Q = Exp(-Xx(1) * Xx(2))
         'P = Sin(Xx(2) * Xx(3))
         'Z1 = 6 * c(1) + 4 * Xx(1) * c(2) - 2 * Xx(2) * Q * c(3)
         'Z2 = 2 * Xx(3) * P * c(1) - 324 * (Xx(2) + 0.1) * c(2) - 2 * Xx(1) * Q * c(3)
         'Z3 = 2 * Xx(2) * P * c(1) + 2 * Cos(Xx(3)) * c(2) + 40 * c(3)
         ' canvas.Print "Z1,Z2,Z3 = " & Z1 & ";" & Z2 & ";" & Z3
         
         
         
         Dim vTpt As Double, zoSum2 As Double
         zoSum2 = 0
         For i = 1 To nEq
            vTpt = 0
            For j = 1 To nEq
             vTpt = vTpt + 2 * C(j) * PartialDif(glVbScriptA(j - 1), stEqs(j), VvarsSt, Xx, i)
            Next j
             Z(i) = vTpt
             zoSum2 = zoSum2 + Z(i) * Z(i)
         Next i
         
        ' canvas.Print "Z(1),Z(2),Z(3)  =" & Z(1) & ";" & Z(2) & ";" & Z(3)
         vTpt = 0
         'Dim Z0A As Double
         'Z0A = Sqr(Z1 ^ 2 + Z2 ^ 2 + Z3 ^ 2)
         Z0 = Sqr(zoSum2)
         '?canvas.Print "Z0A,Z0=" & Z0A & ";" & Z0
'C        STEP 4
         If (Abs(Z0) < 1E-20) Then
          MsgBox ("ZERO GRADIENT: MAY HAVE MIN")
            GoTo 400
         End If
'C        STEP 5
         For j = 1 To nEq 'DO 10 J=1,N
10       Z(j) = Z(j) / Z0
         Next j
         
'C        A IS USED FOR ALPHA AND G0 FOR G(3)
         A(1) = 0
         X0 = 1
         '????????????????????????????????????////
         'Dim G0a As Double
         'G0a = Fxyz_T(Xx(1) - X0 * Z(1), Xx(2) - X0 * Z(2), Xx(3) - X0 * Z(3))   '在點 x(1) - X0 * Z(1), x(2) - X0 * Z(2), x(3) - X0 * Z(3)處方程式1,方程式2,方程式3平方和
         Dim xXtpt() As Double
         ReDim xXtpt(1 To nEq)
         For i = 1 To nEq
            xXtpt(i) = Xx(i) - X0 * Z(i)
         Next i
         
          sumF = 0
          For i = 1 To nEq
          sumF = sumF + (FunctionValues(glVbScriptA(i - 1), stEqs(i), VvarsSt, xXtpt)) ^ 2
          Next i
          G0 = sumF
          '?canvas.Print "g0a,g0=" & G0a & ";" & G0
          '????????????????????????????????????????
'C        STEP 6
         FLAG = True
         If (Abs(G0) < Abs(G1)) Then FLAG = False
500      If Not (FLAG) Then GoTo 600
'C           STEPS 7 AND 8
            X0 = X0 / 2
           ' MsgBox ("x0=" & X0 & "zero= " & Zero)
            If (X0 < Zero) Then MsgBox ("NO LIKELY IMPROVEMENT: MAY HAVE MIN")
               
             If (X0 < Zero) Then GoTo 400
            
            'G0 = F(x(1) - X0 * Z(1), x(2) - X0 * Z(2), x(3) - X0 * Z(3))
            '????????????????????????????????????////
            'Dim G0b As Double
            'G0b = Fxyz_T(Xx(1) - X0 * Z(1), Xx(2) - X0 * Z(2), Xx(3) - X0 * Z(3))
            
         For i = 1 To nEq
            xXtpt(i) = Xx(i) - X0 * Z(i)
         Next i
         'G0 = F(x(1) - X0 * Z(1), x(2) - X0 * Z(2), x(3) - X0 * Z(3))   '在點 x(1) - X0 * Z(1), x(2) - X0 * Z(2), x(3) - X0 * Z(3)處方程式1,方程式2,方程式3平方和
          sumF = 0
          For i = 1 To nEq
          sumF = sumF + (FunctionValues(glVbScriptA(i - 1), stEqs(i), VvarsSt, xXtpt)) ^ 2
          Next i
          G0 = sumF
         ' canvas.Print "g0b,g0=" & G0b & ";" & G0
          '????????????????????????????????????????
            
            If (Abs(G0) < Abs(G1)) Then FLAG = False
            If (Abs(G0) < Abs(G1)) Then GoTo 500
        
600      A(3) = X0
         G(3) = G0
'C        STEP 9
         X0 = X0 / 2
         'G(2) = F(x(1) - X0 * Z(1), x(2) - X0 * Z(2), x(3) - X0 * Z(3))
         For i = 1 To nEq
            xXtpt(i) = Xx(i) - X0 * Z(i)
         Next i
         
          sumF = 0
          For i = 1 To nEq
          sumF = sumF + (FunctionValues(glVbScriptA(i - 1), stEqs(i), VvarsSt, xXtpt)) ^ 2
          Next i
          G(2) = sumF
          A(2) = X0
'C        STEP 10

         H1 = (G(2) - G(1)) / (A(2) - A(1))
         H2 = (G(3) - G(2)) / (A(3) - A(2))
         H3 = (H2 - H1) / (A(3) - A(1))
'C        STEP 11
         X0 = 0.5 * (A(1) + A(2) - H1 / H3)
         'G0 = F(x(1) - X0 * Z(1), x(2) - X0 * Z(2), x(3) - X0 * Z(3))
         
          '????????????????????????????????????////
      
         For i = 1 To nEq
            xXtpt(i) = Xx(i) - X0 * Z(i)
         Next i
         
          sumF = 0
          For i = 1 To nEq
          sumF = sumF + (FunctionValues(glVbScriptA(i - 1), stEqs(i), VvarsSt, xXtpt)) ^ 2
          Next i
          G0 = sumF
          '????????????????????????????????????????
'C        STEP 12
         A0 = X0
         For j = 1 To nEq 'DO 20 J=1,N
           If (Abs(G(j)) < Abs(G0)) Then
               A0 = A(j)
               G0 = G(j)
           End If
20       'CONTINUE
         Next j
         'MsgBox ("pass 20;k= " & K)
         If Abs(A0) < Zero Then
            MsgBox ("NO CHANGE LIKELY: PROBABLY ROUNDING DIFFICULTIES")
            
            GoTo 400
         End If
'C        STEP 13
         For j = 1 To nEq 'DO 30 J=1,N
30       Xx(j) = Xx(j) - A0 * Z(j) '
         ' canvas.Print " K= " & K & ";X(j)= " & Xx(j)
         Next j
         
         'canvas.Print " K= " & K & ";X(1)= " & Xx(1) & ";X(2)= " & Xx(2) & ";X(3)= " & Xx(3)
         
'C        STEP 14
         If (Abs(G0) < Tol) Or (Abs(G0 - G1) < Tol) Then
            MsgBox ("PROCEDURE COMPLETED SUCCESSFULLY")
            
            GoTo 400
         End If
'C        STEP 15
         K = K + 1
         'MsgBox ("k=" & K)
      GoTo 100
'C     STEP 16
Exit Sub
200   MsgBox ("MAXIMUM ITERATIONS EXCEEDED")
      For i = 1 To nEq
       ReDim Preserve XAns(1 To i)
       XAns(i) = 0 'Xx(i)
       Next i
400
       canvas.Print "sta 400 K= " & K & ";X(1)= " & Xx(1) & ";X(2)= " & Xx(2) & ";X(3)= " & Xx(3)
       For i = 1 To nEq
       ReDim Preserve XAns(1 To i)
       XAns(i) = Xx(i)
       Next i
      End Sub
 
       Private Function Fxyz_1(X1, X2, X3) As Double
        Fxyz_1 = 3 * X1 - Cos(X2 * X3) - 0.5
       End Function
       Private Function Fxyz_2(X1, X2, X3) As Double
       Fxyz_2 = X1 * X1 - 81 * (X2 + 0.1) ^ 2 + Sin(X3) + 1.06
       End Function
       
      Private Function Fxyz_3(X1, X2, X3) As Double
         Fxyz_3 = Exp(-X1 * X2) + 20 * X3 + (10 * PI - 3) / 3
      End Function
     Private Function Fxyz_T(X1, X2, X3)
       Fxyz_T = Fxyz_1(X1, X2, X3) ^ 2 + Fxyz_2(X1, X2, X3) ^ 2 + Fxyz_3(X1, X2, X3) ^ 2
     End Function
      
