Note, you may not have the ExtendGraphics package installed to be able to label contour plots. http://support.wolfram.com/applicationpacks/unsupported/extend/

<< ExtendGraphics`LabelContour`
<< Graphics`ImplicitPlot`
<< Graphics`FilledPlot`

f[x_, 
        y_]  =  (x + y - 1/4) Exp[-x^2 - 
            y^2]
g[x_, 
        y]  =  (x - 1/2)^2 + (y - 1/3)^2

    
Plot3D[f[x, y], {x, -2, 2}, {y, -3, 3}, 
      AxesLabel -> {x, y, z}, PlotPoints -> 35, 
      ColorFunction -> Hue, 
      PlotRange -> {-1/2, 1/2}]

    
plot  =  
      ContourPlot[f[x, y], {x, -1, 2}, {y, -1, 2}, 
        Contours -> 9, PlotPoints -> 40, FrameLabel -> {x, y}, 
        ContourShading -> False]
plot2  =  
      ImplicitPlot[g[x, y] == 1, {x, -1, 2}, 
        PlotStyle -> {Hue[h], 
            Dashing[{0.05, 0.05}]}]
plot3  =  LabelContourLines[plot]
Show[plot3, plot2]

Plot[{x, x^2/2 + 1/2, +x^8/8 + 7/8}, {x, 0, 2}, Frame -> True, 
      Axes -> False, FrameTicks -> False, 
      PlotStyle -> {{Hue[0.6]}, {Hue[0]}, {Hue[0]}}]

Simplify[D[g[x, y], x]]

Solve[g[x, y] == 1, x]

Eliminate[{D[f[x, y], x] == LM D[g[x, y], x], 
        D[f[x, y], y] == LM D[g[x, y], y], g[x, y] == 1}, 
      LM]
Solve[{%}, {x, y}];
N[%]
TeXForm[%]

f[0.107593, 1.25313]

f[x_, 
          y_]  =  (x + y - 1/4) Exp[-x^2 - 
              y^2];
g[x_, 
          y]  =  (x - 1/2)^2 + (y - 
                1/3)^2;
plot  =  
      ContourPlot[f[x, y], {x, -1, 3/2}, {y, -1, 2}, 
        Contours -> {-0.76, -0.70, -0.61, -0.52, 0.03, 
            0.12, 0.23, 0.31, 0.40}, PlotPoints -> 40, 
        FrameLabel -> {x, y}, 
        ContourShading -> False]
plot2  =  
      ImplicitPlot[g[x, y] == 1, {x, -1, 3/2}, 
        PlotStyle -> {Hue[h], Dashing[{0.05, 0.05}]}, 
        Frame -> 
          True, 
FrameLabel -> {x, 
            y}]
plot3  =  LabelContourLines[plot]
disk1  =  
      Graphics[Disk[{0.11, 1.25}, 0.04]]
disk2  =  
      Graphics[Disk[{-0.26, -0.31}, 0.04]]
Show[plot2, disk1, disk2, plot3, disk2, 
      PlotRange -> {-1.0, 1.5}, 
      AspectRatio -> 1]

f[x_, y_, z_, t_]  =  
      x^2 + y^2 + z^2 + t^2
Solve[{D[f[x, y, z, t], x] == 0, D[f[x, y, z, t], y] == 0, 
        D[f[x, y, z, t], z] == 0, D[f[x, y, z, t], t] == 0}, {x, 
        y, z, t}]

Clear[x, y, z, t, LM1, LM2]
f[x_, y_, z_, t_]  =  
      x^2 + y^2 + z^2 + t^2
g1[x_, y_, z_, t_]  =  x - z
g2[x_, y_, z_, t_]  =  y^2 + t

    
eq1  =  
      D[f[x, y, z, t], x] == 
        LM1 D[g1[x, y, z, t], x] + 
          LM2 D[g2[x, y, z, t], x]
