## Local maximum and minimum for functions with two variables

Definition. Let f(x,y) defined over a region R which contains the point (a,b).

1. Function f(x,y) has local maximum in (a,b), if f(a,b) \geq f(x,y) for all points (x,y) in an open disk with center (a,b).

2. Function f(x,y) has local minimum (a,b), if f(a,b) \leq f(x,y) for all points (x,y) in an open disk with center (a,b).

We can use the following theorem for finding them.

Theorem 1. If f(x,y) has local max or min in the interior of R where the first partial derivatives exist, then f_{x}(a,b) = 0 and f_{y}(a,b) = 0.

We call these points critical. We say that a critical point is an inflexion point, if in a small disk with center (a,b), there are points with f(x,y) > f(a,b) and f(x,y) < f(a,b).

In order to find out what kind of point is a critical point (i.e. local max or min) we use the following criterion.

Criterion. Let f(x,y) has second order partial derivatives and are continuous in a disk with center (a,b) and that f_{x}(a,b) = f_{y}(a,b) = 0. Then:

i. f has local maximum in (a,b) if f_{xx} < 0 and f_{xx}f_{yy} - f_{xy}^{2} > 0 at (a,b).

ii. f has local maximum in (a,b) if f_{xx} > 0 and f_{xx}f_{yy} - f_{xy}^{2} > 0 at (a,b).

iii. f has inflexion point in (a,b) if f_{xx}f_{yy} - f_{xy}^{2} < 0 at (a,b).

iv. The criterion does not work if f_{xx}f_{yy} - f_{xy}^{2} = 0 at (a,b).

We can easily implement the previous in Sage (we did not include the case where we have infinitely many points).

"""
Find local maxima and minima of functions with two variables in R^2

AUTHORS:

Argyris Kalampakas, [email protected] : initial version 25/5/2015

This program is free software: you can redistribute it and/or modify
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>

"""

#############
# We work at R^2 so we don't have boundary points.
# f : symbolic expression with two variables <<x,y>>
# flag : default=true, enables plotting
# TODO : fmin : default=0, for growing the plain of the function, accepts only positive
# TODO : symbolic expression become variable independent
#############

def local_maxmin(f, flag = true):

print "The function f is: "
print f

# Algorithm for Fermat's theorem

### Algorithm for keeping only the real critical points
type(f)
fx = diff(f,x)
fy = diff(f,y)
criticalpoint_ = solve([fx, fy], x, y)

def are_real_roots(s):
criticalpoint=[]
for i in range(0,len(s)):
if s[i][0].rhs() in RR and s[i][1].rhs() in RR:
criticalpoint.append(criticalpoint_[i])
return criticalpoint
criticalpoint=are_real_roots(criticalpoint_)
###

if criticalpoint:

print "--------------------"
print "The critical points are:"
print [[cr[0], cr[1]] for cr in criticalpoint]

Min = [[cr[0].rhs()] for cr in criticalpoint] + [[cr[1].rhs()] for cr in criticalpoint] + [[f(cr[0].rhs(),cr[1].rhs())] for cr in criticalpoint]
Min = min(Min) + max(Min)
Min = abs(Min[0]),abs(Min[1])
Min = min(Min)

# Change the range of the plain
F = plot3d(f,(x,-Min-2,Min+2),(y,-Min-2,Min+2))

# Algorithm for second derivative test

fxx = diff(fx,x)
fyy = diff(fy,y)
fxy = diff(fx,y)
h = fxx*fyy - fxy**2
if h.has(x) & h.has(y):
H = [[h(cr[0].rhs(), cr[1].rhs())] for cr in criticalpoint]
elif h.has(x):
H = [[h(cr[0].rhs())] for cr in criticalpoint]
elif h.has(y):
H = [[h(cr[1].rhs())] for cr in criticalpoint]
else:
H = h

LM=[]
if type(H) is list:

for i in range(0,len(H)):
i = 0

if H[i][0] > 0:

if fxx.has(x):

