program ALG103; { STEEPEST DESCENT METHOD To approximate a solution P to the minimization problem G(P) = MIN( G(X) : X IN n-DIM ) given an initial approximation X: INPUT: Number n of variables; initial approximation X; tolerance TOL; maximum number of iterations N. OUTPUT: Approximate solution X or a message of failure. } const ZERO = 1.0E-20; type VEC = array [ 1 .. 10] of real; var X,Z,C : VEC; G,A : array [1 .. 4] of real; Z0,X0,G0,H1,H2,H3,A0,TOL : real; N,NN,K,I,FLAG1 : integer; FLAG,OK : boolean; AA : char; OUP : text; NAME : string [ 14 ]; { Change procedures CF and P for a new problem } function CF ( I : integer; var X : VEC ) : real; begin case I of 1 : CF := 3*x[1] - cos(x[2]*x[3]) - 0.5; 2 : CF := x[1]*x[1] - 81*sqr(x[2]+0.1) + sin(x[3]) + 1.06; 3 : CF := exp(-x[1]*x[2]) + 20*x[3] + (10*pi-3)/3 end end; function P ( I : integer; var X : VEC ) : real; begin case I of 1 : P := 2 * ( 3 ) * CF( 1, X ) + 2 * ( 2*x[1] ) * CF( 2, X ) + 2 * ( -x[2]*exp(-x[1]*x[2]) ) * CF( 3, X ); 2 : P := 2 * ( X[3]*sin(X[2]*X[3]) ) * CF( 1, X ) + 2 * ( -162*(X[2]+0.1) ) * CF( 2, X ) + 2 * ( -x[1]*exp(-x[1]*x[2]) ) * CF( 3, X ); 3 : P := 2 * ( X[2]*sin(X[2]*X[3]) ) * CF( 1, X ) + 2 * ( cos(X[3]) ) * CF( 2, X ) + 2 * ( 20 ) * CF( 3, X ) end end; function F ( var X : VEC ) : real; var D : real; I : integer; begin D := 0.0; for I := 1 to N do D := D + sqr( CF( I, X ) ); F := D end; procedure INPUT; begin writeln('This is the Steepest Descent Method.'); OK := false; writeln; writeln('NOTE THAT THE FUNCTION DEFINITIONS ARE VERY COMPLICATED'); writeln ('Have the functions been set up as follows: '); writeln; writeln (' 1. F( var X : VEC ) : real '); writeln (' which is the function to be minimized ( function G ) '); writeln (' or the sum of the squares of the component'); writeln (' functions in a nonlinear system; '); writeln; writeln (' 2. CF( I : integer; var X : VEC ) : real '); writeln (' which is the Ith component function of a nonlinear '); writeln (' system - not used for simple minimization problems '); writeln; writeln (' 3. P( I : integer; var X : VEC ) : real '); writeln (' which is the partial derivative of F with respect '); writeln (' to the Ith variable '); writeln; writeln ('Answer Y or N. '); writeln; readln ( AA ); if ( AA = 'Y' ) or ( AA = 'y' ) then begin while ( not OK ) do begin writeln ('Input the number n of equations. '); readln ( N ); if ( N >= 2 ) then OK := true else writeln ('N must be an integer greater than 1. ') end; OK := false; while ( not OK ) do begin writeln ('Input the Tolerance. '); readln ( TOL ); if ( TOL > 0.0 ) then OK := true else writeln ('Tolerance must be positive. ') end; OK := false; while ( not OK ) do begin writeln ('Input the maximum number of iterations. '); readln ( NN ); if ( NN > 0 ) then OK := true else writeln ('Must be a positive integer. ') end; for I := 1 to N do begin writeln ('Input initial approximation X(', I, ').' ); readln ( X[I] ) end end else writeln ('The program will end so that the functions can be defined.') end; procedure OUTPUT; begin writeln ('Select output destination '); writeln ('1. Screen '); writeln ('2. Text file '); writeln ('Enter 1 or 2 '); readln ( FLAG1 ); if ( FLAG1 = 2 ) then begin write ('Input the file name in the form - '); writeln ('drive:name.ext '); writeln ('for example: A:OUTPUT.DTA '); readln ( NAME ); assign ( OUP, NAME ) end else assign ( OUP, 'CON'); rewrite ( OUP ); writeln ('Select amount of output '); writeln ('1. Answer only '); writeln ('2. All intermediate approximations '); writeln ('Enter 1 or 2 '); readln (FLAG1); writeln(OUP,'STEEPEST DESCENT METHOD FOR NONLINEAR SYSTEMS'); writeln(OUP); if FLAG1 = 2 then begin writeln(OUP,'Iteration, Approximation') end end; begin INPUT; if ( OK ) then begin OUTPUT; K := 1; { STEP 2 } while ( ( OK ) and ( K <= NN ) ) do begin { STEP 3 } G[1] := F( X ); Z0 := 0.0; for I := 1 to N do begin Z[I] := P( I, X ); Z0 := Z0 + sqr( Z[I] ) end; Z0 := sqrt( Z0 ); { STEP 4 } if ( Z0 <= ZERO ) then begin OK := false; writeln (OUP,'Zero qradient - may have a minimum '); end else begin { STEP 5 } for I := 1 to N do Z[I] := Z[I] / Z0; A[1] := 0.0; X0 := 1.0; for I := 1 to N do C[I] := X[I] - X0 * Z[I]; G0 := F( C ); { STEP 6 } FLAG := true; while ( FLAG and OK ) do begin if G0 < G[1] then FLAG := false; { STEPS 7 AND 8 } X0 := 0.5 * X0; if ( X0 <= ZERO ) then begin OK := false; writeln (OUP,'No likely improvement - may '); writeln (OUP,'have a minimum '); end else begin for I := 1 to N do C[I] := X[I]-X0*Z[I]; G0 := F( C ); end end; if ( OK ) then begin A[3] := X0; G[3] := G0; { STEP 9 } X0 := 0.5 * X0; for I := 1 to N do C[I] := X[I]-X0*Z[I]; A[2] := X0; G[2] := F( 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]); { STEP 11 } X0 := 0.5*(A[1]+A[2]-H1/H3); for I := 1 to N do C[I] := X[I]-X0*Z[I]; G0 := F( C ); { STEP 12 } A0 := X0; for I := 1 to N do if ( abs( G[I] ) < abs( G0 ) ) then begin A0 := A[I]; G0 := G[I] end; if ( abs( A0 ) <= ZERO ) then begin OK := false; writeln (OUP,'No change likely '); write (OUP,'- probably rounding error '); writeln (OUP,'problems '); end else begin { STEP 13 } for I := 1 to N do X[I] := X[I]-A0*Z[I]; { STEP 14 } if (FLAG1 = 2) then begin write(OUP,K:3); for I := 1 to N do write(OUP,X[I]:12:8); writeln(OUP); end; if ((abs(G0) < TOL) or (abs(G0-G[1]) < TOL)) then begin OK := false; writeln(OUP,'Iteration number ',K); writeln(OUP,'gives solution'); writeln(OUP); for I := 1 to N do write(OUP,X[I]:12:8); writeln(OUP); writeln(OUP); writeln(OUP,'to within ',TOL); writeln(OUP); writeln (OUP,'Process is complete '); end else { STEP 15 } K := K + 1 end end end end; if ( K > NN ) then begin { STEP 16 } writeln (OUP,'Process does not converge in ',NN ); write (OUP,' iterations ') end; close(OUP) end end.