# Critical Points of Functions of Two and Three Variables

copyright © 2009 by Jonathan Rosenberg based on an earlier M-book, copyright © 2000 by Paul Green and Jonathan Rosenberg

## Contents

## Finding and Classifying Critical Points

We recall that a critical point of a function of several variables is a point at which the gradient of the function is either
the zero vector **0** or is undefined. Although every point at which a function takes a local extreme value is a critical point, the converse is
not true, just as in the single variable case. In this lesson we will be interested in identifying critical points of a function
and classifying them. As in the single variable case, since the first partial derivatives vanish at every critical point,
the classification depends on the values of the second partial derivatives. Let us illustrate what is involved with a specific
example, namely the function defined below, which we have been considering throughout the past several notebooks.

syms x y z f=((x^2-1)+(y^2-4)+(x^2-1)*(y^2-4))/(x^2+y^2+1)^2

f = ((x^2 - 1)*(y^2 - 4) + x^2 + y^2 - 5)/(x^2 + y^2 + 1)^2

To find the critical points of `f` we must set both partial derivatives of `f` to 0 and solve for `x` and `y`. We begin by computing the first partial derivatives of `f`.

fx=diff(f,x) fy=diff(f,y)

fx = (2*x + 2*x*(y^2 - 4))/(x^2 + y^2 + 1)^2 - (4*x*((x^2 - 1)*(y^2 - 4) + x^2 + y^2 - 5))/(x^2 + y^2 + 1)^3 fy = (2*y + 2*y*(x^2 - 1))/(x^2 + y^2 + 1)^2 - (4*y*((x^2 - 1)*(y^2 - 4) + x^2 + y^2 - 5))/(x^2 + y^2 + 1)^3

To find critical points of `f`, we must set the partial derivatives equal to 0 and solve for `x` and `y`. Since the equations in this case are algebraic, we can use ` solve`. (For more complicated functions built in part out of transcendental functions like

`exp`,

`log`, trig functions, etc., it may be be necessary to solve for the critical points numerically. More about this later.)

[xcr,ycr]=solve(fx,fy)

xcr = 0 3^(1/2)/3 -3^(1/2)/3 (2^(1/2)*i)/2 -(2^(1/2)*i)/2 (2^(1/2)*i)/2 -(2^(1/2)*i)/2 ycr = 0 0 0 10^(1/2)/2 10^(1/2)/2 -10^(1/2)/2 -10^(1/2)/2

Incidentally, it is worth pointing out a trick here. By default, the outputs `xcr` and `ycr` of the ` solve` command are column vectors. As we have just typed the

`command, the vectors appear one after another. But a better alternative is:`

**solve**[xcr,ycr]=solve(fx,fy); [xcr,ycr]

ans = [ 0, 0] [ 3^(1/2)/3, 0] [ -3^(1/2)/3, 0] [ (2^(1/2)*i)/2, 10^(1/2)/2] [ -(2^(1/2)*i)/2, 10^(1/2)/2] [ (2^(1/2)*i)/2, -10^(1/2)/2] [ -(2^(1/2)*i)/2, -10^(1/2)/2]

Now `xcr` and `ycr` are stacked together in a matrix, so we can see each x-value together with the corresponding y-value. MATLAB has found seven
critical points, of which the last four are not real. We now compute the second derivatives of f.

fxx=diff(fx,x) fxy=diff(fx,y) fyy=diff(fy,y)

fxx = (2*y^2 - 6)/(x^2 + y^2 + 1)^2 - (4*((x^2 - 1)*(y^2 - 4) + x^2 + y^2 - 5))/(x^2 + y^2 + 1)^3 + (24*x^2*((x^2 - 1)*(y^2 - 4) + x^2 + y^2 - 5))/(x^2 + y^2 + 1)^4 - (8*x*(2*x + 2*x*(y^2 - 4)))/(x^2 + y^2 + 1)^3 fxy = (4*x*y)/(x^2 + y^2 + 1)^2 - (4*y*(2*x + 2*x*(y^2 - 4)))/(x^2 + y^2 + 1)^3 - (4*x*(2*y + 2*y*(x^2 - 1)))/(x^2 + y^2 + 1)^3 + (24*x*y*((x^2 - 1)*(y^2 - 4) + x^2 + y^2 - 5))/(x^2 + y^2 + 1)^4 fyy = (2*x^2)/(x^2 + y^2 + 1)^2 - (4*((x^2 - 1)*(y^2 - 4) + x^2 + y^2 - 5))/(x^2 + y^2 + 1)^3 + (24*y^2*((x^2 - 1)*(y^2 - 4) + x^2 + y^2 - 5))/(x^2 + y^2 + 1)^4 - (8*y*(2*y + 2*y*(x^2 - 1)))/(x^2 + y^2 + 1)^3

