> restart:
> # 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: maximum number of iterations N.
> #
> # OUTPUT:  Approximate solution X or a message of failure.
> F := proc(X,N) local D, I2, f:
> D := 0:
> for I2 from 1 to N do
> D := D+(CF[I2](seq(X[i],i=0..N-1)))*(CF[I2](seq(X[i],i=0..N-1))):
> od:
> f := D:
> RETURN(evalf(f)):
> end:
> print(`This is the Steepest Descent Method.`):
> OK := FALSE:
> while OK = FALSE do
> print(`Input the number n of equations.`):
> N := scanf(`%d`)[1]: print(`n = `):print(N):
> if N >= 2 then
> OK := TRUE:
> else
> print(`N must be an integer greater than 1.`):
> fi:
> od:
> print(`The functions CF represent the components of F(X).`):
> for I2 from 1 to N do
> print(`Input the function CF[I](x1..xn) for I = `,I2):
> CF[I2] := scanf(`%a`)[1]: print(`Function is `):print(CF[I2]):
> od:
> for I2 from 1 to N do
> P[I2] := 0:
> for J from 1 to N do
> P[I2] := P[I2]+2*CF[J]*diff(CF[J],evaln(x || I2)):
> od:
> P[I2] := unapply(P[I2],evaln(x || (1..N))):
> od:
> for I2 from 1 to N do
> CF[I2] := unapply(CF[I2],evaln(x || (1..N))):
> od:
> OK := FALSE:
> while OK = FALSE do
> print(`Input tolerance\n`):
> TOL := scanf(`%f`)[1]:print(`Tolerance = `):print(TOL):
> if TOL > 0 then
> OK := TRUE:
> else
> print(`Tolerance must be positive.`):
> fi:
> od:
> OK := FALSE:
> while OK = FALSE do
> print(`Input the maximum number of iterations.`):
> NN := scanf(`%d`)[1]:print(`Maximum number of iterations = `):print(NN):
> if NN > 0 then
> OK := TRUE:
> else
> print(`Must be a positive integer.`):
> fi:
> od:
> for I2 from 1 to N do
> print(`Input initial approximation X(I) for I = `, I2):
> X[I2-1] := scanf(`%f`)[1]:print(`Input = `):print(X[I2-1]):
> od:
> if OK = TRUE then
> print(`Select output destination`):
> print(`1. Screen`):
> print(`2. Text file`):
> print(`Enter 1 or 2`):
> FLAG1 := scanf(`%d`)[1]:print(`Your response is `):print(FLAG1):
> if FLAG1 = 2 then
> print(`Input the file name in the form - drive:\\name.ext`):
> print(`for example  A:\\OUTPUT.DTA`):
> NAME := scanf(`%s`)[1]:print(`Output file is `):print(NAME):
> OUP := fopen(NAME,WRITE,TEXT):
> else
> OUP := default:
> fi:
> print(`Select amount of output`):
> print(`1. Answer only`):
> print(`2. All intermediate approximations`):
> print(`Enter 1 or 2`):
> FLAG1 := scanf(`%d`)[1]:print(`Your response is `):print(FLAG1):
> fprintf(OUP, `STEEPEST DESCENT METHOD FOR NONLINEAR SYSTEMS\n\n`):
> if FLAG1 = 2 then
> fprintf(OUP, `Iteration, Approximation\n`):
> fi:
> # Step 1
> K := 1:
> # Step 2
> while OK = TRUE and K <= NN do
> # Step 3
> G[0] := evalf(F(X, N)):
> Z0 := 0:
> for I2 from 1 to N do
> Z[I2-1] := evalf(P[I2](seq(X[i],i=0..N-1))):
> Z0 := Z0+(Z[I2-1])*(Z[I2-1]):
> od:
> Z0 := sqrt(Z0):
> # Step 4
> if Z0 <= 1.0e-20 then
> OK := FALSE:
> fprintf(OUP, `0 qradient - may have a minimum\n`):
> else
> # Step 5
> for I2 from 1 to N do
> Z[I2-1] := Z[I2-1] / Z0:
> od:
> A[0] := 0:
> X0 := 1:
> for I2 from 1 to N do
> C[I2-1] := X[I2-1]-X0*Z[I2-1]:
> od:
> G0 := evalf(F(C, N)):
> # Step 6
> FLAG := TRUE:
> if G0 < G[0] then
> FLAG := FALSE:
> fi:
> while FLAG = TRUE and OK = TRUE do
> # Steps 7 and 8
> X0 := 0.5*X0:
> if X0 <= 1.0e-20 then
> OK := FALSE:
> fprintf(OUP, `No likely improvement - may\n`):
> fprintf(OUP, `have a minimum\n`):
> else
> for I2 from 1 to N do
> C[I2-1] := X[I2-1]-X0*Z[I2-1]:
> od:
> G0 := evalf(F(C, N)):
> fi:
> if G0 < G[0] then
> FLAG := FALSE:
> fi:
> od:
> if OK = TRUE then
> A[2] := X0:
> G[2] := G0:
> # Step 9
> X0 := 0.5*X0:
> for I2 from 1 to N do
> C[I2-1] := X[I2-1]-X0*Z[I2-1]:
> od:
> A[1] := X0:
> G[1] := evalf(F(C, 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 I2 from 1 to N do
> C[I2-1] := X[I2-1]-X0*Z[I2-1]:
> od:
> G0 := evalf(F(C, N)):
> # Step 12
> A0 := X0:
> for I2 from 1 to N do
> if abs(G[I2-1]) < abs(G0) then
> A0 := A[I2-1]:
> G0 := G[I2-1]:
> fi:
> od:
> if abs(A0) <= 1.0e-20 then
> OK := FALSE:
> fprintf(OUP, `No change likely\n`):
> fprintf(OUP, `- probably rounding error problems\n`):
> else
> # Step 13
> for I2 from 1 to N do
> X[I2-1] := evalf(X[I2-1]-A0*Z[I2-1]):
> od:
> # Step 14
> if FLAG1 = 2 then
> fprintf(OUP, ` %2d`, K):
> for I2 from 1 to N do
> fprintf(OUP, ` %11.8f`, X[I2-1]):
> od:
> fprintf(OUP, `\n`):
> fi:
> if abs(G0) < TOL or abs(G0-G[0]) < TOL then
> OK := FALSE:
> fprintf(OUP, `Iteration number %d\n`, K):
> fprintf(OUP, `gives solution\n\n`):
> for I2 from 1 to N do
> fprintf(OUP, ` %11.8f`, X[I2-1]):
> od:
> fprintf(OUP, `\n\nto within %.10e\n\n`, TOL):
> fprintf(OUP, `Process is complete\n`):
> else
> # Step 15
> K := K+1:
> fi:
> fi:
> fi:
> fi:
> od:
> if K > NN then
> # Step 16
> fprintf(OUP, `Process does not converge in %d\n`, NN):
> fprintf(OUP, ` iterations\n`):
> fi:
> if OUP <> default then
> fclose(OUP):
> print(`Output file `,NAME,` created successfully`):
> fi:
> fi: