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