Of course there is no need to compute `fyx`, since for reasonable functions it always coincides with `fxy`. We compute the ` Hessian determinant` (called the discriminant in your text) of the second partial derivatives,

`fxx*fyy-fxy^2`, which is essential for the classification, and evaluate it at the critical points.

hessdetf=fxx*fyy-fxy^2

hessdetf = ((2*y^2 - 6)/(x^2 + y^2 + 1)^2 - (4*((x^2 - 1)*(y^2 - 4) + x^2 + y^2 - 5))/(x^2 + y^2 + 1)^3 + (24*x^2*((x^2 - 1)*(y^2 - 4) + x^2 + y^2 - 5))/(x^2 + y^2 + 1)^4 - (8*x*(2*x + 2*x*(y^2 - 4)))/(x^2 + y^2 + 1)^3)*((2*x^2)/(x^2 + y^2 + 1)^2 - (4*((x^2 - 1)*(y^2 - 4) + x^2 + y^2 - 5))/(x^2 + y^2 + 1)^3 + (24*y^2*((x^2 - 1)*(y^2 - 4) + x^2 + y^2 - 5))/(x^2 + y^2 + 1)^4 - (8*y*(2*y + 2*y*(x^2 - 1)))/(x^2 + y^2 + 1)^3) - (4*x*(2*y + 2*y*(x^2 - 1))/(x^2 + y^2 + 1)^3 + 4*y*(2*x + 2*x*(y^2 - 4))/(x^2 + y^2 + 1)^3 - 4*x*y/(x^2 + y^2 + 1)^2 - 24*x*y/(x^2 + y^2 + 1)^4*((x^2 - 1)*(y^2 - 4) + x^2 + y^2 - 5))^2

Another way to do the calculations, which is a little more elegant, is to use the ` jacobian` command. Recall that this computes the gradient. But by using it twice, we can get the matrix of 2nd derivatives and then
take the determinant.

gradf = jacobian(f, [x, y]) hessmatf = jacobian(gradf, [x, y]) hessdetf1 = simplify(det(hessmatf)) simplify(hessdetf - hessdetf1)

gradf = [ (2*x + 2*x*(y^2 - 4))/(x^2 + y^2 + 1)^2 - (4*x*((x^2 - 1)*(y^2 - 4) + x^2 + y^2 - 5))/(x^2 + y^2 + 1)^3, (2*y + 2*y*(x^2 - 1))/(x^2 + y^2 + 1)^2 - (4*y*((x^2 - 1)*(y^2 - 4) + x^2 + y^2 - 5))/(x^2 + y^2 + 1)^3] hessmatf = [ (2*y^2 - 6)/(x^2 + y^2 + 1)^2 - (4*((x^2 - 1)*(y^2 - 4) + x^2 + y^2 - 5))/(x^2 + y^2 + 1)^3 + (24*x^2*((x^2 - 1)*(y^2 - 4) + x^2 + y^2 - 5))/(x^2 + y^2 + 1)^4 - (8*x*(2*x + 2*x*(y^2 - 4)))/(x^2 + y^2 + 1)^3, (4*x*y)/(x^2 + y^2 + 1)^2 - (4*y*(2*x + 2*x*(y^2 - 4)))/(x^2 + y^2 + 1)^3 - (4*x*(2*y + 2*y*(x^2 - 1)))/(x^2 + y^2 + 1)^3 + (24*x*y*((x^2 - 1)*(y^2 - 4) + x^2 + y^2 - 5))/(x^2 + y^2 + 1)^4] [ (4*x*y)/(x^2 + y^2 + 1)^2 - (4*y*(2*x + 2*x*(y^2 - 4)))/(x^2 + y^2 + 1)^3 - (4*x*(2*y + 2*y*(x^2 - 1)))/(x^2 + y^2 + 1)^3 + (24*x*y*((x^2 - 1)*(y^2 - 4) + x^2 + y^2 - 5))/(x^2 + y^2 + 1)^4, (2*x^2)/(x^2 + y^2 + 1)^2 - (4*((x^2 - 1)*(y^2 - 4) + x^2 + y^2 - 5))/(x^2 + y^2 + 1)^3 + (24*y^2*((x^2 - 1)*(y^2 - 4) + x^2 + y^2 - 5))/(x^2 + y^2 + 1)^4 - (8*y*(2*y + 2*y*(x^2 - 1)))/(x^2 + y^2 + 1)^3] hessdetf1 = (x^6*(4*y^4 + 104*y^2 - 196) + x^4*(4*y^6 + 284*y^4 - 652*y^2 + 316) - x^2*(4*y^8 - 32*y^6 + 64*y^4 + 512*y^2 - 84) - x^8*(4*y^2 + 36) + 24*y^2 + 88*y^4 - 40*y^6 - 8)/(x^2 + y^2 + 1)^7 ans = 0

