| Download
Worksheet to illustrate a process for finding global maximum and minimum values for functions of two variables on compact domains
Project: Public Vector Calculus Worksheets
Views: 283# Finding maximum and minimum value of a function on a compact domain in 2D # Author: Thomas Krainer (Penn State Altoona) # # Version: 02/11/2016 # Changes from previous versions: # Syntax fixes for Cloud Version. Changed options for the display of function graph. # New assembly of the region plot to improve plot quality # # Version History: # # Version 11/15/2012 # Original Version: 08/17/2011 # # WARNING: Computations and 3D rendering take a LONG TIME until complete! # x,y,z,m=var('x,y,z,m') u,v,w=var('u,v,w') # Below: Functions that perform all computations def checkinequalities(functionlist,varsubs,operand='>',ndigits=4): """Takes a list of expressions 'functionlist' in the variables x and y, plugs in the values given \ in the dictionary 'varsubs' for x and y, uses the given operand (default is >)\ to compare all resulting expressions with zero. Returns true if varsubs contains parameters.\ Numerical accuracy is adjustable via digits (default is 4).""" if (len(SR(varsubs[x]).variables())+len(SR(varsubs[y]).variables()) > 0): return true else: string='('+str(round(RR(functionlist[0].subs({x:varsubs[x],y:varsubs[y]})),ndigits))+' '+operand+' 0)' for i in range(1,len(functionlist)): string+=' and ('+str(round(RR(functionlist[i].subs({x:varsubs[x],y:varsubs[y]})),ndigits))\ +' '+operand+' 0)' return sage_eval(string) def critpointsonboundary(f,functionlist,returndictionary=False): """Computes the solutions of the Lagrange multiplier equations for f(x,y) subject to all constraints in\ functionlist. Returns only the solutions for which all constraints are >= 0. Return value is by default\ a list giving all x, y, and m(ultiplier) values. Can change this to a dictionary.""" if (returndictionary == True): lagrange={} else: lagrange=[] for function in functionlist: F(x,y,m) = f - m*function # solution_dict=True does not support empty solution set --> reason for extra step through list. solnlist=solve([F.gradient()[0](x,y,m)==0,F.gradient()[1](x,y,m)==0,F.gradient()[2](x,y,m)==0],x,y,m) if (len(solnlist) != 0): solutions=solve([F.gradient()[0](x,y,m)==0,F.gradient()[1](x,y,m)==0,\ F.gradient()[2](x,y,m)==0],x,y,m,solution_dict=True) # Eliminate non-real point solutions. Keep parametrized pieces if returned by solve. tobedeleted=[] for solution in solutions: if (len(SR(solution[x]).variables())+len(SR(solution[y]).variables())\ +len(SR(solution[m]).variables())==0): if (CC(solution[x]).is_real() == False) or (CC(solution[y]).is_real() == False) \ or (CC(solution[m]).is_real() == False): tobedeleted.append(solution) for item in tobedeleted: solutions.remove(item) # Eliminate solutions that do not satisfy the inequalities g(x,y) >= 0 for all g in functionlist tobedeleted=[] for solution in solutions: if (checkinequalities(functionlist,solution,'>=') == False): tobedeleted.append(solution) for item in tobedeleted: solutions.remove(item) else: solutions=[] if (returndictionary == True): lagrange[function]=solutions else: lagrange.extend(solutions) return lagrange def getintersectionpoints(functionlist): """Computes a dictionary of all points (x,y) that satisfy f(x,y)=0 for precisely two elements f(x,y)\ in functionlist.""" if (len(functionlist)<=1): return [] else: g(x,y)=functionlist.pop() listofsolutions=getintersectionpoints(functionlist) for function in functionlist: solns=solve([g==0,function==0],x,y) if (len(solns) != 0): solutions=solve([g==0,function==0],x,y,solution_dict=True) else: solutions=[] # Avoid multiple appearances of the same point. Eliminate non-real solutions tobedeleted=[] for solution in solutions: if (solution in listofsolutions) or (CC(solution[x]).is_real() == False) or \ (CC(solution[y]).is_real() == False): tobedeleted.append(solution) for item in tobedeleted: solutions.remove(item) listofsolutions.extend(solutions) functionlist.append(g) return listofsolutions def intersectionpointsonboundary(functionlist): """Calls getintersectionpoints with functionlist and removes all points from the results of that call\ that do not satisfy g(x,y) >= 0 for all g(x,y) in functionlist.""" allintersections=getintersectionpoints(functionlist) tobedeleted=[] for point in allintersections: if (checkinequalities(functionlist,point,'>=') == False): tobedeleted.append(point) for item in tobedeleted: allintersections.remove(item) return allintersections def critpointsinterior(f,functionlist): """Finds the critical points (x_c,y_c) of f(x,y) that satisfy g(x_c,y_c) > 0 for all g in functionlist,\ and returns a list of dictionaries of these values.""" h(x,y)=f tobereturned=[] solnlist=solve([h.gradient()[0]==0,h.gradient()[1]==0],x,y) # solution_dict=True does not support empty solution set --> reason for extra step through list. if (len(solnlist) != 0): solutions=solve([h.gradient()[0]==0,h.gradient()[1]==0],x,y,solution_dict=True) # Eliminate non-real point solutions. Keep parametrized pieces if returned by solve. tobedeleted=[] for solution in solutions: if (len(SR(solution[x]).variables())+len(SR(solution[y]).variables()) == 0): if (CC(solution[x]).is_real() == False) or (CC(solution[y]).is_real() == False): tobedeleted.append(solution) for item in tobedeleted: solutions.remove(item) # Eliminate solutions that do not satisfy the inequalities g(x,y) > 0 for all g in functionlist tobedeleted=[] for solution in solutions: if (checkinequalities(functionlist,solution,'>') == False): tobedeleted.append(solution) for item in tobedeleted: solutions.remove(item) tobereturned.extend(solutions) return tobereturned def mcleanup(valuelist): """Gets valuelist of dictionaries {x: #, y: #, m: #) and removes the m-parts.""" tobereturned=[] if (len(valuelist) > 0): for value in valuelist: tobereturned.append({x:value[x], y:value[y]}) return tobereturned def listunion(lists): """Takes a list of lists 'lists', and gives back a list of the union (without multiple occurences).""" tobereturned=[] for list in lists: for item in list: if (item in tobereturned): pass else: tobereturned.append(item) return tobereturned def getminmaxfromlist(f,xylist): """Plugs points from xylist into f, and returns triple dictionary (x,y,z) of min/max values\ and occurences. Only one occurence is returned.""" zlist=[] for point in xylist: zlist.append(RR(f.subs({x:point[x], y:point[y]}))) minimum=SR(f.subs({x: xylist[zlist.index(min(zlist))][x], y: xylist[zlist.index(min(zlist))][y]})) maximum=SR(f.subs({x: xylist[zlist.index(max(zlist))][x], y: xylist[zlist.index(max(zlist))][y]})) return ([{x: xylist[zlist.index(min(zlist))][x], y: xylist[zlist.index(min(zlist))][y],z: minimum},{x: xylist[zlist.index(max(zlist))][x], y: xylist[zlist.index(max(zlist))][y],z: maximum}]) def getminmaxondomain(f,functionlist,returnxylist=False): """Returns maximum and minimum of f(x,y) on domain defined by g(x,y) >= 0 for all g in functionlist.\ Output is a tuple of two dictionaries {x: #, y: #, z: #}, where (x,y) marks one occurence of the minimum\ and maximum value, respectively.""" list1=critpointsinterior(f,functionlist) auxlist_21=critpointsonboundary(f,functionlist,returndictionary=True) auxlist_22=[] for function in functionlist: auxlist_22.append(auxlist_21[function]) auxlist_23=listunion(auxlist_22) list2=mcleanup(auxlist_23) # list2=mcleanup(critpointsonboundary(f,functionlist)) list3=intersectionpointsonboundary(functionlist) xylist=listunion([list1,list2,list3]) if (returnxylist==False): return getminmaxfromlist(f,xylist) else: return [getminmaxfromlist(f,xylist),xylist,list1,auxlist_21,list3] # Below: Functions to format plots def viewingbox(functionlist,r=0.1): """Determines boundingbox data [[xmin,xmax],[ymin,ymax]] that contains the domain specified by\ g(x,y) >= 0 for all g in functionlist. r specifies by how much beyond the domain we wish to extend\ the viewing window in x- and y-directions.""" minmax_x = getminmaxondomain(x,functionlist) xgap=minmax_x[1][z]-minmax_x[0][z] minmax_y = getminmaxondomain(y,functionlist) ygap=minmax_y[1][z]-minmax_y[0][z] return [[minmax_x[0][z]-r*xgap,minmax_x[1][z]+r*xgap],[minmax_y[0][z]-r*ygap,minmax_y[1][z]+r*ygap]] def boundingbox(f,functionlist,r=0.1,extremevalues=[]): """Gives back the bounding box as a list [xmin,xmax,ymin,ymax,zmin,zmax].""" [[xmin,xmax],[ymin,ymax]]=viewingbox(functionlist,r) if (len(extremevalues) == 0): extremevalues=getminmaxondomain(f,functionlist) zgap=RR(extremevalues[1][z])-RR(extremevalues[0][z]) zmin=RR(extremevalues[0][z])-r*zgap zmax=RR(extremevalues[1][z])+r*zgap return [RR(xmin),RR(xmax),RR(ymin),RR(ymax),RR(zmin),RR(zmax)] preamble='<h1>Finding maximum and minimum values in 2D</h1>' preamble+='<p>This program illustrates the process of finding the maximum and minimum values of functions ' preamble+='of two variables $f(x,y)$ on closed<br>and bounded domains $D$ in ${\\mathbb R}^2$. ' preamble+='The domain $D$ is specified by a set of defining functions $g_1(x,y),\\ldots,g_N(x,y)$ in ' preamble+='the sense that<br>$(x,y) \\in D$ if and only if $g_j(x,y) \\geq 0$ for all $j = 1,\\ldots,N$.</p>' preamble+='<p>It is assumed that all functions are continuously differentiable, that $D$ is closed and ' preamble+='bounded, and that for each defining function $g(x,y)$<br>we have $\\nabla g \\neq 0$ whenever $g = 0$, ' preamble+='and at points where more than one of the defining function vanishes the respective gradients<br>' preamble+='are supposed to be independent.<br>' preamble+='The interior of $D$ is then the set of points where all functions $g_j$ are strictly positive, and ' preamble+='the boundary consists of pieces of zero sets<br>of the defining functions. The condition on the ' preamble+='gradients of the defining functions insures that the boundary is a collection of<br>' preamble+='smooth curves that intersect each other transversally.</p>' preamble+='<p>In the program below, you can specify the function $f(x,y)$ to be discussed and a list ' preamble+='$[g_1(x,y),\\ldots,g_N(x,y)]$ to specify the domain.<br>Depending on the input, it may take quite a while ' preamble+='for the program to finish. Proper functionality strongly depends on the symbolic<br>solver. ' preamble+='It is advisable to stick to algebraic functions of simple structure, and to avoid cases ' preamble+='where the critical point equations <br>' preamble+='or Lagrange multiplier equations are solved by parametrized ' preamble+='curves rather than merely a finite collection of points.' html(preamble) @interact def MaxMinonDomains(f=input_box(default=(1/5)*((x-1/2)^2+(y-3)^2)*((x+1)^2+(y+2)^2)/(x^2+y^2+1),label='$f(x,y) =$',type=SR,width=150),\ functionlist=input_box(default=[4-x*y,4+x*y,16-x^2-y^2,(x-1/2)^2+(4/9)*(y-7/10)^2-1,(1/4)*(x+4)^2-y], label='Defining functions: ',type=list,width=150),\ show_computations=('Show computations:',True),\ show_plots=('Show plots:',True),\ rescale_regionplot=('Rescale region plot:',False),\ rescale_graphplot=('Rescale graph plot:',True)): # Suggested default values: TAKES A LONG TIME! # # functionlist=[4-x*y,4+x*y,16-x^2-y^2,(x-1/2)^2+(4/9)*(y-7/10)^2-1,(1/4)*(x+4)^2-y] # f(x,y)=(1/5)*((x-1/2)^2+(y-3)^2)*((x+1)^2+(y+2)^2)/(x^2+y^2+1) # Make sure input is properly formatted auxfunction(x,y)=f f=auxfunction for i in range(0,len(functionlist)): auxfunction(x,y)=functionlist[i] functionlist[i]=auxfunction # Get data about max/min of f and critical points, etc. maxminf_full=getminmaxondomain(f,functionlist,returnxylist=True) maxminf=maxminf_full[0] xylist=maxminf_full[1] interiorpoints=maxminf_full[2] boundarypoints=maxminf_full[3] boundaryintersections=maxminf_full[4] # Assemble written output output='<h3>Function and domain</h3>' output+='Function: ' output+='$$ f(x,y) = %s. $$ <br>'%(latex(SR(f(x,y).simplify_full()))) output+='Defining functions for the domain: ' output+='$$ \\left\\{\\begin{aligned} ' output+='g_{1}(x,y) &= %s '%(latex(SR(functionlist[0](x,y).simplify_full()))) for i in range(1,len(functionlist)): output+=' \\\ g_{%s}(x,y) &= %s '%(latex(i+1),latex(SR(functionlist[i](x,y).simplify_full()))) output+='\\end{aligned}\\right. $$' output1='<h3>Critical points in the interior</h3>' output1+='To find the critical points in the interior we have to solve the following equations: ' output1+='$$ \\left\\{\\begin{aligned} \\nabla f(x,y) &= \\langle 0,0 \\rangle \\\ ' output1+='g(x,y) &> 0 \\; \\textrm{for all defining functions} \\; g(x,y).\\end{aligned}\\right. $$ <br>' output1+='For this example this amounts to solving the following: ' output1+='$$ \\left\\{\\begin{aligned} %s &= 0 \\\ %s &= 0 '%(latex(SR(f.gradient()[0].simplify_full())),latex(SR(f.gradient()[1].simplify_full()))) for function in functionlist: output1+='\\\ %s &> 0 '%(latex(SR(function.simplify_full()))) output1+='\\end{aligned}\\right. $$ <br>' output1+='The solutions are all points in the following set ($\\{\\}$ is the empty set): ' output1+='$$ (x,y) \in \\left\\{' if (len(interiorpoints)==0): output1+=' \\right\\} $$ ' else: output1+='\\left(%s,%s\\right)'%(latex(SR(interiorpoints[0][x]).simplify_full()),\ latex(SR(interiorpoints[0][y]).simplify_full())) for i in range(1,len(interiorpoints)): output1+=',\\, \\left(%s,%s\\right)'%(latex(SR(interiorpoints[i][x]).simplify_full()),\ latex(SR(interiorpoints[i][y]).simplify_full())) output1+=' \\right\\} $$ ' output1+='<h3>Critical points on the boundary</h3>' output1+='For every defining function $g(x,y)$ we solve the Lagrange multiplier equations and ' output1+='retain those <br>solutions $(x,y)$ that are inside the domain, i.e., $h(x,y) \\geq 0$ for all ' output1+='other defining functions $h(x,y)$: ' output1+='$$ \\left\\{\\begin{aligned} \\nabla f(x,y) &= m\\cdot\\nabla g(x,y) \\\ ' output1+='g(x,y) &= 0 \\\ h(x,y) &\\geq 0 \\; \\textrm{for all other defining functions} \\; h(x,y). ' output1+='\\end{aligned}\\right. $$ <br>' for function in functionlist: output1+='Defining function: $g(x,y) = %s$ '%(latex(SR(function).simplify_full())) output1+='$$ \\left\\{\\begin{aligned} %s &= %s \\\ %s &= %s \\\ %s &= 0 '%(latex(SR(f.gradient()[0].simplify_full())),latex(m*SR(function.gradient()[0].simplify_full())),latex(SR(f.gradient()[1].simplify_full())),latex(m*SR(function.gradient()[1].simplify_full())),latex(SR(function).simplify_full())) auxlist=[] for auxfunction in functionlist: if (auxfunction == function): pass else: auxlist.append(auxfunction) for auxfunction in auxlist: output1+='\\\ %s &\\geq 0 '%(latex(SR(auxfunction.simplify_full()))) output1+='\\end{aligned}\\right. $$ <br>' output1+='Solutions: $(x,y;m) \\in \\left\\{' if (len(boundarypoints[function]) == 0): output1+=' \\right\\}$ <br><br>' else: output1+='\\left(%s,%s;%s\\right)'%(latex(SR(boundarypoints[function][0][x]).simplify_full()),\ latex(SR(boundarypoints[function][0][y]).simplify_full()),\ latex(SR(boundarypoints[function][0][m]).simplify_full())) for i in range(1,len(boundarypoints[function])): output1+=',\\, \\left(%s,%s;%s\\right)'%(latex(SR(boundarypoints[function][i][x]).simplify_full()),\ latex(SR(boundarypoints[function][i][y]).simplify_full()),\ latex(SR(boundarypoints[function][i][m]).simplify_full())) output1+=' \\right\\}$ <br><br>' output1+='<h3>Intersection points of boundary curves</h3>' output1+='Maximum and minimum values could be attained at the intersection points of any two ' output1+='boundary curves. <br>For all pairs of defining functions $g(x,y)$ and $h(x,y)$ we need ' output1+='to find the solutions $(x,y)$ <br>to the following equations: ' output1+='$$ \\left\\{\\begin{aligned} g(x,y) &= 0 \\\ h(x,y) &= 0 \\\ k(x,y) &\\geq 0 ' output1+='\\; \\textrm{for all other defining functions} \\; k(x,y). ' output1+='\\end{aligned}\\right. $$ <br>' output1+='For this example this amounts to solving %s systems of equations. '%(latex(binomial(len(functionlist),2))) output1+='<br><br>Solutions: $(x,y) \\in \\left\\{' if (len(boundaryintersections)==0): output1+=' \\right\\}$ <br><br>' else: output1+='\\left(%s,%s\\right)'%(latex(SR(boundaryintersections[0][x]).simplify_full()),\ latex(SR(boundaryintersections[0][y]).simplify_full())) for i in range(1,len(boundaryintersections)): output1+=',\\, \\left(%s,%s\\right)'%(latex(SR(boundaryintersections[i][x]).simplify_full()),\ latex(SR(boundaryintersections[i][y]).simplify_full())) output1+=' \\right\\}$ <br><br>' if (show_computations == True): output+=output1 output+='<h3>Evaluation of the function to determine maximum and minimum</h3>' output+='<table border="1" cellspacing=""5"" cellpadding=""5"">' output+='<tr><th> $(x,y)$ </th><th> $f(x,y)$ </th><th> Max/Min </th></tr>' for point in xylist: if (f(point[x],point[y]) == maxminf[0][z]): output+='<tr><td> $\\left(%s,%s\\right)$ </td><td> $%s$ </td><td> Minimum </td></tr>'%(latex(SR(point[x].simplify_full())),latex(SR(point[y].simplify_full())),latex(SR(f(point[x],point[y]).simplify_full()))) else: if (f(point[x],point[y]) == maxminf[1][z]): output+='<tr><td> $\\left(%s,%s\\right)$ </td><td> $%s$ </td><td> Maximum </td></tr>'%(latex(SR(point[x].simplify_full())),latex(SR(point[y].simplify_full())),latex(SR(f(point[x],point[y]).simplify_full()))) else: output+='<tr><td> $\\left(%s,%s\\right)$ </td><td> $%s$ </td><td></td></tr>'%(latex(SR(point[x].simplify_full())),latex(SR(point[y].simplify_full())),latex(SR(f(point[x],point[y]).simplify_full()))) output+='</table>' html(output) # Assemble plots [xmin,xmax,ymin,ymax,zmin,zmax] = boundingbox(f,functionlist,extremevalues=maxminf) Graph=implicit_plot3d(w-f(u,v),(u,xmin,xmax),(v,ymin,ymax),(w,zmin,zmax),plot_points=150,color=(135/255,180/255,250/255),opacity=0.9,region = lambda u,v,w : checkinequalities(functionlist,{x:u, y:v},'>=')) regionstring='['+str(SR(functionlist[0]))+' >= 0' for counter in range(1,len(functionlist)): regionstring+=', '+str(SR(functionlist[counter]))+' >= 0' regionstring+=']' regionsage = sage_eval(regionstring,{'x':x,'y':y}) Region=region_plot(regionsage,(x,xmin,xmax),(y,ymin,ymax),incol='lightblue',bordercol='black',plot_points=150) Regionbox=Region.get_minmax_data() Regionsize=vector((Regionbox['xmax']-Regionbox['xmin'],Regionbox['ymax']-Regionbox['ymin'])).norm() PointsOfInterest=Graphics() for point in xylist: if (len(SR(point[x]).variables())+len(SR(point[y]).variables()) == 0): if (f(point[x],point[y]) == maxminf[0][z]): PointsOfInterest+=circle((point[x],point[y]),0.007*Regionsize,fill=True,\ color='green',edgecolor='green') Graph+=point3d((point[x],point[y],f(point[x],point[y])),size=10,color='green') else: if (f(point[x],point[y]) == maxminf[1][z]): PointsOfInterest+=circle((point[x],point[y]),0.007*Regionsize,fill=True,\ color='red',edgecolor='red') Graph+=point3d((point[x],point[y],f(point[x],point[y])),size=10,color='red') else: PointsOfInterest+=circle((point[x],point[y]),0.007*Regionsize,fill=True,\ color='blue',edgecolor='blue') Graph+=point3d((point[x],point[y],f(point[x],point[y])),size=10,color='blue') if (show_plots == True): if (rescale_regionplot == False): show(Region+PointsOfInterest,aspect_ratio=1) else: show(Region+PointsOfInterest) if (rescale_graphplot == False): show(Graph,aspect_ratio=1) else: show(Graph)
Finding maximum and minimum values in 2D
This program illustrates the process of finding the maximum and minimum values of functions of two variables f(x,y) on closed
and bounded domains D in R2. The domain D is specified by a set of defining functions g1(x,y),…,gN(x,y) in the sense that
(x,y)∈D if and only if gj(x,y)≥0 for all j=1,…,N.
It is assumed that all functions are continuously differentiable, that D is closed and bounded, and that for each defining function g(x,y)
we have ∇g=0 whenever g=0, and at points where more than one of the defining function vanishes the respective gradients
are supposed to be independent.
The interior of D is then the set of points where all functions gj are strictly positive, and the boundary consists of pieces of zero sets
of the defining functions. The condition on the gradients of the defining functions insures that the boundary is a collection of
smooth curves that intersect each other transversally.
In the program below, you can specify the function f(x,y) to be discussed and a list [g1(x,y),…,gN(x,y)] to specify the domain.
Depending on the input, it may take quite a while for the program to finish. Proper functionality strongly depends on the symbolic
solver. It is advisable to stick to algebraic functions of simple structure, and to avoid cases where the critical point equations
or Lagrange multiplier equations are solved by parametrized curves rather than merely a finite collection of points.
Interact: please open in CoCalc