Numerical Methods: Mathematica Notebooks
(c) John H. Mathews, 1996

To accompany the text:

Numerical Methods:  for Mathematics, 
Science, and Engineering, 2nd Ed, 1992
Prentice Hall, Inc.
Englewood Cliffs, New Jersey 07632  U.S.A.
Prentice Hall, Inc.; USA, Canada, Mexico ISBN 0-13-624990-6
Prentice Hall, International Editions:      ISBN 0-13-625047-5

This free software is compliments of the author.
E-mail address:  in%"mathews@fullerton.edu"
http://titan.fullerton.edu/~mathews/

CONTENTS

Chapter 8.	Numerical Optimization
	Algorithm 8.1	Golden Search for a Minimum	
	Algorithm 8.2	Nelder-Mead's Minimization Method	
	Algorithm 8.3	Local Minimum Search Using Quadratic Interpolation	
	Algorithm 8.4	Steepest Descent or Gradient Method

Algorithm 8.1  (Golden Search for a Minimum).  
To numerically approximate the minimum of f(x) on the interval  [a,b]  by using a golden search.  Proceed with the method only if f(x) is a unimodal function on the interval  [a,b].
Section	8.1	Minimization of a Function   Page  413

Off[General::spell1];
Clear[GoldenSearch,f0,a0,b0,delta,epsilon,p,Yp];
GoldenSearch[f0_,a0_,b0_,delta_,epsilon_] :=
Module[{f,a,b,c,d,h,k,rone,rtwo,Ya,Yb,Yc,Yd},
  Set@@{f[x_],f0};
  a = a0;
  b = b0;
  rone = (Sqrt[5]-1)/2 //N;
  rtwo = rone^2;
  h = b - a;
  Ya = f[a] //N; Print[{N[a,14],N[Ya,14]}];
  Yb = f[b] //N; Print[{N[b,14],N[Yb,14]}];
  c = a + rtwo h;
  d = a + rone h;
  Yc = f[c] //N; Print[{N[c,14],N[Yc,14]}];
  Yd = f[d] //N; Print[{N[d,14],N[Yd,14]}];
  k = 1;
  While[Abs[Yb-Ya]>epsilon || h>delta,
    If[Yc < Yd,
      Module[{},
        b = d;
        Yb = Yd;
        d = c;
        Yd = Yc;
        h = b - a;
        c = a + rtwo h;
        Yc = f[c];
        Print[{N[c,14],N[Yc,14]}]; ],
      Module[{},
        a = c;
        Ya = Yc;
        c = d;
        Yc = Yd;
        h = b - a;
        d = a + rone h;
        Yd = f[d];
        Print[{N[d,14],N[Yd,14]}]; ]; ];
    k = k+1 ];
  p = a;
  Yp = Ya;
  If[Yb < Ya, p = b; Yp = Yb];
  Print[" "];
  Print["The function is f[x] = ",f[x]];
  Print["A local minimum of f[x] is:"]; 
  Print["  p  = ",N[p,14]];
  Print["f[p] = ",N[Yp,14]]; ];
On[General::spell1];

Chapter 8, Example 8.2, Page 404.

Find the minimum of the function   f(x)  =  x^2  -  sin(x)  over  [0,1].

Clear[f,x];
f[x_] = x^2 - Sin[x]

Plot[f[x],{x,0,1}];

GoldenSearch[f[x],0,1,10^-10,10^-15]

Let us compare this answer with Mathematica's subroutine FindMinimum.

sol = FindMinimum[f[x],{x,0.5}];
N[sol,14]