The fact that `hessdetf - hessdetf1` simplified to 0 showed that both methods give the same thing.

The next step is to make a table of critical points, along with values of the Hessian determinant and the second partial derivative
with respect to `x`. Since we are not interested in then non-real critical points, we can truncate the vectors `xcr` and `ycr` after the first three entries.

xcr = xcr(1:3); ycr = ycr(1:3); for k = 1:3 [xcr(k), ycr(k), subs(hessdetf, [x,y], [xcr(k), ycr(k)]), ... subs(fxx, [x,y], [xcr(k), ycr(k)])] end

ans = [ 0, 0, -8, -2] ans = [ 3^(1/2)/3, 0, 405/64, 27/16] ans = [ -3^(1/2)/3, 0, 405/64, 27/16]

Notice that all three of the real critical points are on the x-axis, and the first one is at the origin. Since the Hessian
determinant is negative at the origin, we conclude that the critical point at the origin is a saddle point. The other two
are local extrema and, since `fxx` is positive at both of them, they are local minima.

Let now reconstruct the contour and gradient plot we obtained for `f` in the last lesson.

[xx, yy] = meshgrid(-3:.1:3,-3:.1:3); ffun = @(x,y) eval(vectorize(f)); fxfun = @(x,y) eval(vectorize(gradf(1))); fyfun = @(x,y) eval(vectorize(gradf(2))); contour(xx, yy, ffun(xx,yy), 30) hold on [xx, yy] = meshgrid(-3:.25:3,-3:.25:3); quiver(xx, yy, fxfun(xx,yy), fyfun(xx,yy), 0.6) axis equal tight, hold off

From the plot it is evident that there are two symmetrically placed local minima on the x-axis with a saddle point between
them. We recognize the minima from the fact that they are surrounded by closed contours with the gradient arrows pointing
**outward**. The pinching in of contours near the origin indicates a saddle point. The contour through the origin crosses itself, forming
a transition between contours enclosing the two minima separately and contours enclosing them together.

## Problem 1:

Find and classify all critical points of the function

wich you studied in the previous lesson, and compare your results with the contour and gradient plot of `g` that you have already made.

## Eigenvalues of the Hessian

We now take another look at the classification of critical points, with a view to extending it to functions of three or more
variables. Let us define the function f2 and compute its first and second derivatives using MATLAB's symbolic ` jacobian` operator.

f2=x^3-3*x^2+5*x*y-7*y^2 gradf2=jacobian(f2,[x,y]) hessmatf2=jacobian(gradf2,[x,y])

f2 = x^3 - 3*x^2 + 5*x*y - 7*y^2 gradf2 = [ 3*x^2 - 6*x + 5*y, 5*x - 14*y] hessmatf2 = [ 6*x - 6, 5] [ 5, -14]

The last output above is a symmetric matrix, known as the Hessian matrix, whose entries are the second partial derivatives
of `f`. In order to find the critical points, we must again set the first partial derivatives equal to 0.

[xcr2,ycr2]=solve(gradf2(1),gradf2(2)); [xcr2,ycr2]

ans = [ 0, 0] [ 59/42, 295/588]

There are two critical points. Let us now evaluate the Hessian matrix at the critical points.

H1=subs(hessmatf2,[x,y], [xcr2(1),ycr2(1)]) H2=subs(hessmatf2,[x,y], [xcr2(2),ycr2(2)])

