# Worksheet#15
#
# Systems of non-polynomial equations.
#
--------------------------------------------------------------------------------
# We will try to find the critical points of a function f(x,y). First we load the function and
# then plot its graph.
--------------------------------------------------------------------------------
> 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)
--------------------------------------------------------------------------------
> with(plots);
[animate, animate3d, conformal, contourplot, cylinderplot, densityplot,
display, display3d, fieldplot, fieldplot3d, gradplot, gradplot3d,
implicitplot, implicitplot3d, loglogplot, logplot, matrixplot, odeplot,
pointplot, polarplot, polygonplot, polygonplot3d, polyhedraplot, replot,
setoptions, setoptions3d, spacecurve, sparsematrixplot, sphereplot,
surfdata, textplot, textplot3d, tubeplot]
--------------------------------------------------------------------------------
> ?contourplot
--------------------------------------------------------------------------------
> contourplot(f(x,y),x=-3..3,y=-3..3,title=`with default settings`);
--------------------------------------------------------------------------------
** Maple V Graphics **
--------------------------------------------------------------------------------
# This is a bit rough so we refine the grid to 100 by 100. The default is 25 by 25.
--------------------------------------------------------------------------------
> contourplot(f(x,y),x=-3..3,y=-3..3,grid=[100,100],title=`finer grid`);
--------------------------------------------------------------------------------
#
** Maple V Graphics **
--------------------------------------------------------------------------------
#
# That's much better. Even though this is a contour plot, that is the 2 dimensional graph of
# the contour lines f(x,y)=constants, Maple thinks this is a 3d plot. Using the projection
# menu I got the following picture. It looks great in colour. Try the various menu options.
#
#
** Maple V Graphics **
--------------------------------------------------------------------------------
# To find the critical points we compute the partial derivatives of f(x,y) with respect to x
# and y and then set them equal to zero. The solutions are the critical points.
--------------------------------------------------------------------------------
> f1:=D[1](f);f2:=D[2](f);
2 2
f1 := (x,y) -> - 2 cos(y - x - 1) x + sin(2 y - x)
2 2
f2 := (x,y) -> cos(- y + x + 1) + 4 sin(- 2 y + x) y
--------------------------------------------------------------------------------
# To find the critical points we must solve the simultaneous equations f1(x,y)=0, f2(x,y)=0.
# We first plot these curves.
--------------------------------------------------------------------------------
> with(plots):
> F1:=implicitplot(f1(x,y),x=-3..3,y=-3..3,grid=[100,100]):
--------------------------------------------------------------------------------
> F2:=implicitplot(f2(x,y),x=-3..3,y=-3..3,grid=[100,100]):
--------------------------------------------------------------------------------
> display([F1,F2],title=`f1(x,y)=0 and f2(x,y)=0`);
--------------------------------------------------------------------------------
#
** Maple V Graphics **
--------------------------------------------------------------------------------
# Every intersectionpoint is a critical point, there are lots of them. Using the mouse to
# point at intersection points in the plot window we find the following intersection points
# (x,y)=(0.37,-0.43), (0.20,0.56), (1.06,0.13). There are many more. Next we use the fsolve
# command to find solutions to 10 places.
--------------------------------------------------------------------------------
> fsolve({f1,f2},{x,y},{x=0.36..0.38,y=-0.44..-0.42});
Error, (in fsolve) should use exactly all the indeterminates
--------------------------------------------------------------------------------
> fsolve({f1(x,y),f2(x,y)},{x,y},{x=0.36..0.38,y=-0.44..-0.42});
{y = -.4317738060, x = .3728572391}
--------------------------------------------------------------------------------
> fsolve({f1(x,y),f2(x,y)},{x,y},{x=0.19..0.21,y=0.55..0.57});
Error, (in fsolve/genroot) cannot converge to a solution
--------------------------------------------------------------------------------
# What's this? We know there is a solution near (0.20,0.56), which fsolve should be able to
# find. Perhaps using the pointer in the plot window is not too accurate.
--------------------------------------------------------------------------------
> fsolve({f1(x,y),f2(x,y)},{x,y},{x=0.15..0.25,y=0.51..0.59});
{y = .5602946206, x = .2230969126}
--------------------------------------------------------------------------------
# Now we can see why our first attempt did not work; the actual x-coordinate was outside
# the range of x specified. We could repeat this for more intersection points, but you get the
# idea. To be sure we do have solutions we should check the residual values, that is the
# values of f1(x,y) and f2(x,y) at the above solutions. We use tha assign command to turn
# the equations above into assignment of values and then use the print command.
--------------------------------------------------------------------------------
> assign("");
--------------------------------------------------------------------------------
> print(x,y,f1(x,y),f2(x,y),f(x,y));
-9 -9
.2230969126, .5602946206, .2*10 , -.4*10 , .4490307669
--------------------------------------------------------------------------------
>