Algorithm 8.2  (Nelder-Mead's Minimization Method).  
To approximate a local minimum of f(x1,x2,...,xn), where f is a continuous
function of n real variables, and given the n+1 initial starting points  
Vk = (vk,1,...,vk,n)  for  k = 0, 1, ...,n.
Section	8.1	Minimization of a Function   Page  414

Remark. Subroutines need not have parameters, they can be a name only.

Off[General::spell1];
Clear[Nelder,V0,Y,epsilon,n,Lo,Hi,Li,Ho,
      Cv,Ev,M,R,YC,YE,YR,max,min,count];
Nelder[V0_,epsilon_] :=
Module[{n},
 size := Module[{j,k,sum},
  norm = 0;
  For[j=0,j<=n,j++,
   sum = 0;
   For[k=1,k<=n,k++,
    sum = sum + (V[[Lo+1,k]]-V[[j+1,k]])^2];
    If[sum > norm, norm = sum]; ];        
  norm = Sqrt[norm];
  Return[norm]; ];
 order := Module[{j},
  Lo = 0;
  Hi = 0;
  For[j=1,j<=n,j++,
   If[Y[[j+1]]<Y[[Lo+1]], Lo = j];
   If[Y[[j+1]]>Y[[Hi+1]], Hi = j]; ];
  Li = Hi;
  Ho = Lo;
  For[j=0,j<=n,j++,
   If[j!=Lo && Y[[j+1]]<Y[[Li+1]], Li = j];
   If[j!=Hi && Y[[j+1]]>Y[[Ho+1]], Ho = j]; ]; ];
 newpoints := Module[{j,sum},
  sum = -V[[Hi+1]];
  For[j=0,j<=n,j++,
   sum = sum + V[[j+1]]; ];
  M = sum/n;
  R = 2 M - V[[Hi+1]];
  YR = f[R] ];
 shrink := Module[{j},
  For[j=0,j<=n,j++,
   If[j+1 != Lo+1,
    V[[j+1]] = (V[[j+1]]+V[[Lo+1]])/2;
    Y[[j+1]] = f[V[[j+1]]]; ]; ]; ];
 replace := Module[{},
  V[[Hi+1]] = R; Y[[Hi+1]] = YR ];
 improve := Module[{},
  If[YR < Y[[Ho+1]],
   Module[{},
    If[Y[[Li+1]] < YR,
     replace,
     Module[{},
      Ev = 2R - M;
      YE = f[Ev];
      If[YE < Y[[Li+1]],
       V[[Hi+1]] = Ev;
       Y[[Hi+1]] = YE,
       replace ]; ]; ]; ],
   Module[{},
    If[YR < Y[[Hi+1]], replace ];
    Cv = (V[[Hi+1]] + M)/2;
    YC = f[Cv];
    If[YC < Y[[Hi+1]],
     V[[Hi+1]] = Cv; Y[[Hi+1]] = YC,
     shrink ]; ]; ]; ];
 V = V0;
 n = Length[V[[1]]];
 Y = Table[0,{n+1}];
 For[k=0,k<=n,k++,
  Y[[k+1]] = f[V[[k+1]]]; ];
 min = 2;
 max = 55;
 count = 0;
 order; 
 While[(Y[[Hi+1]]>Y[[Lo+1]]+epsilon && count < max) ||
  count < min,
  newpoints;
  improve;
  count = count+1;
  order; 
  Print[{N[Y[[Lo+1]],14],N[V[[Lo+1]],14] }]; ];
 size;
 Print[" "];
 Print["f[{x,y}] = ",f[{x,y}]];
 Print["A local minimum of f[{x,y}] is:"]; 
 Print["  P  = ",N[V[[Lo+1]],14]];
 Print["F[P] = ",N[Y[[Lo+1]],14]]; ];
On[General::spell1];

Chapter 9, Example 8.4, Page 409.

Use the Nelder-Mead method to find the minimum of  

f(x,y)  =  100(y - x^2)^2 + (1 - x)^2,

Rosenbrock's parabolic valley, circa 1960.

Clear[f,x,y,V];
f[{x_,y_}] := 100(y - x^2)^2 + (1 - x)^2
f[{x,y}]
V = {{0.9,0.9},
     {0.9,1.1},
     {1.1,0.9}}

Plot3D[f[{x,y}],{x,0.5,1.5},{y,0.5,1.5},
       AxesLabel->{"x","y","z"},
       ViewPoint->{1,3,1}];

Nelder[V,10^-10]

Let us compare this answer with Mathematica's subroutine FindMinimum.

sol = FindMinimum[f[{x,y}],{x,1.8},{y,1.2}];
N[sol,14]

Algorithm 8.3 (Local Minimum Search Using Quadratic
Interpolation).  To find a local minimum of the function  f(x)  over the interval  [a,b],  by starting with one initial approximation  p0  and then 
searching the intervals  [a, p0]  and  [p0, b].
Section	8.1	Minimization of a Function   Page  416

Off[General::spell1];
Clear[QuadraticSearch,f0,x0,delta,epsilon];
QuadraticSearch[f0_,x0_,delta_,epsilon_] :=
Module[{f},
 Set@@{f[x_],f0};
 p0 = x0;
 Jmax = 50;
 Kmax = 50;
 big = 10^9 //N;
 k = 0;
 Err = 1.0;
 cond = 0;
 h = 1.0;
 If[Abs[p0] > 10000.0, h = Abs[p0]/10000 //N];
 While[k<Kmax && Err>epsilon && cond != 5, 
    If[f'[p0]//N>0, h=-Abs[h]];
    p1 = p0 + h;
    p2 = p0 + 2 h;
    Pmin = p0;
    Y0 = f[p0]; Print[{N[p0,14],N[Y0,14]}];
    Y1 = f[p1]; Print[{N[p1,14],N[Y1,14]}];
    Y2 = f[p2]; Print[{N[p2,14],N[Y2,14]}];
    Ymin = Y0;
    cond = 0;
    j = 0;
    While[j<Jmax && Abs[h]>delta && cond==0,
      If[Y0 <= Y1,
        Module[{},
          p2 = p1;
          Y2 = Y1;
          h = h/2;
          p1 = p0 + h;
          Y1 = f[p1];
          Print[{N[p1,14],N[Y1,14]}]; ],
        Module[{},
          If[Y2 < Y1,
            Module[{},
              p1 = p2;
              Y1 = Y2;
              h = 2 h;
              p2 = p0 + 2*h;
              Y2 = f[p2];
              Print[{N[p2,14],N[Y2,14]}]; ],
            cond = -1 ] ] ];
      j = j+1;
      If[Abs[h]>big || Abs[P0]>big, cond=5] ];
    If[cond==5,
      Pmin = p1;
      Ymin = f[p1];
      Print[{N[Pmin,14],N[Ymin,14]}];
      Goto[getout] ];
    d = 4 Y1 - 2 Y0 - 2 Y2 //N;
    If[d < 0,
      hmin = h (4 Y1 - 3 Y0 - Y2)/d,
      hmin = h/3; cond=4 ];
    Pmin = p0 + hmin;
    Ymin = f[Pmin]; Print[{N[Pmin,14],N[Ymin,14]}];
    h  = Abs[h];
    h0 = Abs[hmin];  
    h1 = Abs[hmin - h];
    h2 = Abs[hmin - 2 h];
    If[h0 < h, h = h0];
    If[h1 < h, h = h1];
    If[h2 < h, h = h2];
    If[h == 0, h = hmin];
    If[h < delta, cond=1];
    If[Abs[h]>big || Abs[Pmin]>big, cond=5];
    E0 = Abs[Y0 - Ymin];
    E1 = Abs[Y1 - Ymin];
    E2 = Abs[Y2 - Ymin];
    If[E0!=0 && E0<Err, Err=E0];
    If[E1!=0 && E1<Err, Err=E1];
    If[E2!=0 && E2<Err, Err=E2];
    If[E0==0 && E1==0 && E2==0, Err=0];
    If[Err<epsilon, cond=2];
    p0 = Pmin;
    k = k+1;
    Label[getout] ];
  If[cond == 2 && h < delta, cond=3];   
  Print[" "];
  Print["The function is f[x] = ",f[x]];
  Print["A local minimum of f[x] is:"]; 
  Print["  p  = ",N[Pmin,14]];
  Print["f[p] = ",N[Ymin,14]]; ];
On[General::spell1];

Chapter 8, Exercise 4, Page 420.

Find the minimum of the function   f(x)  =  3x^2  -  2x  +  5  over  [0,1].

Clear[f,x];
f[x_] = 3x^2 - 2x + 5

Plot[f[x],{x,0,1}];

Remark: For some reason, the first execution of the subroutine takes a long time to get past the fourth iterate.  However, it runs much faster the second time it is executed!  I can't explain it.  Maybe Mathematica learns?

QuadraticSearch[f[x],0.5,0.001,0.001]

Let us compare this answer with Mathematica's routine.

sol = FindMinimum[f[x],{x,0.5}];
N[sol,14]

Algorithm 8.4  (Steepest Descent  or Gradient Method).
To numerically approximate a local minimum of  f(X),  where  f  is 
a continuous function of  n  real variables and  X = (x1, x2, ..., xn)  
by starting with one point  P0  and using the gradient method.
Section	8.1	Minimization of a Function   Page  418

Off[General::spell1];
Clear[GradSearch,X0,delta,epsilon];
GradSearch[X0_,delta_,epsilon_] :=
Module[{},
  Qmin := Block[{},
  big = 10^8 //N;
  jmax = 15;
  P0 = P0 //N;
  P1 = P0 + h S0 //N;
  P2 = P0 + 2 h S0 //N;
  Y0 = f[P0] //N; Print[{N[Y0,14],N[P0,14]}];
  Y1 = f[P1] //N; Print[{N[Y1,14],N[P1,14]}];
  Y2 = f[P2] //N; Print[{N[Y2,14],N[P2,14]}];
  cond = 0;
  j = 0;
  While[j < jmax && cond == 0,
    length = norm[P0];
    If[Y0 <= Y1,
      Module[{},
        P2 = P1;
        Y2 = Y1;
        h  = h/2;
        P1 = P0 + h S0;
        Y1 = f[P1]; ],
      Module[{},
        If[Y2 < Y1,
          Module[{},
            P1 = P2;
            Y1 = Y2;
            h  = 2 h;
            P2 = P0 + 2 h S0;
            Y2 = f[P2]; ],
          cond = -1] ] ];
      j = j+1;
      If[h < delta, cond = 1];
      If[Abs[h]>big || length>big, cond=5] ];
    count = count + j;
    If[cond==5, Pmin=P1; Ymin=Y1; Goto[getout] ];
    d = 4 Y1 - 2 Y0 - 2 Y2;
    If[d < 0,
      hmin = h (4 Y1 - 3 Y0 - Y2)/d,
      hmin = h/3; cond = 4 ];
    Pmin = P0 + hmin S0;
    Ymin = f[Pmin];
    h0 = Abs[hmin];
    h1 = Abs[hmin - h];
    h2 = Abs[hmin - 2 h];
    If[h0 < h, h = h0];
    If[h1 < h, h = h1];
    If[h2 < h, h = h2];
    If[h == 0, h = hmin];
    If[h  < delta, cond = 1];
    E0 = Abs[Y0 - Ymin];
    E1 = Abs[Y1 - Ymin];
    E2 = Abs[Y2 - Ymin];
    If[E0!=0 && E0<Err, Err = E0];
    If[E1!=0 && E1<Err, Err = E1];
    If[E2!=0 && E2<Err, Err = E2];
    If[E0==0 && E1==0 && E2==0, Err = 0];
    If[Err < epsilon, cond = 2];
    If[cond==2 && h < delta, cond = 3];
    Label[getout] ];
  P0 = X0;
  max = 60;
  h = 1;
  length = norm[P0];
  If[length > 10^4//N, h = length/10^4];
  Err = 1;
  count = 0;
  cond  = 0;
  While[count < max && cond != 5 &&
    (h > delta || Err > epsilon),
  Module[{},
    S0 = s[P0];
    Qmin;
    P0 = Pmin;
    Y0 = Ymin;
     Print[{N[Ymin,14],N[Pmin,14]}];
    count = count + 1 ] ];
  Print[" "];
  Print["f[{x,y}] = ",f[{x,y}]];
  Print["A local minimum of f[{x,y}] is:"]; 
  Print["  P  = ",N[Pmin,14]];
  Print["F[P] = ",N[Ymin,14]]; ];
On[General::spell1];

Chapter 9, Example 8.4, Page 409.

Use the gradient method to find the minimum of  

f(x,y)  =  100(y - x^2)^2 + (1 - x)^2,

Rosenbrock's parabolic valley, circa 1960.

Clear[f,x,y,grad,g,s,norm,V,P0];
f[{x_,y_}] =  100(y - x^2)^2 + (1 - x)^2
grad = Map[D[f[{x,y}],#]&,{x,y}];
g[{x_,y_}] = grad
norm[V_] := Sqrt[V^2/.List->Plus]//N;
s[{x_,y_}] = - g[{x,y}]/norm[grad];
P0 = {0.99,1.01}

Plot3D[f[{x,y}],{x,0.5,1.5},{y,0.5,1.5},
       AxesLabel->{"x","y","z"},
       ViewPoint->{1,3,1}];

GradSearch[P0,10^-10,10^-10]

Let us compare this answer with Mathematica's routine.

sol = FindMinimum[f[{x,y}],{x,1.8},{y,1.2}];
N[sol,14]