H1 = [ -6, 5] [ 5, -14] H2 = [ 17/7, 5] [ 5, -14]

At this point, however, instead of computing the Hessian determinant, which is the determinant of the Hessian matrix, we compute
the **eigenvalues** of the matrix. The determinant of a matrix is the *product* of its eigenvalues, but the eigenvalues carry a lot more information
than the determinant alone.

eig(H1) eig(H2)

ans = - 41^(1/2) - 10 41^(1/2) - 10 ans = - (25*29^(1/2))/14 - 81/14 (25*29^(1/2))/14 - 81/14

We can obtain numerical versions with the ` double` command:

double(eig(H1)) double(eig(H2))

ans = -16.4031 -3.5969 ans = -15.4021 3.8307

The significance of the eigenvalues of the Hessian matrix is that if all of them are positive at a critical point, the function
has a local minimum there; if all are negative, the function has a local maximum; if they have mixed signs, the function has
a saddle point; and if at least one of them is 0, the critical point is degenerate. These statements **all generalize to functions of more than two variables**, unlike the rule about the Hessian determinant. For functions of two variables, the determinant, being the product of the
two eigenvalues, takes a negative value if and only if they have opposite signs. In the present case, we see that the critical
point at the origin is a local maximum of `f2`, and the second critical point is a saddle point.

The point of this reformulation is that the Hessian matrix and its eigenvalues are just as easily computed and interpreted in the three-dimensional case. In this case, there will be three eigenvalues, all of which are positive at a local minimum and negative at a local maximum. A critical point of a function of three variables is degenerate if at least one of the eigenvalues of the Hessian matrix is 0, and has a saddle point in the remaining case, when at least one eigenvalue is positive, at least one is negative, and none is 0.

## Problem 2:

Find and classify the critical points of the function

You may have trouble substituting the symbolic values of the coordinates of the critical points into the Hessian matrix; if
so, convert them to numerical values, using ` double`.

## A Graphical/Numerical Method

For some functions (especially if they involve transcendental functions such as `exp`, `log`, `sin`, `cos`, etc.), the formulas for the partial derivatives may be too complicated to use solve to find the critical points. In this
case, a graphical or numerical method may be necessary. Here is an example:

k=x*sin(y)+y^2*cos(x) kx=diff(k,x) ky=diff(k,y)

k = y^2*cos(x) + x*sin(y) kx = sin(y) - y^2*sin(x) ky = x*cos(y) + 2*y*cos(x)

In this case solving analytically for the solutions of `kx = 0`, `ky = 0` is going to be hopeless. Also, we can't solve numerically without knowing how many solutions there are and roughly where
they are located. The way out is therefore to draw plots of the solutions to `kx = 0` and to `ky = 0`. Then we see where these curves intersect and can use the graphical information to get starting values for the use of a
numerical solver. In our particular example, we can try:

[x1,y1]=meshgrid(-8:.1:8,-8:.1:8); kxfun=@(x,y) eval(vectorize(kx)); kyfun=@(x,y) eval(vectorize(ky)); contour(x1,y1,kxfun(x1,y1),[0,0],'r'), hold on contour(x1,y1,kyfun(x1,y1),[0,0],'b'), hold off

Here the plot of `kx = 0` is shown in red and the plot of `ky = 0` is shown in blue. Thus, for example, we have critical points near (6.2, 2.4) and near (6.2, 3.2). We can locate them more
accurately using the M-file `newton2d.m`, which uses a two-dimensional version of Newton's method to solve more accurately:

```
type newton2d.m
```

