Chapter 13 Multiple Local Minima and Maxima

CHAPTER 13 13.13 Because of multiple local minima and maxima, there is no really simple means to test whether a single maximum occurs within an interval without actually performing a search. However, if we assume that the function has one maximum and no minima within the interval, a check can be included. Here is a VBA program to implement the Golden section search algorithm for maximization and solve Example 13.1. Option Explicit Sub GoldMax() Dim ier As Integer Dim xlow As Double, xhigh As Double Dim xopt As Double, fopt As Double xlow = 0 xhigh = 4 Call GoldMx(xlow, xhigh, xopt, fopt, ier) If ier = 0 Then MsgBox "xopt = " & xopt MsgBox "f(xopt) = " & fopt Else MsgBox "Does not appear to be maximum in [xl, xu]" End If End Sub Sub GoldMx(xlow, xhigh, xopt, fopt, ier) Dim iter As Integer, maxit As Integer, ea As Double, es As Double Dim xL As Double, xU As Double, d As Double, x1 As Double Dim x2 As Double, f1 As Double, f2 As Double Const R As Double = (5 ^ 0.5 - 1) / 2 ier = 0 maxit = 50 es = 0.001 xL = xlow xU = xhigh iter = 1 d = R * (xU - xL) x1 = xL + d x2 = xU - d f1 = f(x1) f2 = f(x2) If f1 > f2 Then xopt = x1 fopt = f1 Else xopt = x2 fopt = f2 End If If fopt > f(xL) And fopt > f(xU) Then Do d = R * d If f1 > f2 Then xL = x2 x2 = x1 x1 = xL + d f2 = f1 f1 = f(x1) Else xU = x1 x1 = x2 x2 = xU - d f1 = f2 f2 = f(x2) End If iter = iter + 1 If f1 > f2 Then xopt = x1 fopt = f1 Else xopt = x2 fopt = f2 End If If xopt <> 0 Then ea = (1 - R) * Abs((xU - xL) / xopt) * 100 If ea <= es Or iter >= maxit Then Exit Do Loop Else ier = 1 End If End Sub Function f(x) f = -(2 * Sin(x) - x ^ 2 / 10) End Function 13.14 The easiest way to set up a maximization algorithm so that it can do minimization is the realize that minimizing a function is the same as maximizing its negative. Therefore, the following algorithm minimizes or maximizes depending on the value of a user input variable, ind, where ind = -1 and 1 correspond to minimization and maximization, respectively. Option Explicit Sub GoldMinMax() Dim ind As Integer 'Minimization (ind = -1); Maximization (ind = 1) Dim xlow As Double, xhigh As Double Dim xopt As Double, fopt As Double xlow = 0 xhigh = 4 Call GoldMnMx(xlow, xhigh, -1, xopt, fopt) MsgBox "xopt = " & xopt MsgBox "f(xopt) = " & fopt End Sub Sub GoldMnMx(xlow, xhigh, ind, xopt, fopt) Dim iter As Integer, maxit As Integer, ea As Double, es As Double Dim xL As Double, xU As Double, d As Double, x1 As Double Dim x2 As Double, f1 As Double, f2 As Double Const R As Double = (5 ^ 0.5 - 1) / 2 maxit = 50 es = 0.001 xL = xlow xU = xhigh iter = 1 d = R * (xU - xL) x1 = xL + d x2 = xU - d f1 = f(ind, x1) f2 = f(ind, x2) If f1 > f2 Then xopt = x1 fopt = f1 Else xopt = x2 fopt = f2 End If Do d = R * d If f1 > f2 Then xL = x2 x2 = x1 x1 = xL + d f2 = f1 f1 = f(ind, x1) Else xU = x1 x1 = x2 x2 = xU - d f1 = f2 f2 = f(ind, x2) End If iter = iter + 1 If f1 > f2 Then xopt = x1 fopt = f1 Else xopt = x2 fopt = f2 End If If xopt <> 0 Then ea = (1 - R) * Abs((xU - xL) / xopt) * 100 If ea <= es Or iter >= maxit Then Exit Do Loop fopt = ind * fopt End Sub Function f(ind, x) f = ind * (1.1333 * x ^ 2 - 6.2667 * x + 1) End Function 13.15 Because of multiple local minima and maxima, there is no really simple means to test whether a single maximum occurs within an interval without actually performing a search. However, if we assume that the function has one maximum and no minima within the interval, a check can be included. Here is a VBA program to implement the Quadratic Interpolation algorithm for maximization and solve Example 13.2. Option Explicit Sub QuadMax() Dim ier As Integer Dim xlow As Double, xhigh As Double Dim xopt As Double, fopt As Double xlow = 0 xhigh = 4 Call QuadMx(xlow, xhigh, xopt, fopt, ier) If ier = 0 Then MsgBox "xopt = " & xopt MsgBox "f(xopt) = " & fopt Else MsgBox "Does not appear to be maximum in [xl, xu]" End If End Sub Sub Dim Dim Dim Dim QuadMx(xlow, xhigh, xopt, fopt, ier) iter As Integer, maxit As Integer, ea As Double, es As Double x0 As Double, x1 As Double, x2 As Double f0 As Double, f1 As Double, f2 As Double xoptOld As Double ier = 0 maxit = 50 es = 0.01 x0 = xlow x2 = xhigh x1 = (x0 + x2) / 2 f0 f1 f2 If = f(x0) = f(x1) = f(x2) f1 > f0 Or f1 > f2 Then xoptOld = x1 Do xopt = f0 * (x1^2 - x2^2) + f1 * (x2^2 - x0^2) + f2 * (x0^2 - x1^2) xopt = xopt / (2*f0 * (x1 - x2) + 2*f1 * (x2 - x0) + 2*f2 * (x0 - x1)) fopt = f(xopt) iter = iter + 1 If xopt > x1 Then x0 = x1 f0 = f1 x1 = xopt f1 = fopt Else x2 = x1 f2 = f1 x1 = xopt f1 = fopt End If If xopt <> 0 Then ea = Abs((xopt - xoptOld) / xopt) * 100 xoptOld = xopt If ea <= es Or iter >= maxit Then Exit Do Loop Else ier = 1 End If End Sub Function f(x) f = -(2 * Sin(x) - x ^ 2 / 10) End Function 13.16 Here is a VBA program to implement the Newton-Raphson method for maximization. Option Explicit Sub NRMax() Dim xguess As Double Dim xopt As Double, fopt As Double xguess = 2.5 Call NRMx(xguess, xopt, fopt) MsgBox "xopt = " & xopt MsgBox "f(xopt) = " & fopt End Sub Sub Dim Dim Dim Dim NRMx(xguess, xopt, fopt) iter As Integer, maxit As Integer, ea As Double, es As Double x0 As Double, x1 As Double, x2 As Double f0 As Double, f1 As Double, f2 As Double xoptOld As Double maxit = 50 es = 0.01 Do xopt = xguess - df(xguess) / d2f(xguess) fopt = f(xopt) If xopt <> 0 Then ea = Abs((xopt - xguess) / xopt) * 100 xguess = xopt If ea <= es Or iter >= maxit Then Exit Do Loop End Sub Function f(x) f = -(2 * Sin(x) - x ^ 2 / 10) End Function Function df(x) df = 2 * Cos(x) - x / 5 End Function Function d2f(x) d2f = -2 * Sin(x) - 1 / 5 End Function 13.17 Here is a VBA program to implement the Newton-Raphson method for maximization.  5 − 1 ( 4 − 2 ) = 1.23606 d1 =   2   x1 = 2 + d1 = 3.23606 x2 = 4 − d1 = 2.76394 f ( x1 ) = −4.69808 f ( x2 ) = −5.55333 f ( x2 ) < f ( x1 ) ⇒ x1 is new xU  5 −1  ( 3.23606 − 2) = 0.763927 d2 =   2    x1 = 2 + d 2 = 2.7639 x2 = 3.23606 − d 2 = 2.472133 f ( x1 ) = −5.55331 f ( x2 ) = −4.82656 f ( x1 ) < f ( x2 ) ⇒ x2 is new x L  5 −1  ( 3.23606 − 2.472133 ) = 0. 4721 d3 =   2    x1 = 2.472133 + d3 = 2.9442 x2 = 3.23606 − d3 = 2.7639 f ( x1 ) = −4.9353 f ( x2 ) = −5.55331 f ( x2 ) < f ( x1 ) ⇒ x1 is new xU  5 −1  ( 2.9442 − 2.472133 ) = 0.29175 d4 =   2    x1 = 2.472133 + d 4 = 2.7638 x2 = 2.9442 − d 4 = 2.6524 f ( x1 ) = −5.55331 f ( x2 ) = −5.4082 ∴ at time t = 2.76 , minimum pressure is –5.55331