eq2  =  
      D[f[x, y, z, t], y] == 
        LM1 D[g1[x, y, z, t], y] + 
          LM2 D[g2[x, y, z, t], y]
eq3  =  
      D[f[x, y, z, t], z] == 
        LM1 D[g1[x, y, z, t], z] + 
          LM2 D[g2[x, y, z, t], z]
eq4  =  
      D[f[x, y, z, t], t] == 
        LM1 D[g1[x, y, z, t], t] + 
          LM2 D[g2[x, y, z, t], t]
eq5  =  g1[x, y, z, t] == 2
eq6  =  g2[x, y, z, t] == 4
Eliminate[{eq1, eq2, eq3, eq4, eq5, eq6}, {LM1, 
        LM2}]

    
Solve[{D[f[x, y, z, t], x] == 
          LM1 D[g1[x, y, z, t], x] + 
            LM2 D[g2[x, y, z, t], x], 
D[f[x, y, z, t], 
            y] == 
          LM1 D[g1[x, y, z, t], y] + 
            LM2 D[g2[x, y, z, t], y], 
D[f[x, y, z, t], 
            z] == 
          LM1 D[g1[x, y, z, t], z] + 
            LM2 D[g2[x, y, z, t], z], 
D[f[x, y, z, t], 
            t] == 
          LM1 D[g1[x, y, z, t], t] + 
            LM2 D[g2[x, y, z, t], t], 
g1[x, y, z, 
            t] == 2, g2[x, y, z, t] == 4}, {x, y, z, t, LM1, 
        LM2}]

f[1, Sqrt[7/2], -1, 1/2]

N[%]

Here is how Mathematica can solve the TV problem in Chapter 2, both unconstrained and constrained.

Clear[x, y, f, LM1]
f[x_, 
        y_]  =  (339  -  0.01 x  - 0.003 y) x  +  (399  -  
              0.004 x  -  0.01 y) y  -  (400000 + 195  x + 
            225  y)

    
Solve[{D[f[x, y], x] == 0, D[f[x, y], y] == 0}, {x, 
        y}]

    
Clear[x, y, LM1]
g1[x_, y_]  =  x
Solve[{D[f[x, y], x] == LM1 D[g1[x, y], x], 
        D[f[x, y], y] == LM1 D[g1[x, y], y], 
        g1[x, y] == 5000}, {x, y, LM1}]

    
Clear[x, y, LM1]
g1[x_, y_]  =  y
Solve[{D[f[x, y], x] == LM1 D[g1[x, y], x], 
        D[f[x, y], y] == LM1 D[g1[x, y], y], 
        g1[x, y] == 8000}, {x, y, LM1}]

    
Clear[x, y, LM1]
g1[x_, y_]  =  x + y
Solve[{D[f[x, y], x] == LM1*D[g1[x, y], x], 
        D[f[x, y], y] == LM1*D[g1[x, y], y], 
        g1[x, y] == 10000}, {x, y, LM1}]

    
Clear[x, y, LM1]
g1[x_, y_]  =  x
Solve[{D[f[x, y], x] == LM1*D[g1[x, y], x], 
        D[f[x, y], y] == LM1*D[g1[x, y], y], g1[x, y] == 0}, {x, 
        y, LM1}]

    
Clear[x, y, LM1]
g1[x_, y_]  =  y
Solve[{D[f[x, y], x] == LM1*D[g1[x, y], x], 
        D[f[x, y], y] == LM1*D[g1[x, y], y], g1[x, y] == 0}, {x, 
        y, LM1}]

<< Miscellaneous`RealOnly`

f[x_, y_, 
        z_]  = ( 
            40  -  x/100 +  a y/100)  x  -  (45  -  
              5  x/100) y  + 3  z/100