if fxx(criticalpoint[i][0].rhs()) > 0:
LM.append(point3d((criticalpoint[i][0].rhs(), criticalpoint[i][1].rhs(),
f(criticalpoint[i][0].rhs(),criticalpoint[i][1].rhs())),color='black', size=15))
print "--------------------"
print criticalpoint[i]
print "It is local minimum, the value of the function at this point is:"
print f(criticalpoint[i][0].rhs(),criticalpoint[i][1].rhs())
elif fxx(criticalpoint[i][0].rhs()) < 0:
LM.append(point3d((criticalpoint[i][0].rhs(), criticalpoint[i][1].rhs(),
f(criticalpoint[i][0].rhs(),criticalpoint[i][1].rhs())),color='red', size=15))
print "--------------------"
print criticalpoint[i]
print "It is local maximum, the value of the function at this point is:"
print f(criticalpoint[i][0].rhs(),criticalpoint[i][1].rhs())
else:
LM.append(point3d((criticalpoint[i][0].rhs(), criticalpoint[i][1].rhs(),
f(criticalpoint[i][0].rhs(),criticalpoint[i][1].rhs())),color='yellow', size=15))
print "--------------------"
print criticalpoint[i]
print "Hessian is zero. We do not have a method for this point, to examine what it is."
## We have to create higher-order derivative test

criticalpoint.pop(H.index(H[i]))
H.pop(H.index(H[i]))

else:

if fxx > 0:
LM.append(point3d((criticalpoint[i][0].rhs(), criticalpoint[i][1].rhs(),
f(criticalpoint[i][0].rhs(),criticalpoint[i][1].rhs())),color='black', size=15))
print "--------------------"
print criticalpoint[i]
print "It is local minimum, the value of the function at this point is:"
print f(criticalpoint[i][0].rhs(),criticalpoint[i][1].rhs())
elif fxx < 0:
LM.append(point3d((criticalpoint[i][0].rhs(), criticalpoint[i][1].rhs(),
f(criticalpoint[i][0].rhs(),criticalpoint[i][1].rhs())),color='red', size=15))
print "--------------------"
print criticalpoint[i]
print "It is local maximum, the value of the function at this point is:"
print f(criticalpoint[i][0].rhs(),criticalpoint[i][1].rhs())
else:
LM.append(point3d((criticalpoint[i][0].rhs(), criticalpoint[i][1].rhs(),
f(criticalpoint[i][0].rhs(),criticalpoint[i][1].rhs())),color='yellow', size=15))
print "--------------------"
print criticalpoint[i]
print "Hessian is zero. We do not have a method for this point, to examine what it is."
## We have to create higher-order derivative test

criticalpoint.pop(H.index(H[i]))
H.pop(H.index(H[i]))

elif H[i][0] < 0:
LM.append(point3d((criticalpoint[i][0].rhs(), criticalpoint[i][1].rhs(),
f(criticalpoint[i][0].rhs(),criticalpoint[i][1].rhs())),color='green', size=15))
print "--------------------"
print criticalpoint.pop(H.index(H[i]))
H.pop(H.index(H[i]))
else:
LM.append(point3d((criticalpoint[i][0].rhs(), criticalpoint[i][1].rhs(),
f(criticalpoint[i][0].rhs(),criticalpoint[i][1].rhs())),color='yellow', size=15))
print "--------------------"
print criticalpoint.pop(H.index(H[i]))
print "Hessian is zero. We do not have a method for this point, to examine what it is."
H.pop(H.index(H[i]))

else:
if H > 0:

if fxx.has(x):

