222 - Topic 09 - Multivariate Extrema -- Sage

Extrema for Multivariate Functions

John Travis

Mississippi College

This sage worksheet takes a given multivariate function and compute various derivatives, critical points and the other stuff necessary to locate critical points and aid in describing any possible extreme values. Students using this for preparing homework should copy the results of each cell into a word processor and prepare a report with all conclusions. In other words, let Sage do the calculus and the student do the interpretation/application. If pasting formulas into a Word document, the formulas might look funny. That's ok.

Don't forget to get a nice perspective on your surface with all the points plotted (at the bottom) before pasting it into your document.

%auto # To start things off, enter the function f that you want to use. var('x,y') # f = x^2 - 3*x*y - y^2 # f = x^2*y+y^3-12*y # f = 2*x^3-3*x^2*y+3*x^2-y^3+3*y-4 # f = y^3-3*y*x^2-3*y^2-3*x^2+1 # f = x*sin(y) # this one has infinitely many critical points so the solver has fits # f = (1/2-x^3+y^2)*e^(1-x^2-y^2) f = (12*x*y-3/2*x^2*y^2)/(2*x+2*y) # this one has issues with d show(f)

$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{3 \, {\left(x^{2} y^{2} - 8 \, x y\right)}}{4 \, {\left(x + y\right)}}$

(Be certain to evaluate each cell--in order--from the top. The results from each active cell should print below the cell.)

fx = derivative(f,x) fy = derivative(f,y) grad = [fx,fy] print 'The gradient vector' show(grad)

The gradient vector

