# Lesson 11. Critical Points, continued # ============================ # > restart: with(plots): with(linalg): # We were classifying the critical points of the following function. > f:= (x,y) -> sin(y-x^2-1) + cos(2*y^2-x); > eq1:= D[1](f)(x,y)=0; eq2:= D[2](f)(x,y)=0; # From these we derived the following equations: > eq3:= factor(eq1 + 2*x*eq2); > eq4:= factor(4*y*eq1 + eq2); # We had found and classified the critical points where sin(2 y^2-x) = 0 # and cos(y-x^2-1)=0. # These are on the intersections of the parabolas y - x^2 - 1 = (n+1/2) # Pi and 2 y^2 - x = m Pi for integers m and n. # Now we will look at critical points where -1 + 8 x y = 0. Among # them is one point we found with "fsolve": > s2:= fsolve({eq1,eq2},{x,y},x=-1.5 .. -0.5, y=-0.5 .. 0.5); # Here is a picture of the parabolas and the hyperbola -1 + 8 x y = 0 # together with the wiggly curves defined by eq1 and eq2. > display({ implicitplot({eq1, eq2}, x=-3 .. 3, y = -3 .. 3, \ colour=red),\ plot({seq(x^2 + 1 + (nn-1/2)*Pi, nn= -3 .. 0)},\ x=-3 .. 3, colour=green), > plot({seq([2*y^2 + mm*Pi, y, y=-3 .. 3], mm = -6 .. 0)}, colour=green),\ plot(1/(8*x), x=-3 .. 3, colour=blue) },\ view=[-3 .. 3, -3 .. 3], axes = BOXED); # Note that if -1 + 8 x y = 0, eq3 and eq4 are both satisfied but eq1 and # eq2 may not be. > subs(y=1/(8*x), [eq1, eq2]); > h:= unapply(lhs("[1]), x); # So the critical points with -1 + 8 x y = 0 correspond to the zeros of # h(x). There's not much hope of finding these exactly. > plot(h(x), x=-4 .. 4); # It looks a bit confused for x near 0. Actually it looks worse if you look # closer. > plot(h(x), x = - 0.2 .. 0.2); # Of course x near 0 means |y| is large. > g:= unapply(subs(x=1/(8*y), h(x)), y); > plot(g(y), y=1 .. 8); # We'll look at g(y) for large y. When y is large, the 1/y # factor makes the cos term small, so the sin term # dominates. Inside that term, the 2 y^2 dominates # the 1/(8 y). So g(y) is approximately - sin(- 2 y^2) # = sin(2 y^2). In fact, # |g(y) - sin(2 y^2)| <= 1/(4 |y|) + 1/(8 |y|) = 3/(8 |y|) # So we should look for a zero of g(y) near # 2 y^2 = n Pi, i.e. y = (+|-) sqrt(n Pi/2), for any positive # integer n. # It can be shown that the error in this approximation should # be of order 1/n as n -> infinity. Let's see, for n from 1 to 20, # whether g(sqrt(n Pi/2) - 1/n) and g(sqrt(n Pi/2) + 1/n) # have different signs. If so, there will be a zero of g somewhere # in the interval between sqrt(n Pi/2) - 1/n and ... + 1/n. > plot({seq([[nn, g(sqrt(nn*Pi/2)+1/nn)], \ [nn, g(sqrt(nn*Pi/2)-1/nn)]], nn=1..20)}); # Apparently it's true for n > 1. # Now what about the classification of the critical point? > H:= hessian(f(x,y), [x,y]); # When y is large and x = 1/(8 y), this is approximately: > Happ:= matrix([[ -2*cos(-y+1) - cos(2*y^2),\ 4*cos(2*y^2)*y],\ [ 4*cos(2*y^2)*y,\ -16*cos(2*y^2)*y^2]]); > det(Happ); > Happ[1,1]; # Here's what we can expect: # a local max near y=sqrt(n Pi/2) when # cos(sqrt(n Pi/2)-1) > 0 and cos(n Pi) = +1 # a local min when cos(sqrt(n Pi/2)-1) < 0 and cos(n Pi) = -1 # a saddle point when cos(sqrt(n Pi/2)-1) and cos(n Pi) have opposite # signs. > predict:= n -> if # not type(n,numeric) then 'predict(n)' elif\ evalf(cos(sqrt(n*Pi/2)-1)*(-1)^n) < 0\ then 'saddle'\ elif (-1)^n = 1 then 'max'\ else 'min' fi; > predict(5); # Usually cos(sqrt(n Pi/2)-1) and cos(sqrt((n+1) Pi/2)-1) have the # same sign, so local max or min's and saddle points alternate. # cos(sqrt(n Pi/2)-1) and cos(sqrt((n+1) Pi/2)-1) have different signs if # sqrt(n Pi/2) -1 < (m + 1/2) Pi < sqrt((n+1) Pi/2) -1 for some integer # m. This means n is the greatest integer # <= 2 Pi (m+1/2 + 1/Pi)^2. > ys:= [seq(fsolve(g(y), y, \ y= sqrt(nn*Pi/2)- 1/nn .. sqrt(nn*Pi/2) + 1/nn), \ nn = 2 .. 22)]; > xs:= map(y -> 1/(8*y), "); > evalf(seq([nn, \ eigenvals(subs(x=xs[nn-1], y=ys[nn-1], eval(H))), \ predict(nn) ], nn=2 .. 22)); > seq(floor(2 *Pi *(mm+1/2 + 1/Pi)^2), mm=0..5); >