if fxx(criticalpoint[0][0].rhs()) > 0:
LM.append(point3d((criticalpoint[0][0].rhs(), criticalpoint[0][1].rhs(),
f(criticalpoint[0][0].rhs(),criticalpoint[0][1].rhs())),color='black', size=15))
print "--------------------"
print criticalpoint[0]
print "It is local minimum, the value of the function at this point is:"
print f(criticalpoint[0][0].rhs(),criticalpoint[0][1].rhs())
elif fxx(criticalpoint[0][0].rhs()) < 0:
LM.append(point3d((criticalpoint[0][0].rhs(), criticalpoint[0][1].rhs(),
f(criticalpoint[0][0].rhs(),criticalpoint[0][1].rhs())),color='red', size=15))
print "--------------------"
print criticalpoint[0]
print "It is local maximum, the value of the function at this point is:"
print f(criticalpoint[0][0].rhs(),criticalpoint[0][1].rhs())
else:
LM.append(point3d((criticalpoint[0][0].rhs(), criticalpoint[0][1].rhs(),
f(criticalpoint[0][0].rhs(),criticalpoint[0][1].rhs())),color='yellow', size=15))
print "--------------------"
print criticalpoint[0]
print "Hessian is zero. We do not have a method for this point, to examine what it is."
## We have to create higher-order derivative test

else:

if fxx > 0:
LM.append(point3d((criticalpoint[0][0].rhs(), criticalpoint[0][1].rhs(),
f(criticalpoint[0][0].rhs(),criticalpoint[0][1].rhs())),color='black', size=15))
print "--------------------"
print criticalpoint[0]
print "It is local minimum, the value of the function at this point is:"
print f(criticalpoint[0][0].rhs(),criticalpoint[0][1].rhs())
elif fxx < 0:
LM.append(point3d((criticalpoint[0][0].rhs(), criticalpoint[0][1].rhs(),
f(criticalpoint[0][0].rhs(),criticalpoint[0][1].rhs())),color='red', size=15))
print "--------------------"
print criticalpoint[0]
print "It is local maximum, the value of the function at this point is:"
print f(criticalpoint[0][0].rhs(),criticalpoint[0][1].rhs())
else:
LM.append(point3d((criticalpoint[0][0].rhs(), criticalpoint[0][1].rhs(),
f(criticalpoint[0][0].rhs(),criticalpoint[0][1].rhs())),color='yellow', size=15))
print "--------------------"
print criticalpoint[0]
print "Hessian is zero. We do not have a method for this point, to examine what it is."
## We have to create higher-order derivative test

elif H < 0:
LM.append(point3d((criticalpoint[0][0].rhs(), criticalpoint[0][1].rhs(),
f(criticalpoint[0][0].rhs(),criticalpoint[0][1].rhs())),color='green', size=15))
print "--------------------"
print criticalpoint.pop()
else:
LM.append(point3d((criticalpoint[0][0].rhs(), criticalpoint[0][1].rhs(),
f(criticalpoint[0][0].rhs(),criticalpoint[0][1].rhs())),color='yellow', size=15))
print "--------------------"
print criticalpoint.pop()
print "Hessian is zero. We do not have a method for this point, to examine what it is."

if flag:
print "--------------------"
print "Green points are saddle points."
print "Red points are local maximum."
print "Black points are local minimum."
print "Yellow points are for those that we do not have a method to examine what it is."
PLOT = F
while LM:
PLOT = PLOT + LM.pop(0)

show(PLOT)

else:
print "--------------------"
print "This function either has infinite many critical points or does not have any. In the first case we can not compute them with our code."

return None

# Define the function
sage:var('x,y')
sage:f = x**3 + y**3 - x
sage:local_maxmin(f,false)

The function f is:
x^3 + y^3 - x
--------------------
The critical points are:
[[x == -1/3*sqrt(3), y == 0], [x == 1/3*sqrt(3), y == 0]]
--------------------
[x == -1/3*sqrt(3), y == 0]
Hessian is zero. We do not have a method for this point, to examine what it is.
--------------------
[x == 1/3*sqrt(3), y == 0]
Hessian is zero. We do not have a method for this point, to examine what it is.


## Finding extrema points of functions subject to some constraints (Lagrange multipliers)

Differential_Calculus_of_many_variables (last edited 2015-07-08 11:19:04 by drazioti)