alg103.txt

Numerical Anaysis 8th Edition Burden and Faires (Maple Source)

> 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.
> alg103 := proc() global F,CF; local OK, N, I, P, J, TOL, NN, X, FLAG1, NAME, OUP, K, G, Z0, Z, A, X0, C, G0, FLAG, H1, H2, H3, A0;
> F := proc(X,N) local D, I, f;
> D := 0;
> for I from 1 to N do
> D := D+(CF[I](seq(X[i],i=0..N-1)))*(CF[I](seq(X[i],i=0..N-1)));
> od;
> f := D;
> RETURN(evalf(f));
> end;
> printf(`This is the Steepest Descent Method.\n`);
> OK := FALSE;
> while OK = FALSE do
> printf(`Input the number n of equations.\n`);
> N := scanf(`%d`)[1];
> if N >= 2 then
> OK := TRUE;
> else
> printf(`N must be an integer greater than 1.\n`);
> fi;
> od;
> printf(`The functions CF represent the components of F(X).\n`);
> for I from 1 to N do
> printf(`Input the function CF[%d](x1..x%d).\n`,I, N);
> CF[I] := scanf(`%a`)[1];
> od;
> for I from 1 to N do
> P[I] := 0;
> for J from 1 to N do
> P[I] := P[I]+2*CF[J]*diff(CF[J],evaln(x || I));
> od;
> P[I] := unapply(P[I],evaln(x || (1..N)));
> od;
> for I from 1 to N do
> CF[I] := unapply(CF[I],evaln(x || (1..N)));
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input tolerance\n`);
> TOL := scanf(`%f`)[1];
> if TOL > 0 then
> OK := TRUE;
> else
> printf(`Tolerance must be positive.\n`);
> fi;
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input the maximum number of iterations.\n`);
> NN := scanf(`%d`)[1];
> if NN > 0 then
> OK := TRUE;
> else
> printf(`Must be a positive integer.\n`);
> fi;
> od;
> for I from 1 to N do
> printf(`Input initial approximation X(%d).\n`, I);
> X[I-1] := scanf(`%f`)[1];
> od;
> if OK = TRUE then
> printf(`Select output destination\n`);
> printf(`1. Screen\n`);
> printf(`2. Text file\n`);
> printf(`Enter 1 or 2\n`);
> FLAG1 := scanf(`%d`)[1];
> if FLAG1 = 2 then
> printf(`Input the file name in the form - drive:\\name.ext\n`);
> printf(`for example  A:\\OUTPUT.DTA\n`);
> NAME := scanf(`%s`)[1];
> OUP := fopen(NAME,WRITE,TEXT);
> else
> OUP := default;
> fi;
> printf(`Select amount of output\n`);
> printf(`1. Answer only\n`);
> printf(`2. All intermediate approximations\n`);
> printf(`Enter 1 or 2\n`);
> FLAG1 := scanf(`%d`)[1];
> 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 I from 1 to N do
> Z[I-1] := evalf(P[I](seq(X[i],i=0..N-1)));
> Z0 := Z0+(Z[I-1])*(Z[I-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 I from 1 to N do
> Z[I-1] := Z[I-1] / Z0;
> od;
> A[0] := 0;
> X0 := 1;
> for I from 1 to N do
> C[I-1] := X[I-1]-X0*Z[I-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 I from 1 to N do
> C[I-1] := X[I-1]-X0*Z[I-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 I from 1 to N do
> C[I-1] := X[I-1]-X0*Z[I-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 I from 1 to N do
> C[I-1] := X[I-1]-X0*Z[I-1];
> od;
> G0 := evalf(F(C, N));
> # Step 12
> A0 := X0;
> for I from 1 to N do
> if abs(G[I-1]) < abs(G0) then
> A0 := A[I-1];
> G0 := G[I-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 I from 1 to N do
> X[I-1] := evalf(X[I-1]-A0*Z[I-1]);
> od;
> # Step 14
> if FLAG1 = 2 then
> fprintf(OUP, ` %2d`, K);
> for I from 1 to N do
> fprintf(OUP, ` %11.8f`, X[I-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 I from 1 to N do
> fprintf(OUP, ` %11.8f`, X[I-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):
> printf(`Output file %s created successfully`,NAME);
> fi;
> fi;
> RETURN(0);
> end;
> alg103();