g[x_, y_, z_]  =  4  x + y + z

    
D[f[x, y, z], x]
D[f[x, y, z], y]
D[f[x, y, z], z]

    
D[g[x, y, z], x]
D[g[x, y, z], y]
D[g[x, y, z], z]

    
sol  =  
        Solve[{D[f[x, y, z], x] == LM D[g[x, y, z], x], 
            D[f[x, y, z], y] == LM D[g[x, y, z], y], 
            D[f[x, y, z], z] == LM D[g[x, y, z], z], 
            g[x, y, z] == 10}, {x, y, z, 
            LM}];
Simplify[sol]
N[%]
TeXForm[sol]

z = ((2 (-39438 - 6962 a + 5 (a^2)))/(((5 + a))^2))
y = -((2 (5467 + 
                  1994 a))/(((5 + a))^2))
x = (4503/(5 + a))

Plot[x, {a, -10, 10}, 
      PlotRange -> {-10000, 10000}]
Plot[y, {a, -10, 10}]
Plot[z, {a, -10, 10}]
Plot[f[x, y, z], {a, -10, 10}]

Sensx[a_]  =  Simplify[D[x, a] a/x]
Sensy[a_]  =  Simplify[D[y, a] a/y]
Sensz[a_]  =  Simplify[D[z, a] a/z]
Sensf[a_]  =  
      Simplify[D[f[x, y, z], a] a/f[x, y, z]]
TeXForm[Sensf[α]]

Plot[Sensx[a], {a, -10, 10}, Frame -> True, 
      FrameLabel -> {"a", "Sensitivity of Optimal x"}, 
      PlotRange -> {-60, 50}]
Plot[Sensy[a], {a, -10, 10}, Frame -> True, 
      FrameLabel -> {"a", "Sensitivity of Optimal y"}, 
      PlotRange -> {-60, 50}]
Plot[Sensz[a], {a, -10, 10}, Frame -> True, 
      FrameLabel -> {"a", "Sensitivity of Optimal z"}, 
      PlotRange -> {-60, 50}]
Plot[Sensf[a], {a, -10, 10}, Frame -> True, 
      FrameLabel -> {"a", "Sensitivity of Maximum f"}, 
      PlotRange -> {-60, 50}]

Sensf[6. ]

Example 2.5

Here I create the plots for the examples in Chapter 2. Remember, you may not have the LabelContour package on your computer, so you won;t be able to label the contours.

I created this mainly for the last plot, to show what it looks like when you get a maximum at the intersection of two contraint curves.

Clear[x, y, z]

f[x_, 
        y_]  =  (339 - 0.01  x - 0.003  y) 
          x + (399 - 0.004  x - 0.01  y) 
          y - (400000 + 195  x + 225  y)
region[x_]  =  
      If[x < 5000, If[x < 2000, 8000, 10000 - x], 
        0]

constraint  =  
      FilledPlot[region[x], {x, 0, 8000}, 
        PlotRange -> {0, 
            12000}]

plot1  =  
      ContourPlot[f[x, y], {x, 0, 8000}, {y, 0, 12000}, 
        Contours -> {300000, 400000, 532308, 550000}, 
        PlotPoints -> 40, FrameLabel -> {x, y}, 
        ContourShading -> False]
plot3  =  
      LabelContourLines[plot1]

Show[constraint, plot3]

f[x_, 
        y_]  =  (339 - 0.01  x - 0.003  y) 
          x + (399 - 0.004  x - 0.01  y) 
          y - (400000 + 195  x + 225  y)
region[x_]  =  
      If[x < 3000, If[x < 2000, 8000, 10000 - x], 
        0]

constraint  =  
      FilledPlot[region[x], {x, 0, 8000}, 
        PlotRange -> {0, 
            12000}]

plot1  =  
      ContourPlot[f[x, y], {x, 0, 8000}, {y, 0, 12000}, 
        Contours -> {300000, 400000, f[3000, 7000], 550000}, 
        PlotPoints -> 40, FrameLabel -> {x, y}, 
        ContourShading -> False]
plot3  =  
      LabelContourLines[plot1]

Show[constraint, plot3]