|
(* 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[z1,z2,z3])^2;
];
f = d;
f
];
OK = 0;
n = 3;
For[i = 1,i <= n,i++,
TEMP = Input["This is the Steepest Descent Method.\n
If n is not 3 the program must be altered.\n
Input the function CF("<>ToString<>") in terms of x1,x2,x3\n"];
CF[x1_,x2_,x3_] := Evaluate[TEMP];
];
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<>")\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[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,N[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,N[X[i-1],10]];
];
Write[OUP,"\n"];
Write[OUP,"\n"];
Write[OUP,"to within \n",TOL];
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];
];
]; |
评分
-
1
查看全部评分
-
|