Numerical Methods: Mathematica 3.0 Notebooks
(c) John H. Mathews, 1998
To accompany the text:
Numerical Methods:  for Mathematics, 
Science, and Engineering, 2nd Ed, 1992
Prentice Hall, Inc.
Simon & Schuster,  One Lake Street
Upper Saddle River, New Jersey 07458  U.S.A.
Prentice Hall, Inc.; USA, Canada, Mexico ISBN 0-13-624990-6
Prentice Hall, International Editions:      ISBN 0-13-625047-5
www.Prenhall.com
This free software is compliments of the author.
E-mail address:      in%"mathews@fullerton.edu"

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.

Clear["Global`*"];
Off[General::"spell1"];
Clear[a, a0, b, b0, c, count, Cv, d, delta, epsilon, Ev, f, f0, 
        GoldenSearch, GradSearch, h, Hi, Ho, j, k, Li, Lo, max, min, M, n, 
        Nelder, p, QuadraticSearch, rone, rtwo, R, sum, V0, x0, X0, Y, Ya, 
        Yb, Yc, Yd, Yp, YC, YE, YR]On[General::"spell1"];

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

Execute this cell to activate Golden Search for a Minimum.
GoldenSearch[f0_,a0_,b0_,delta_,epsilon_]
f0_ is the function f(x),
a0_ is the left endpoint,
b0_ is the right endpoint,
delta_ is the tolerance for the abscissa,
epsilon_ is the tolerance for the ordinate.

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 = N[1\/2\ \((\@5 - 1)\)]; 
		rtwo = rone\^2; 
		h = b - a; 
        
		Ya = N[f[a]]; 
		Print[{N[a, 14], N[Ya, 14]}]; 
		
        Yb = N[f[b]]; 
		Print[{N[b, 14], N[Yb, 14]}]; 
		
        c = a + rtwo\ h; 
		d = a + rone\ h; 
		Yc = N[f[c]]; 
		
        Print[{N[c, 14], N[Yc, 14]}]; 
		Yd = N[f[d]]; 
		
        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}, PlotStyle -> RGBColor[1, 0, 0]]; 

    Print["y = f[x] = ", f[x]]; \)

GoldenSearch[f[x], 0, 1, 0.0000000001, 0.0000000001];

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

\(sol = FindMinimum[f[x], {x, 0.5}]; 
Print["f[x] = ", f[x]]; 

    Print["  p  = ", 
      N[sol[[2, 1, 2]], 14]]; 

    Print["f[p] = ", 
      N[sol[[1]], 14]]; \)

Algorithm 8.2  Nelder-Mead's Minimization Method.  
To approximate a local minimum of  f ((Subscript[x, 1]), (Subscript[x, 2]), ⋯, (Subscript[x, n])),  where  f(X)  is a continuous function of  n  real variables, and given the  n + 1  initial starting points  (Subscript[V, k]) = ((Subscript[v, k, 1]), (Subscript[v, k, 2]), ⋯, (Subscript[v, k, n]))  for  k = 0, 1, ⋯, n.
Section  8.1   Minimization of a Function   Page  414

