(*          STEEPEST DESCENT ALGORITHM 10.3
*
*  To approximate a solution P to the minimization problem
*    G(P) = MIN( G(X) : X in R(n) )
*        given an initial approximation x
*
*   INPUT: Number n of variables; initial approximation x;
*	   tolerance TOL; Maximun number of iterations NN.
*
*  OUTPUT: Approximate solution X or a message of failure
*
*)
F[z1_,z2_,z3_,n_] := Module[{d,i,f},
   d = 0;
   For[i =1, i <= n, i++,
      d = d+(CF[i][z1,z2,z3])^2;
   ];
   f = d;
   f
];
OK = 0;
n = 3;
For[i = 1,i <= n,i++,
   TEMP[i] = Input["This is the Steepest Descent Method.\n
    If n is not 3 the program must be altered.\n
    Input the function CF("<>ToString[i]<>") in terms of x1,x2,x3\n"];
   CF[i][x1_,x2_,x3_] := Evaluate[TEMP[i]];
];
   PT = 0;
   DER = D[TEMP[1],x1];
   PT = PT+2*TEMP[1]*DER;
   DER = D[TEMP[2],x1];
   PT = PT+2*TEMP[2]*DER;
   DER = D[TEMP[3],x1];
   PT = PT+2*TEMP[3]*DER;
   P[1][x1_,x2_,x3_] := Evaluate[PT];
   PT = 0;
   DER = D[TEMP[1],x2];
   PT = PT+2*TEMP[1]*DER;
   DER = D[TEMP[2],x2];
   PT = PT+2*TEMP[2]*DER;
   DER = D[TEMP[3],x2];
   PT = PT+2*TEMP[3]*DER;
   P[2][x1_,x2_,x3_] := Evaluate[PT];
   PT = 0;
   DER = D[TEMP[1],x3];
   PT = PT+2*TEMP[1]*DER;
   DER = D[TEMP[2],x3];
   PT = PT+2*TEMP[2]*DER;
   DER = D[TEMP[3],x3];
   PT = PT+2*TEMP[3]*DER;
   P[3][x1_,x2_,x3_] := Evaluate[PT];
OK = 0;
While[OK == 0,  
   TOL = Input["Input the tolerance.\n"];
   If[TOL > 0,
      OK = 1,
      Input["Tolerance must be positive\n
      \n
      Press 1 [enter] to continue\n"];
   ];
];
OK = 0;
While[OK == 0,
   NN = Input["Input maximum number of iterations\n
     no decimal points\n"];
   If[NN <= 0,
     Input["Must be a positive integer.\n
     Enter 1 [enter] to continue\n"],
     OK = 1;
   ];
];
For[i = 1, i <= n, i++,
   X[i-1] = 
   Input["Input the initial approximation X("<>ToString[i]<>")\n"];
];
If[OK == 1,
   FLAG1 = Input["Select output destination\n
      1. Screen\n
      2. Text file\n
    Enter 1 or 2\n"];
   If[FLAG1 == 2,
      NAME = InputString["Input the file name\n
         For example:   output.dta\n"];
      OUP = OpenWrite[NAME,FormatType->OutputForm],
      OUP = "stdout";
   ];
   FLAG1 = Input["Select amount of output\n
	1. Answer only\n
	2. All intermediate approximations\n
        enter 1 or 2\n"];
   Write[OUP,"STEEPEST DESCENT METHOD FOR NONLINEAR SYSTEMS\n"];
   If[FLAG1 == 2,
      Write[OUP,"Iteration,  Approximation\n"];
   ];
(* Step 1  *)
   K = 1;
(* Step 2  *)
   While[OK == 1 && K <= NN,
(* Step 3  *)
      G[0] = N[F[X[0],X[1],X[2],n]];
      Z0 = 0;
      For[i = 1, i <= n, i++,
	 Z[i-1] = N[P[i][X[0],X[1],X[2]]];
   	 Z0 = Z0+(Z[i-1])^2;
      ];
      Z0 = Sqrt[Z0];
(* Step 4  *)
      If[Z0 <= 10^-20,
	 OK = 0;
	 Write[OUP,"0 gradient - may have a minimum\n"],
(* Step 5  *)
	 For[i = 1, i <= n, i++,
	    Z[i-1] = Z[i-1]/Z0;
	 ];
	 A[0] = 0;
	 X0 = 1;
	 For[i = 1, i <= n, i++,
	    c[i-1] = X[i-1]-X0*Z[i-1];
	 ];
	 G0 = N[Evaluate[F[c[0],c[1],c[2],n]]];
(* Step 6  *)
	 FLAG = 1;
	 If[G0 < G[0],
	    FLAG = 0;
	 ];
	 While[FLAG == 1 && OK == 1,
(* Steps 7 and 8  *)
	    X0 = 0.5*X0;
	    If[Abs[X0] <= 10^-20,
	       OK = 0;
	       Write[OUP,"No likely improvement - may have a minimum\n"],
	       For[i = 1, i <= n, i++,
		  c[i-1] = X[i-1]-X0*Z[i-1];
	       ];
	       G0 = N[Evaluate[F[c[0],c[1],c[2],n]]];
	    ];
	    If[G0 < G[0],
	       FLAG = 0;
	    ];
         ];
	 If[OK == 1,
	    A[2] = X0;
	    G[2] = G0;
(* Step 9  *)  
	    X0 = 0.5*X0;
	    For[i = 1, i <= n, i++,
	       c[i-1] = X[i-1]-X0*Z[i-1];
	    ];
	    A[1] = X0;
	    G[1] = N[Evaluate[F[c[0],c[1],c[2],n]]];
(* Step 10  *)
	    H1 = (G[1]-G[0])/(A[1]-A[0]);
	    H2 = (G[2]-G[1])/(A[2]-A[1]);
	    H3 = (H2-H1)/(A[2]-A[0]);
(* Step 11  *)
	    X0 = 0.5*(A[0]+A[1]-H1/H3);
	    For[i = 1, i <= n, i++,
	       c[i-1] = X[i-1]-X0*Z[i-1];
	    ];
	    G0 = N[Evaluate[F[c[0],c[1],c[2],n]]];
(* Step 12  *)
	    A0 = X0;
	    For[i = 1, i <= n, i++,
	       If[Abs[G[i-1]] < Abs[G0],
	 	  A0 = A[i-1];
		  G0 = G[i-1];
	       ];
    	    ];
	    If[Abs[A0] <= 10^-20,
	       OK = 0;
	       Write[OUP,"No change likely\n"];
	       Write[OUP,"- probably rounding error problems\n"],
(* Step 13  *)
	       For[i = 1, i <= n, i++,
		  X[i-1] = Evaluate[X[i-1]-A0*Z[i-1]];
	       ];
(* Step 14  *)
	       If[FLAG1 == 2,
		  Write[OUP,K];
	   	  For[i = 1, i <= n, i++,
		     Write[OUP,SetPrecision[X[i-1],10]];
	   	  ];
	  	  Write[OUP,"\n"];
	       ];
	       If[Abs[G0] < TOL || Abs[G0-G[0]] < TOL,
		  OK = 0;
		  Write[OUP,"Iteration number \n",K];
		  Write[OUP,"gives solution \n"];
		  Write[OUP,"\n"];
		  For[i = 1, i <= n, i++,
		     Write[OUP,SetPrecision[X[i-1],10]];
	 	  ];
	          Write[OUP,"\n"];
		  Write[OUP,"\n"];
		  Write[OUP,"to within \n",SetPrecision[TOL,9]];
		  Write[OUP,"\n"];
		  Write[OUP,"Process is complete\n"],
(* Step 15  *)
		  K = K+1;
               ];
            ];
         ];
      ];
   ];
(* Step 16  *)
   If[K > NN,
      Write[OUP,"Process does not converge in ",NN," iterations\n"];
   ];
   If[OUP == "OutputStream[",NAME," 3]",
      Print["Output file: ",NAME," created successfully\n"];
      Close[OUP];
   ];
];
Quit[];