Parametric Find Maximum

Question

I'd like to plot the maximum of a curve in function of a parameter. Meaning I looking for maxima that depend on a parameter.

Here is what I tried to do but FindMaximum is not parametric, so it did not work :

a=1/4

F[x_, y_, z_] = (1 + a*y)*x^2*z^2 - b*x*y
mux[x, y, z] = D[F[x, y, z], x]
muy[x, y, z] = D[F[x, y, z], y]
muz[x, y, z] = D[F[x, y, z], z]
MUx[x, y] = (mux[x, y, z] - muz[x, y, z]) /. z -> 1 - x - y
MUy[x, y] = (muy[x, y, z] - muz[x, y, z]) /. z -> 1 - x - y
Hxy[x, y] = D[MUx[x, y], y]
Hyx[x, y] = D[MUy[x, y], x]
Hxy[x, y] - Hyx[x, y]
Hyy[x, y] = D[MUy[x, y], y]
Hxx[x, y] = D[MUx[x, y], x]

matrice = {{x*(1 + x)*Hxx[x, y] + Hxy[x, y]*x*y, 
   x*(1 + x)*Hxy[x, y] + Hyy[x, y]*x*y}, {y*(1 + y)*Hxy[x, y] + 
    Hxx[x, y]*x*y, y*(1 + y)*Hyy[x, y] + Hxy[x, y]*x*y}}

Eig = Eigenvalues[matrice]

f1[b_] = FindMaximum[{Eig[[1]], 
   0 < x < 1 && 0 < y < 1 && x + y < 1}, {x, y}][[1]]
f2[b_] = FindMaximum[{Eig[[2]], 
   0 < x < 1 && 0 < y < 1 && x + y < 1}, {x, y}][[1]]
Plot[f1[b], {b, 0, 10}]

I thought about using ParametricNdsolve and the derivative at the maximum would be 0, but nothing tells me that they are 0. Meaning, given that I'm looking for a maximum that is in some interval, it could be that the maximum is at the border, where the derivative is not zero.

How could I proceed plz ?

Answer 1

The problem is with the function definitions for f1 and f2. The use of Set (=) instead of SetDelayed (:=) means that FindMaximum is evaluated for the definition; however b does not have a numeric value, so FindMaximum fails. Then Part ([[1]]) extracts the first argument of FindMaximum and this becomes the function definition.

The big problem is that one cannot simply replace = by :=. Function definitions in Mathematica rewrite code. A definition of the form f1[b_] := <code> replaces the literal occurrences of b in <code> by the argument that matches the pattern b_. in Eig[[1]], b does not occur literally; instead, it is a global variable in the value of Eig. But Mathematica evaluates Eig after it has replaced the occurrences of b.

The simplest way around this is to use Block[] to temporarily assign b a value. (This is in effect what Table does.) We should also use ?NumericQ to prevent the functions from calling FindMaximum when b does not have a numeric value.

ClearAll[f1, f2];
f1[b0_?NumericQ] := Block[{b = b0}, 
   FindMaximum[{Eig[[1]], 0 < x < 1 && 0 < y < 1 && x + y < 1}, {x, y}][[1]]];
f2[b0_?NumericQ] := Block[{b = b0}, 
   FindMaximum[{Eig[[2]], 0 < x < 1 && 0 < y < 1 && x + y < 1}, {x, y}][[1]]];

For some reason, I get larger maxima than in @J.A's answer. Spot checks verify that the values below are indeed values of f1[b].

Plot[f1[b1], {b1, 0, 10}, MaxRecursion -> 1] // AbsoluteTiming

Answer 2

So I eventually used Table, which allows a "manual parametrization" :

a = 1/2
F[x_, y_, z_] = (1 + a*y)*x^2*z^2 - b*x*y
mux[x, y, z] = D[F[x, y, z], x]
muy[x, y, z] = D[F[x, y, z], y]
muz[x, y, z] = D[F[x, y, z], z]
MUx[x, y] = (mux[x, y, z] - muz[x, y, z]) /. z -> 1 - x - y
MUy[x, y] = (muy[x, y, z] - muz[x, y, z]) /. z -> 1 - x - y
Hxy[x, y] = D[MUx[x, y], y]
Hyx[x, y] = D[MUy[x, y], x]
Hxy[x, y] - Hyx[x, y]
Hyy[x, y] = D[MUy[x, y], y]
Hxx[x, y] = D[MUx[x, y], x]

matrice = {{x*(1 + x)*Hxx[x, y] + Hxy[x, y]*x*y, 
   x*(1 + x)*Hxy[x, y] + Hyy[x, y]*x*y}, {y*(1 + y)*Hxy[x, y] + 
    Hxx[x, y]*x*y, y*(1 + y)*Hyy[x, y] + Hxy[x, y]*x*y}}

Eig = Eigenvalues[matrice]

ListLinePlot[
 Table[FindMaximum[{Eig[[1]], 
     0 <= x <= 1 && 0 <= y <= 1 && x + y <= 1}, {x, y}][[1]], {b, 0,
    10, 0.1}]]
ListLinePlot[
 Table[FindMaximum[{Eig[[2]], 
     0 <= x <= 1 && 0 <= y <= 1 && x + y <= 1}, {x, y}][[1]], {b, 0,
    10, 0.1}]]

And I got the following graphs :