# Lesson 10. Critical Points
# ==================
#
> restart: with(plots):
# We wanted to examine the critical points of the following
# function.
> f:=(x,y) -> sin(y-x^2-1) + cos(2*y^2-x);
2 2
f := (x,y) -> sin(y - x - 1) + cos(2 y - x)
# We made the following contour plot.
> contourplot(f(x,y), x=-3 .. 3,y=-3 .. 3, \
grid=[50,50], axes=BOXED, shading=ZHUE);
** Maple V Graphics **
# The critical points are the points where the partial derivatives of
# f with respect to x and y are 0.
** Maple V Graphics **
> eq1:= D[1](f)(x,y)=0; eq2:= D[2](f)(x,y)=0;
2 2
eq1 := - 2 cos(- y + x + 1) x - sin(- 2 y + x) = 0
2 2
eq2 := cos(- y + x + 1) + 4 sin(- 2 y + x) y = 0
# The critical points are the intersections of these curves.
> display({ \
implicitplot(eq1, x=-3 .. 3, y=-3 .. 3, colour=red), \
implicitplot(eq2, x=-3 .. 3, y=-3 .. 3, colour=blue)});
# We found some solutions using "fsolve".
> s1:= fsolve({eq1,eq2},{x,y});
s1 := {y = -1.646743751, x = 2.281937309}
> s2:= fsolve({eq1,eq2},{x,y},x=-1.5 .. -0.5,
> y=-0.5 .. 0.5);
s2 := {x = -.9452256199, y = -.1322435590}
#
# To classify a critical point, we use the second derivative test.
# The several variable version of this uses the Hessian matrix,
# which is the matrix of second partial derivatives of the function.
#
# Almost all of the linear algebra commands are in the
# "linalg" package.
> with(linalg):
Warning: new definition for norm
Warning: new definition for trace
# The "matrix" command produces a matrix. You
# enter the matrix as a list of lists, each representing
# a row.
> matrix([[a,b],[c,d]]);
[ a b ]
[ ]
[ c d ]
# We could make the Hessian matrix of a function of
# two variables as follows:
> matrix([[ D[1,1](g)(x,y), D[1,2](g)(x,y) ], \
[ D[2,1](g)(x,y), D[2,2](g)(x,y) ]]);
[ D[1, 1](g)(x, y) D[1, 2](g)(x, y) ]
[ ]
[ D[1, 2](g)(x, y) D[2, 2](g)(x, y) ]
# It's easier, however, to use the ready-made "hessian"
# command. The first input is an expression in several variables,
# and the second is a list of the variables.
> hessian(g(x,y), [x,y]);
[ 2 2 ]
[ d d ]
[ ----- g(x, y) ------- g(x, y) ]
[ 2 dy dx ]
[ dx ]
[ ]
[ 2 2 ]
[ d d ]
[ ------- g(x, y) ----- g(x, y) ]
[ dy dx 2 ]
[ dy ]
# Here it is for our f:
> H:= hessian(f(x,y),[x,y]);
H := 2 2 2 2
[4 sin(- y + x + 1) x - 2 cos(- y + x + 1) - cos(- 2 y + x),
2 2
- 2 sin(- y + x + 1) x + 4 cos(- 2 y + x) y]
2 2
[- 2 sin(- y + x + 1) x + 4 cos(- 2 y + x) y,
2 2 2 2
sin(- y + x + 1) - 16 cos(- 2 y + x) y + 4 sin(- 2 y + x)]
# The entries are too long to display nicely as a matrix,
# but this is a 2 by 2 matrix. I assigned it to the variable H.
#
# Matrices and vectors have some rather odd behaviour in Maple.
# Here's what happens if we take another look at H.
> H;
H
# It just returns the name H. If we want to see the entries of the
# matrix, we must use "eval".
> eval(H);
2 2 2 2
[4 sin(- y + x + 1) x - 2 cos(- y + x + 1) - cos(- 2 y + x),
2 2
- 2 sin(- y + x + 1) x + 4 cos(- 2 y + x) y]
2 2
[- 2 sin(- y + x + 1) x + 4 cos(- 2 y + x) y,
2 2 2 2
sin(- y + x + 1) - 16 cos(- 2 y + x) y + 4 sin(- 2 y + x)]
# Now to classify a critical point [x0, y0], we want to substitute x0
# for x and y0 for y in H. More precisely, in eval(H) (because we
# want the entries of H and not just its name).
> subs(s1, eval(H));
[20.82895153 sin(7.853981633) - 2 cos(7.853981633) - cos(-3.141592653),
- 4.563874618 sin(7.853981633) - 6.586975004 cos(-3.141592653)]
[ - 4.563874618 sin(7.853981633) - 6.586975004 cos(-3.141592653),
sin(7.853981633) - 43.38823970 cos(-3.141592653) + 4 sin(-3.141592653)]
# Why doesn't this calculate the sin and cos values
# numerically? Because "subs" is a bit peculiar --- it
# does the substitution, but very little simplification
# after. We need "eval" again. But this time (I'm
# not sure why) we actually need "map" to perform
# "eval" on each of the entries of the matrix.
# (While we're here, notice that some of the numbers
# here look familiar!)
> H1:=map(eval,");
[ 21.82895153 2.023100386 ]
H1 := [ ]
[ 2.023100386 44.38823970 ]
> subs(s2, eval(H));
[3.573805890 sin(2.025695032) - 2 cos(2.025695032) - cos(-.9802023377),
1.890451240 sin(2.025695032) - .5289742360 cos(-.9802023377)]
[1.890451240 sin(2.025695032) - .5289742360 cos(-.9802023377),
sin(2.025695032) - .2798137424 cos(-.9802023377) + 4 sin(-.9802023377)]
> H2:=map(eval,");
[ 3.532257777 1.403641139 ]
H2 := [ ]
[ 1.403641139 -2.579950242 ]
# In second-year Calculus you learn the classification of
# a critical point for a function of two variables based on
# the determinant of the Hessian matrix and its upper left
# entry.
# In Linear Algebra (Math 307, if not 221) you learn that (for any
# number of variables) the classification of a critical point is done
# by checking whether the Hessian matrix at the critical point is
# positive definite (a local minimum), negative definite (a local
# maximum), or indefinite (a saddle point). This in turn depends
# on the eigenvalues of the matrix. Maple can calculate the
# eigenvalues with the "eigenvals" command.
> eigenvals(H1);
21.64895749, 44.56823374
# These are both positive, so the matrix H1 is positive
# definite and the critical point is a local minimum.
> eigenvals(H2);
-2.886877627, 3.839185165
# One is positive and the other negative, so the matrix
# H1 is indefinite and the critical point is a saddle point.
#
# If we have a critical point where all eigenvalues are negative,
# the matrix is negative definite and the critical point is a local
# maximum.
#
# If any of the eigenvalues is 0, the critical point is "degenerate"
# and the second derivative test is inconclusive.
#
# Now let's examine our equations again. Even though "solve"
# couldn't handle them, a little more thought can get results.
> eq1; eq2;
2 2
- 2 cos(- y + x + 1) x - sin(- 2 y + x) = 0
2 2
cos(- y + x + 1) + 4 sin(- 2 y + x) y = 0
# We could eliminate the cos terms by adding the first equation
# and 2x times the second.
# Maple actually does let you add equations, or multiply them by
# expressions.
> eq1 + 2*x*eq2;
2 2
- 2 cos(- y + x + 1) x - sin(- 2 y + x)
2 2
+ 2 x (cos(- y + x + 1) + 4 sin(- 2 y + x) y) = 0
# Maple hasn't performed the simplification because it
# hasn't distributed the 2 x over the sum. We can tell
# Maple to do so. Actually, let's use "factor" which will
# also take out the common factor sin(-2 y^2 + x) from
# the two terms that will be left.
> eq3:=factor(");
2
eq3 := sin(- 2 y + x) (- 1 + 8 x y) = 0
# Similarly, we can get rid of the sin terms.
> 4*y*eq1+eq2;
2 2 2
4 y (- 2 cos(- y + x + 1) x - sin(- 2 y + x)) + cos(- y + x + 1)
2
+ 4 sin(- 2 y + x) y = 0
> eq4:=factor(");
2
eq4 := - cos(- y + x + 1) (- 1 + 8 x y) = 0
# We could solve eq3 and eq4 by hand. Can Maple do it?
> solve(eq3, x);
2 1
2 y , ---
8 y
# As usual, "solve" doesn't find all the solutions. It finds x=1/(8 y)
# which makes -1 + 8 x y = 0, but only x = 2 y^2 to make sin(-2
# y^2+x) = 0. Of course the more complete answer would be x = 2
# y^2 + m Pi, where m is any integer.
#
# By the way, the number Pi is capitalized in Maple input. "pi" is
# also printed the same way in output, but it's treated as just
# another name. This is a common cause of mistakes.
> solve(eq4, y);
2 1
- 1/2 Pi + x + 1, ---
8 x
# And for equation 4 we should have y = 1/(8 x) or
# y = x^2 + 1 + (n-1/2) Pi, where n is any integer.
#
# Taking equations 3 and 4 together, there are two families of
# solutions. In the first, x = 2 y^2 + m Pi and y = x^2 + 1 + (n-1/2)
# Pi. This makes both the cos and sin terms in the equations 0.
> case1:=\
subs({-2*y^2+x = m*Pi, -y+x^2 + 1 = (n-1/2)*Pi},\
[eq1, eq2]);
case1 :=
[- 2 cos((n - 1/2) Pi) x - sin(m Pi) = 0, cos((n - 1/2) Pi) + 4 sin(m Pi) y = 0]
# Maple doesn't know that sin(m Pi) and cos((n-1/2) Pi) are 0.
# Actually it's hardly fair to expect that, because we haven't told
# Maple that m and n are integers. Anyway, we can tell Maple
# these facts.
> sin(m*Pi):= 0: sin((n-1/2)*Pi):= (-1)^(n+1):
> cos(m*Pi):= (-1)^m: cos((n-1/2)*Pi):= 0:
# What does this do? Each Maple function has a "remember
# table", where certain known values of the function are stored.
# When you evaluate a function, e.g. f(0), Maple first looks in the
# remember table to see if a value for f(0) is stored there, and if
# so it uses that value. You can put a value for f(0) in the
# remember table by assigning a value to f(0).
#
> case1;
[0 = 0, 0 = 0]
# For each m, x = 2 y^2 + m Pi is the equation of a parabola
# centred on the x axis.
# For each n, y = x^2 + 1 + (n-1/2) Pi is the equation of a parabola
# centred on the y axis.
# These intersect in up to 4 points.
> display({ plot(x^2 + 1 - 3/2*Pi, x=-3.2 .. 3.2),
> plot([2*y^2-Pi, y, y= -2 .. 2]) });
** Maple V Graphics **
# Now what is the classification of these points?
> subs({-2*y^2+x = m*Pi, -y+x^2 + 1 = (n-1/2)*Pi},\
eval(H));
2
[4 sin((n - 1/2) Pi) x - 2 cos((n - 1/2) Pi) - cos(m Pi),
- 2 sin((n - 1/2) Pi) x + 4 cos(m Pi) y]
[- 2 sin((n - 1/2) Pi) x + 4 cos(m Pi) y,
2
sin((n - 1/2) Pi) - 16 cos(m Pi) y + 4 sin(m Pi)]
> Hmn:=map(eval,");
[ (n + 1) 2 m (n + 1) m ]
[ 4 (-1) x - (-1) - 2 (-1) x + 4 (-1) y ]
Hmn := [ ]
[ (n + 1) m (n + 1) m 2 ]
[ - 2 (-1) x + 4 (-1) y (-1) - 16 (-1) y ]
# It's better to use the Calculus version of the second derivative
# test here instead of trying to find the eigenvalues.
> det(Hmn);
(n + 1) 2 m 2 m (n + 1) (n + 1) m
- 64 (-1) x (-1) y - (-1) (-1) + 16 (-1) x (-1) y
> factor(");
m (n + 1) 2
- (-1) (-1) (- 1 + 8 x y)
# The determinant is the product of the eigenvalues.
# If one of m and n is even and the other odd (and
# 8 x y is not 1) the determinant is negative so the
# eigenvalues have opposite sign, the Hessian is
# indefinite, and the critical point is a saddle point.
#
# If m and n are both even or both odd (and 8 x y is
# not 1), the determinant is positive so the eigenvalues
# have the same sign, the Hessian is either negative definite or
# positive definite, and the critical point is either a local minimum
# or a local maximum.
# Which it is can be decided by looking at one of the
# diagonal entries. These have the same sign as
# (-1)^(n+1) = -(-1)^m. So for m and n both even
# the Hessian is negative definite and the critical point
# is a local maximum, while for m and n both odd
# the Hessian is positive definite and the critical point
# is a local minimum.
# Are our solutions (which we saved as s1 and s2) in this family?
> subs(s1,
> [-2*y^2+x = m*Pi, -y+x^2 + 1 = (n-1/2)*Pi]);
[-3.141592653 = m Pi, 7.853981633 = (n - 1/2) Pi]
# Clearly m=-1. I'll solve for n.
> solve("[2],n);
2.999999999
# m=-1 and n=3 are both odd, so this should be a local minimum
# (as we already saw it was).
{2.499999999 = n - .5000000000, -.9999999995 = m}
> subs(s2,
> [-2*y^2+x = m*Pi, -y+x^2 + 1 = (n-1/2)*Pi]);
[-.9802023377 = m Pi, 2.025695032 = (n - 1/2) Pi]
# Well, certainly m is not an integer. Neither is n.
# So this must be in the second family of solutions.
> subs(s2, 8*x*y-1);
0
>