$\newcommand{\Bold}[1]{\mathbf{#1}}\left[-\frac{3 \, {\left(x y^{2} - 4 \, y\right)}}{2 \, {\left(x + y\right)}} + \frac{3 \, {\left(x^{2} y^{2} - 8 \, x y\right)}}{4 \, {\left(x + y\right)}^{2}}, -\frac{3 \, {\left(x^{2} y - 4 \, x\right)}}{2 \, {\left(x + y\right)}} + \frac{3 \, {\left(x^{2} y^{2} - 8 \, x y\right)}}{4 \, {\left(x + y\right)}^{2}}\right]$

The gradient vector

fxx = diff(fx,x) fyy = diff(fy,y) fxy = diff(fx,y) print 'The collection of second partial derivatives' show([fxx,fyy,fxy])

The collection of second partial derivatives

$\newcommand{\Bold}[1]{\mathbf{#1}}\left[-\frac{3 \, y^{2}}{2 \, {\left(x + y\right)}} + \frac{3 \, {\left(x y^{2} - 4 \, y\right)}}{{\left(x + y\right)}^{2}} - \frac{3 \, {\left(x^{2} y^{2} - 8 \, x y\right)}}{2 \, {\left(x + y\right)}^{3}}, -\frac{3 \, x^{2}}{2 \, {\left(x + y\right)}} + \frac{3 \, {\left(x^{2} y - 4 \, x\right)}}{{\left(x + y\right)}^{2}} - \frac{3 \, {\left(x^{2} y^{2} - 8 \, x y\right)}}{2 \, {\left(x + y\right)}^{3}}, -\frac{3 \, {\left(x y - 2\right)}}{x + y} + \frac{3 \, {\left(x^{2} y - 4 \, x\right)}}{2 \, {\left(x + y\right)}^{2}} + \frac{3 \, {\left(x y^{2} - 4 \, y\right)}}{2 \, {\left(x + y\right)}^{2}} - \frac{3 \, {\left(x^{2} y^{2} - 8 \, x y\right)}}{2 \, {\left(x + y\right)}^{3}}\right]$

The collection of second partial derivatives

fyx = diff(fy,x) show(fyx-fxy)

$\newcommand{\Bold}[1]{\mathbf{#1}}0$

solns = solve(grad,x,y,solution_dict=True) print(html('The collection of solutions to grad(f)=0')) show((sol[x],sol[y]) for sol in solns)

The collection of solutions to grad(f)=0

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|<generator|\phantom{\verb!x!}\verb|object|\phantom{\verb!x!}\verb|<genexpr>|\phantom{\verb!x!}\verb|at|\phantom{\verb!x!}\verb|0x7fa86e6f9eb0>|$

The collection of solutions to grad(f)=0

# To display the results nicely to determine extreme nature... numerical = 'true' for sol in solns: if numerical: html('At critical point $%s$,'%str(latex((sol[x],sol[y])))) d = (fxx(x=sol[x],y=sol[y])*fyy(x=sol[x],y=sol[y])-fxy(x=sol[x],y=sol[y])^2).n() html(' d = $%s$ '%str(latex(d))+'with $f_{xx} = %s.'%str(latex((fxx(x=sol[x],y=sol[y])).n()))) else: html('At critical point $%s$,'%str(latex((sol[x],sol[y])))) d = (fxx(x=sol[x],y=sol[y])*fyy(x=sol[x],y=sol[y])-fxy(x=sol[x],y=sol[y])^2) html(' d = $%s$ '%str(latex(d))+'with $f_{xx} = %s.'%str(latex((fxx(x=sol[x],y=sol[y])))))

At critical point $\left(0, 0\right)$ ,

Traceback (click to the left of this block for traceback)
...
ValueError: power::eval(): division by zero

At critical point ,

Traceback (most recent call last):            d = (fxx(x=sol[x],y=sol[y])*fyy(x=sol[x],y=sol[y])-fxy(x=sol[x],y=sol[y])^2).n()
  File "", line 1, in <module>
    
  File "/tmp/tmpSltbGG/___code___.py", line 5, in <module>
    exec compile(u"for sol in solns:\n    if numerical:\n        html('At critical point $%s$,'%str(latex((sol[x],sol[y]))))\n        d = (fxx(x=sol[x],y=sol[y])*fyy(x=sol[x],y=sol[y])-fxy(x=sol[x],y=sol[y])**_sage_const_2 ).n()\n        html('     d = $%s$ '%str(latex(d))+'with $f_{xx} = %s.'%str(latex((fxx(x=sol[x],y=sol[y])).n())))\n    else:\n        html('At critical point $%s$,'%str(latex((sol[x],sol[y]))))\n        d = (fxx(x=sol[x],y=sol[y])*fyy(x=sol[x],y=sol[y])-fxy(x=sol[x],y=sol[y])**_sage_const_2 )\n        html('     d = $%s$ '%str(latex(d))+'with $f_{xx} = %s.'%str(latex((fxx(x=sol[x],y=sol[y])))))" + '\n', '', 'single')
  File "", line 4, in <module>
    
  File "sage/symbolic/expression.pyx", line 5437, in sage.symbolic.expression.Expression.__call__ (build/cythonized/sage/symbolic/expression.cpp:30903)
  File "sage/symbolic/ring.pyx", line 981, in sage.symbolic.ring.SymbolicRing._call_element_ (build/cythonized/sage/symbolic/ring.cpp:11110)
  File "sage/symbolic/expression.pyx", line 5291, in sage.symbolic.expression.Expression.substitute (build/cythonized/sage/symbolic/expression.cpp:30025)
ValueError: power::eval(): division by zero

# Determine the ranges of the extrema components to get bounds for a nice graph minx = min(sol[x] for sol in solns) maxx = max(sol[x] for sol in solns) miny = min(sol[y] for sol in solns) maxy = max(sol[y] for sol in solns) print 'range of x values: ',(minx,maxx) print 'range of y values: ',(miny,maxy)

range of x values:   (-2*sqrt(3), 2*sqrt(3))
range of y values:   (-2, 2)

range of x values:   (-2*sqrt(3), 2*sqrt(3))
range of y values:   (-2, 2)

# Use this cell if the above cell returns a range of values that is not definite # minx = # maxx = # miny = -pi # maxy = pi # Solns = needs to be a finite string of values...how to get rid of the possible z?? terms?

print 'A nice plot of the given surface with critical points.' Surf = plot3d(f,(x,minx-1,maxx+1),(y,miny-1,maxy+1),color='yellow') delx = (maxx-minx)/6 dely = (maxy-miny)/6 Points = [] G = Graphics() for sol in solns: z = f(x=sol[x],y=sol[y]) Points.append((sol[x],sol[y],z)) G += point3d((sol[x],sol[y],z),size=20) G += text3d('(%s,'%str(sol[x])+'%s,'%str(sol[y])+'%s)'%str(z),(sol[x]+delx,sol[y]+dely,z+1),horizontal_alignment='left', vertical_alignment='bottom') # Dots = point3d(Points,size=10) (Surf+G).show()

A nice plot of the given surface with critical points.

A nice plot of the given surface with critical points.

# The whole thing in one cell var('x,y') @interact def _(f = input_box(2*x^3-3*x^2*y+3*x^2-y^3+3*y-4,label='$f(x,y) = $')): fx = diff(f,x) fy = diff(f,y) grad = [fx,fy] fxx = diff(fx,x) fyy = diff(fy,y) fxy = diff(fx,y) solns = solve(grad,x,y,solution_dict=True) for sol in solns: html('At critical point $%s$,'%str(latex((sol[x],sol[y])))) d = (fxx(x=sol[x],y=sol[y])*fyy(x=sol[x],y=sol[y])-fxy(x=sol[x],y=sol[y])^2) html(' d = $%s$ '%str(latex(d))+'with $f_{xx} = %s.'%str(latex(fxx(x=sol[x],y=sol[y])))) solns = solve(grad,x,y,solution_dict=True) minx = min(sol[x] for sol in solns) maxx = max(sol[x] for sol in solns) miny = min(sol[y] for sol in solns) maxy = max(sol[y] for sol in solns) print html('A nice plot of $f(x,y) = %s$ with critical points.'%str(latex(f))) Surf = plot3d(f,(x,minx-1,maxx+1),(y,miny-1,maxy+1),color='yellow') delx = (maxx-minx)/6 dely = (maxy-miny)/6 Points = [] G = Graphics() for sol in solns: z = f(x=sol[x],y=sol[y]) Points.append((sol[x],sol[y],z)) G += point3d((sol[x],sol[y],z),size=20) G += text3d('(%s,'%str(sol[x])+'%s,'%str(sol[y])+'%s)'%str(z),(sol[x]+delx,sol[y]+dely,z+1),horizontal_alignment='left', vertical_alignment='bottom') # Dots = point3d(Points,size=10) (Surf+G).show()

$f(x,y) = $

Click to the left again to hide and once more to show the dynamic interactive window

var('x,y,z') implicit_plot3d((x^2 + y^2 + z^2), (x, -2, 2), (y, -2, 2), (z, -2, 2), plot_points=60, contour=[1,3,5], region=lambda x,y,z: x<=0.2 or y>=0.2 or z<=0.2).show()

222 - Topic 09 - Multivariate Extrema

1713 days ago by Professor222

Click to the left again to hide and once more to show the dynamic interactive window