# How to find the symbolic maximum of a 3D function?

I have a function of two real variables,

f(x,y) = a * (x * y - y^2) * exp( -b * x^2) * exp(b * x * y )

with the real parameters a and b, and five restrictions:

x>0; y>0; a>0; b>0; y<=x (y equal to or less than x)

In the 3d_plot, with a=b=1, i see there is only one local maximum, which is also the global maxium. How can i get the formula, for x and y coordinate of the maximum, as function of the parameters a and b?

I tried assume() for the restrictions, but that does not work and is made for the other way, to tell what is already there, not what should be changed.

To get rid of the negative function values, i used

r(x,y) = max_symbolic(0, f(x,y))

and that works, but that is only a workaround for the not implemented restrictions.

edit retag close merge delete

Sort by » oldest newest most voted

You have two types of conditions:

• $a>0$, $b>0$ : those are assumptions on the parameters, not part of a condition on candidate solutions

• $x>0$, $y>0$, $y<=x$ : those are restrictions imposed to candidate solutions.

Therefore :

# Problem-specific setup
var("x,y,a,b", domain="real")
f(x,y)=a*(x*y-y^2)*e^(-b*x^2)*e^(b*x*y)
assume(a>0,b>0)
# Extra conditions
Conds=[x>0, y>0, y<=x]
# Derivatives : must be all null to have an extremum/saddle point
Ders=[f(*f.arguments()).diff(u) for u in f.arguments()]
# Second derivatives : must be all negative to have a maximum
SDers=[f(*f.arguments()).diff(u,2) for u in f.arguments()]
Exts=solve(Ders,[x,y], to_poly_solve="force", solution_dict=True)
# Candidate solutions = maxima among extrema
CSols=[u for u in Exts if all([bool(s.subs(u)<=0) for s in SDers])]
# Solutions : maxima fulfilling extra conditions
Sols=[u for u in CSols if all([v.subs(u) for v in Conds])]

sage: Sols
[]


Note that

sage: CSols
[{x: 0, y: 0}]


The (only) maximum $f$ over $\mathbb{R}^2$ is excluded by your restrictions.

HTH,

EDIT : I've found an easier way yo show this. Consider your problem in polar coordinates, or, more lazily, by posing y=k*x (i. e. k=arctan(theta)).

sage: var("a, b, x, y, k", domain="real")
(a, b, x, y, k)
sage: f(x, y) = a * (x * y - y^2) * exp( -b * x^2) * exp(b * x * y )


We can get rid of the special cases of the axes by noting :

sage: f(x,0)
0
sage: f(0,y)
-a*y^2


The former is constant, the latter is maximal (and zero) in y=0.

Now, what is the maximum of f along the line of equation y=k*x ?

sage: Sx=f(x,k*x).diff(x).solve(x);Sx
[x == -sqrt(-1/(b*k - b)), x == sqrt(-1/(b*k - b)), x == 0]


(The third case has already been disposed of...). Now, what are the values of k maximizing f . Well, the derivatives of these maxima with respect to k :

sage: Dk=[f(x,k*x).subs(u).diff(k).simplify_full() for u in Sx[:2]];Dk
[a*e^(-k/(k - 1) + 1/(k - 1))/b, a*e^(-k/(k - 1) + 1/(k - 1))/b]


are constant :

sage: [u.log().expand_log().factor().exp().simplify_log() for u in Dk]
[a*e^(-1)/b, a*e^(-1)/b]


Therefore no bloody maximum, notwithstanding the numerical phantasms of Sage (or Mathematica, for that matter...).

HTH,

more

Thanks for the answer, but at (0, 0) is a saddle point, with value zero. In the 3D plot with a=1, b=1, i clearly see one maximum at appromimately (4, 3). I need that maximum, as function of a and b.

The reason why i need it, is that f(x,y) is the output power of a machine and i want to know the maximum power point (MPP).

http://www.true-random.com/files/unsorted/P_out_a_eq_b_eq_1_2020-09-06a.png (Plot of f(x,y) starting from (0,0)) http://www.true-random.com/files/unsorted/P_out_a_eq_b_eq_1_2020-09-06b.png (Plot of max_symbolic(0, f(x,y)) starting from (0,0))

( 2020-09-08 04:09:56 -0600 )edit

Youu are right ! Damn... The error is probably =Exts=solve(Ders,[x,y], to_poly_solve="force", solution_dict=True), which is extremely naive...

Working on it.

( 2020-09-09 06:47:27 -0600 )edit

Both Sage and Mathematica have a hard time finding simultaneous zeroes of $\frac{\partial f}{\partial x}$ and $\frac{\partial f}{\partial y}$. However,

f11(x,y)=f(x,y).subs({a:1,b:1})
implicit_plot(f11(x,y).diff(x), (-5,5), (-5,5))+implicit_plot(f11(x,y).diff(y), (-5,5), (-5,5), color="red")


is ... inspiring ...

RealLife (TM) interrupt. Back to it later.

( 2020-09-09 08:35:34 -0600 )edit

Thanks for the new answers and working on it!

Setting a=b=1 is only a workaround for plotting. a and b are machine parameters and independet from another.

( 2020-09-09 14:22:47 -0600 )edit

I think that the 3D-plots are misleading : what we "see" is a "maximum along the first diagonal ; furthermore, the position of this "maximum", as in:

plot3d(lambda x,y:max(0,f11(x,y)),(0,20),(0,20),aspect_ratio=[1,1,20],plot_points=100)


strongly depends on the resolution parameter. An artifact...

BTW, f(x,x)==0

Furthermore, the zeroes of the derivatives are given by :

sage: R0.<a,b>=QQ[]
sage: R1=FractionField(R0)
sage: R2.<x,y>=R1[]
sage: J=R2.ideal(*[(f(x,y).diff(u).factor()/(a*e^(-b*x^2 + b*x*y))).polynomial(ring=R2) for u in (x,y)])
sage: J.dimension()
0
sage: J.variety()
verbose 0 (2274: multi_polynomial_ideal.py, variety) Warning: falling back to very slow toy implementation.
[{x: 0, y: 0}]


( 2020-09-10 01:41:13 -0600 )edit