Ask Your Question

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

asked 2020-09-07 23:24:41 +0200

RolfF gravatar image

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 flag offensive close merge delete

1 Answer

Sort by ยป oldest newest most voted

answered 2020-09-08 10:19:52 +0200

Emmanuel Charpentier gravatar image

updated 2020-09-12 13:32:11 +0200

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")
# 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()]
# Extrema/saddle points
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.


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)
sage: f(0,y)

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...).


edit flag offensive delete link 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). (Plot of f(x,y) starting from (0,0)) (Plot of max_symbolic(0, f(x,y)) starting from (0,0))

RolfF gravatar imageRolfF ( 2020-09-08 11:09:56 +0200 )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.

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2020-09-09 13:47:27 +0200 )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,

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.

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2020-09-09 15:35:34 +0200 )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.

RolfF gravatar imageRolfF ( 2020-09-09 21:22:47 +0200 )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()
sage: J.variety()
verbose 0 (2274:, variety) Warning: falling back to very slow toy implementation.
[{x: 0, y: 0}]

(Result already given...)

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2020-09-10 08:41:13 +0200 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools


Asked: 2020-09-07 23:23:58 +0200

Seen: 617 times

Last updated: Sep 12 '20