home *** CD-ROM | disk | FTP | other *** search
/ Collection of Education / collectionofeducationcarat1997.iso / SCIENCE / PRAXIS.ZIP / PRAXIS.I < prev    next >
Text File  |  1987-07-15  |  23KB  |  548 lines
  1. (****************************************************************************)
  2. (*              p r o c e d u r e      p r a x i s                          *)
  3. (*                                                                          *)
  4. (*   p r a x i s  is for the minimization of a function in several          *)
  5. (*   variables. the algorithm used is a modification of a conjugate         *)
  6. (*   gradient method developed by powell. changes are due to brent,         *)
  7. (*   who gives an algol w program.                                          *)
  8. (*   users who are interested in more of the details should read            *)
  9. (*          - powell, m. j. d., 1962. an efficient method for finding       *)
  10. (*            the minimum of a function of several variables without        *)
  11. (*            calculating derivatives, computer journal 7, 155-162.         *)
  12. (*          - brent, r. p., 1973. algorithms for minimization without       *)
  13. (*            derivatives. prentice hall, englewood cliffs.                 *)
  14. (*   if you have any further comments, questions, etc. please feel free     *)
  15. (*   and contact me.                                                        *)
  16. (*                                                                          *)
  17. (*                           karl gegenfurtner     02/17/86                 *)
  18. (*                                                 turbo - p version 1.0    *)
  19. (****************************************************************************)
  20. (*  the use of praxis is fairly simple. there are only two parameters:      *)
  21. (*      - x is of type vector and contains on input an initial guess        *)
  22. (*        of the solution. on output x holds the final solution of the      *)
  23. (*        system of equations.                                              *)
  24. (*      - n is of type integer and gives the number of unknown parameters   *)
  25. (*  the result of praxis is the least calculated value of the function f    *)
  26. (*  f is one of the global parameters used by praxis:                       *)
  27. (*    - f(x,n) is declared forward and is the function to be minimized      *)
  28. (*  all other globals are optional, i.e. you can change them or else        *)
  29. (*  praxis will use some default values, which are adequate for most        *)
  30. (*  problems under consideration.                                           *)
  31. (*    - prin controls the printout from praxis                              *)
  32. (*           0:  no printout at all                                         *)
  33. (*           1:  only initial and final values                              *)
  34. (*           2:  detailed map of the minimization process                   *)
  35. (*           3:  also prints eigenvalues and eigenvectors of the            *)
  36. (*               direction matices in use (for insiders only).              *)
  37. (*    - t is the tolerance for the precision of the solution                *)
  38. (*      praxis returns if the criterion                                     *)
  39. (*             2 * ||x(k)-x(k-1)|| <= sqrt(macheps) * ||x(k)|| + t          *)
  40. (*             is fulfilled more than ktm times, where                      *)
  41. (*    - ktm is another global parameter. it's default value is 1 and        *)
  42. (*      a value of 4 leads to a very(!) cautious stopping criterion         *)
  43. (*    - macheps is the relative machine precision and is                    *)
  44. (*      1.0 e-15 using an 8087 and turbo-87 and                             *)
  45. (*      1.0 e-06 without 8087.                                              *)
  46. (*    - h is a steplength parameter and shoul be set to the expected        *)
  47. (*      distance to the solution. an exceptional large or small value       *)
  48. (*      of h leads to slower 3 on the first few iterations.       *)
  49. (*    - scbd is a scaling parameter and should be set to about 10.          *)
  50. (*      the default is 1 and with that value no scaling is done at all      *)
  51. (*      it's only necessary when the parameters are scaled very different   *)
  52. (*    - illc is a boolean variable and should be set to true if the         *)
  53. (*      problem is known to be illconditioned.                              *)
  54. (****************************************************************************)
  55.  
  56.  
  57. function power(a,b: real):real;
  58. begin    power := exp(b*ln(a));
  59. end;     (* power *)
  60.  
  61.  
  62. procedure minfit(n:integer;eps,tol:real;var ab:matrix;var q:vector);
  63. label     1,
  64.           2,
  65.           3,
  66.           4;
  67. var       l, kt,
  68.           l2,
  69.           i, j, k: integer;
  70.           c, f, g,
  71.           h, s, x,
  72.           y, z:    real;
  73.           e:       vector;
  74. begin     (* householders reduction to bidiagonal form *)
  75.           x:= 0; g:= 0;
  76.           for i:= 1 to n do
  77.           begin e[i]:= g; s:= 0; l:= i+1;
  78.                 for j:= i to n do
  79.                     s:= s + sqr(ab[j,i]);
  80.                 if s<tol then g:= 0 else
  81.                 begin f:= ab[i,i];
  82.                       if f<0
  83.                          then g:= sqrt(s)
  84.                          else g:= - sqrt(s);
  85.                       h:= f*g-s;  ab[i,i]:= f - g;
  86.                       for j:= l to n do
  87.                       begin f:= 0;
  88.                             for k:= i to n do
  89.                                 f:= f + ab[k,i]*ab[k,j];
  90.                             f:= f/h;
  91.                             for k:= i to n do
  92.                                 ab[k,j]:= ab[k,j] + f*ab[k,i];
  93.                       end; (* j *)
  94.                 end; (* if *)
  95.                 q[i]:= g; s:= 0;
  96.                 if i<=n
  97.                    then for j:= l to n do
  98.                             s:= s + sqr(ab[i,j]);
  99.                 if s<tol then g:= 0 else
  100.                 begin f:= ab[i,i+1];
  101.                       if f<0
  102.                          then g:= sqrt(s)
  103.                          else g:= - sqrt(s);
  104.                       h:= f*g-s;  ab[i,i+1]:= f-g;
  105.                       for j:= l to n do e[j]:= ab[i,j]/h;
  106.                       for j:= l to n do
  107.                       begin s:= 0;
  108.                             for k:= l to n do s:= s + ab[j,k]*ab[i,k];
  109.                             for k:= l to n do ab[j,k]:= ab[j,k] + s*e[k];
  110.                       end; (* j *)
  111.                 end; (* if *)
  112.                 y:= abs(q[i])+abs(e[i]);
  113.                 if y > x then x:= y;
  114.           end; (* i *)
  115.           (* accumulation of right hand transformations *)
  116.           for i:= n downto 1 do
  117.           begin if g<>0.0 then
  118.                 begin h:= ab[i,i+1]*g;
  119.                       for j:= l to n do ab[j,i]:= ab[i,j]/h;
  120.                       for j:= l to n do
  121.                       begin s:= 0;
  122.                             for k:= l to n do s:= s + ab[i,k]*ab[k,j];
  123.                             for k:= l to n do ab[k,j]:= ab[k,j] + s*ab[k,i];
  124.                       end; (* j *)
  125.                 end; (* if *)
  126.                 for j:= l to n do
  127.                 begin ab[j,i]:= 0;
  128.                       ab[i,j]:= 0;
  129.                 end;
  130.                 ab[i,i]:= 1; g:= e[i]; l:= i;
  131.           end; (* i *)
  132.           (* diagonalization to bidiagonal form *)
  133.           eps:= eps*x;
  134.           for k:= n downto 1 do
  135.           begin kt:= 0;
  136. 1: kt:= kt + 1;
  137.                 if kt > 30 then
  138.                 begin e[k]:= 0;
  139.                       writeln('+++ qr failed');
  140.                 end;
  141.                 for l2:= k downto 1 do
  142.                 begin l:= l2;
  143.                       if abs(e[l])<=eps then goto 2;
  144.                       if abs(q[l-1])<=eps then goto 4;
  145.                 end; (* l2 *)
  146. 4:   c:= 0; s:= 1;
  147.                 for i:= l to k do
  148.                 begin f:= s*e[i]; e[i]:= c*e[i];
  149.                       if abs(f)<=eps then goto 2;
  150.                       g:= q[i];
  151.                       if abs(f) < abs(g)
  152.                          then h:= abs(g)*sqrt(1+sqr(f/g))
  153.                          else if f <> 0
  154.                                  then h:= abs(f)*sqrt(1+sqr(g/f))
  155.                                  else h:= 0;
  156.                       q[i]:= h;
  157.                       if h = 0 then
  158.                       begin h: