# 1.10 Lesson 10

```# 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, \

** 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
>

```