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