Execute this cell to activate Nelder-Mead's Minimization Method.
Nelder[V0_,epsilon_]
V0_ is the starting point,
epsilon_ is the tolerance.

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 = \@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]] = 
                  1\/2\ 
                    \((V[[j + 1]] + 
                        V[[Lo + 1]])\); 
                
						
                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 = 2\ R - 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 = 1\/2\ 
                    \((V[[Hi + 1]] + M)
                      \); 
						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]], 12], 
              N[V[[Lo + 1]], 12]}]; ]; 

        		size; 
		Print[" "]; 
		
        Print["f[{x,y}] = ", f[{x, y}]]; 
		
        Print["A local minimum of f[{x,y}] is:"]; 
		
        Print["  P  = ", 
          N[V[[Lo + 1]], 12]]; 
		
        Print["F[P] = ", 
          N[Y[[Lo + 1]], 12]]; ]; 

    \(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; 

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

\(V = {{0.9, 0.9}, {0.9, 1.1}, {1.1, 0.9}}; 
Nelder[V, 0.0000000001]; 
    \)

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

\(sol = FindMinimum[f[{x, y}], {x, 1.8}, {y, 1.2}]; 

    Print["f[{x,y}] = ", f[{x, y}]]; 
    Print["  P  = ", 
      N[{sol[[2, 1, 2]], 
          sol[[2, 2, 2]]}, 12]]; 
    Print["F[P] = ", 
      N[sol[[1]], 12]]; \)

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  (Subscript[p, 0])  and then searching the intervals  [a, (Subscript[p, 0])],  and  [(Subscript[p, 0]), b].
Section  8.1   Minimization of a Function   Page  416

Execute this cell to activate Local Minimum Search Using Quadratic Interpolation.
QuadraticSearch[f0_,x0_,delta_,epsilon_]
f0_ is the function f(x),
x0_ is the starting value,
delta_ is the tolerance for the abscissa,
epsilon_ is the tolerance for the ordinate.

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 = N[(10^9)];
		k = 0;
		Err = 1.0;
		cond = 0;
		h = 1.0;
		If[Abs[p0] > 10000.0, h = N[Abs[p0]/10000]];
		While[k < Kmax && Err > epsilon && cond ≠ 5,
			If[(N > 0)[f'[p0]],
				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 = N[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;
			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) = 3 (x^2) - 2 x + 5  over  [0, 1].

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

    Plot[f[x], {x, 0, 1}, PlotStyle -> RGBColor[1, 0, 0]]; 

    Print["y = f[x] = ", f[x]]; \)

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}]; 
Print["f[x] = ", f[x]]; 

    Print["  p  = ", 
      N[sol[[2, 1, 2]], 14]]; 

    Print["f[p] = ", 
      N[sol[[1]], 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 = ((Subscript[x, 1]), (Subscript[x, 2]), ⋯, (Subscript[x, n]))  by starting with one point  (Subscript[P, 0])  and using the gradient method.
Section   8.1   Minimization of a Function   Page  418

Execute this cell to activate Steepest Descent  or Gradient Method. 
GradSearch[X0_,delta_,epsilon_]
X0_ is the starting point,
delta_ is the tolerance for the point,
epsilon_ is the tolerance for the function value.

Off[General::"spell1"]; Clear[GradSearch, X0, delta, epsilon]; \(GradSearch[X0_, delta_, epsilon_] := 
	
      Module[{}, 
		
        Qmin := 
			
          Block[{}, 
				big = N[10\^8]; 
				jmax = 15; 
				
            P0 = N[P0]; 
				P1 = N[P0 + h\ S0]; 
				
            P2 = N[P0 + 2\ h\ S0]; 
				Y0 = N[f[P0]]; 
				
            Print[{N[Y0, 14], N[P0, 14]}]; 
				Y1 = N[f[P1]]; 
				
            Print[{N[Y1, 14], N[P1, 14]}]; 
				Y2 = N[f[P2]]; 
				
            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[N[length > 10\^4], 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, 12], N[Pmin, 12]}]; 
				count = count + 1]]; 
        
		Print[" "]; 
		Print["f[{x,y}] = ", f[{x, y}]]; 

        		Print["A local minimum of f[{x,y}] is:"]; 
		
        Print["  P  = ", N[Pmin, 12]]; 
		
        Print["F[P] = ", N[Ymin, 12]]; ]; 
\(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; 

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

\(grad = {∂\_x\ f[{x, y}], ∂\_y\ f[{x, y}]}; 

    g[{x_, y_}] = grad; 
norm[V_] := N[\@\(V\^2 /. List -> Plus\)]; 

    \(s[{x_, y_}] = \(-\(g[{x, y}]\/norm[grad]\)\); \)\)\(P0 = {0.99, 1.01}; 
GradSearch[P0, 0.0000000001, 0.0000000001]; \)

Let us compare this answer with Mathematica's routine.

\(sol = FindMinimum[f[{x, y}], {x, 1.8}, {y, 1.2}]; 

    Print["f[{x,y}] = ", f[{x, y}]]; 
    Print["  P  = ", 
      N[{sol[[2, 1, 2]], 
          sol[[2, 2, 2]]}, 14]]; 
    Print["F[P] = ", 
      N[sol[[1]], 14]]; \)