function [xsol, ysol] = newton2d(varargin) %NEWTON2D 2-dimensional Newton's method % The command % [xsol, ysol] = newton2d(eq1, eq2, xstart, ystart) % searches for a root of the simultaneous (nonlinear) % equations eq1=0 and eq2=0, starting from the initial % guesses x=xstart, y=ystart. Here eq1 and eq2 should be % entered as strings in the variables 'x' and 'y' or % as symbolic expressions in x and y. Alternatively, % [xsol, ysol] = newton2d(eq1, eq2, xvar, yvar, xstart, ystart) % allows one to specify the names of the independent variables. % The optional argument 'maxiterations' at the end specifies the % maximum number of iterations to try. The default is 20. % % Example: [xsol, ysol] = newton2d('x^2 + y^2 - 4', 'x - y', 1, 1) % gives the output % xsol = % % 1.4142 % % % ysol = % % 1.4142 % % % written by Jonathan Rosenberg, 8/10/99 % revised 8/22/05 if nargin<4 error('not enough input arguments -- need at least the two equations and starting values for x and y'); end if nargin>7, error('too many input arguments'); end eq1=varargin{1}; eq2=varargin{2}; if ischar(eq1) & ~ischar(eq2) error('both equations must be of same form: strings or symbolic') end if ~ischar(eq1) & ischar(eq2) error('both equations must be of same form: strings or symbolic') end % Default value of maxiterations is 20. maxiterations=20; % Default values of xvar and yvar are 'x', 'y'. if nargin<6 xvar = 'x'; yvar = 'y'; xstart = varargin{3}; ystart = varargin{4}; end if nargin>5 xvar = varargin{3}; yvar = varargin{4}; xstart = varargin{5}; ystart = varargin{6}; end if nargin==5, maxiterations=varargin{5}; end if nargin==7, maxiterations=varargin{7}; end % Start by computing partial derivatives. f11 = diff(sym(eq1), xvar); f12 = diff(sym(eq1), yvar); f21 = diff(sym(eq2), xvar); f22 = diff(sym(eq2), yvar); % Vector functions for evaluation F1 = inline(vectorize(eq1),xvar,yvar); F2 = inline(vectorize(eq2),xvar,yvar); A11 = inline(vectorize(f11),xvar,yvar); A12 = inline(vectorize(f21),xvar,yvar); A21 = inline(vectorize(f12),xvar,yvar); A22 = inline(vectorize(f22),xvar,yvar); X = [xstart, ystart]; for counter=1:maxiterations A = [A11(X(1), X(2)), A12(X(1), X(2)); A21(X(1), X(2)), A22(X(1), X(2))]; F = [F1(X(1), X(2)), F2(X(1), X(2))]; if max(abs(F)) < eps, xsol=X(1); % process has converged ysol=X(2); return; end if condest(A) > 10^10, error('matrix of partial derivatives singular or almost singular -- try again with different starting values'); end X = X - F/A; end xsol=X(1); ysol=X(2);

[xsol,ysol]=newton2d(kx, ky, 6.2, 2.4); [xsol,ysol] [xsol,ysol]=newton2d(kx, ky, 6.2, 3.2); [xsol,ysol]

ans = 6.3954 2.4237 ans = 6.2832 3.1416

Thus we have critical points near (6.3954, 2.4237) and (6.2832, 3.1416). The second of these looks suspiciously like (2*pi, pi), and in fact we have

subs(kx, [x,y], sym('[2*pi,pi]')) subs(ky, [x,y], sym('[2*pi,pi]'))

ans = 0 ans = 0

so this hunch is in fact correct. To classify these critical points, we can compute

hessdetk = simplify(diff(kx, x)*diff(ky, y) - diff(kx, y)^2) double(subs(hessdetk, [x,y], [6.3954, 2.4237])) double(subs(diff(kx, x), [x,y], [6.3954, 2.4237])) double(subs(hessdetk, [x,y], [2*pi, pi]))

hessdetk = - (cos(y) - 2*y*sin(x))^2 - y^2*cos(x)*(2*cos(x) - x*sin(y)) ans = 11.2762 ans = -5.8374 ans = -20.7392

so the first point, (6.3954, 2.4237), has positive Hessian determinant, negative 2nd derivative with respect to `x`, and is thus a local maximum. Similarly the point (2*pi, pi) is a saddle point. One could also check these calculations
with a contour plot:

ezcontour(k, [6, 6.5, 2, 4])

The picture indeed shows contours looping around the local maximum at (6.3954, 2.4237) and the characteristic shape of a saddle near (6.2832, 3.1416).

## Additional Problems:

- 1. Find and classify all critical points of the function

Relate your results to a simultaneous contour and gradient plot of the function.

- 2. Find and classify all critical points of the function

MATLAB will report many critical points, but only a few of them are real.

- 3. Find and classify all critical points of the function
`h(x, y) = y^2*exp(x^2) - x - 3*y`. You will need the graphical/numerical method to find the critical points.

- 4. Find and classify all critical points of the function
`h(x, y) = y^2*exp(x) - x - 3*y`. Use both the analytical and the graphical/numerical methods to find the critical points, and compare the results.