# Lesson 9. Solving systems of equations in several variables
# ====================================
> restart;
# Last time, "solve" had found solutions for the following system of two polynomials in two
# variables.
> p1:= (x,y) -> x^4 - x*y^2 - y^3 + 1 ;
> p2:= (x,y) -> x^3 - 5* x^2 * y + y^4 - 2 ;
> sol:=solve({p1(x,y)=0, p2(x,y)=0}, {x,y});
# We now want to get numerical values.
# The values of a "RootOf" can be obtained with the "allvalues" command. This will
# include complex values. If Maple finds an exact expression for the roots, it will use that,
# otherwise it will use numerical approximations.
> allvalues(RootOf(t^3+t^2-1));
# "evalf" will find numerical values for these.
> evalf(");
> allvalues(RootOf(t^5+t^2-1));
> evalf(RootOf(t^3+t^2-1));
# If the expression contains the same RootOf several times, ordinarily "allvalues" does not
# assume the same values should be used. Thus if there are n roots and the RootOf appears
# k times, "allvalues" will give n^k values.
> allvalues(RootOf(t^5+t^2-1) + 1/RootOf(t^5+t^2-1));
# The 'd' option in "allvalues" changes this so that each value uses only one value of the
# RootOf. Thus with this option in the last example we only get 5 answers. Quotes are
# used around the 'd' so this will work even if d has been assigned a value.
> allvalues(RootOf(t^5+t^2-1) + 1/RootOf(t^5+t^2-1),'d');
> v:=allvalues(sol,'d');
# How many of these were there? "nops" counts the number of operands of an expression.
# In the case of a list or set, that's the number of members.
> nops([v]);
# That's what we should expect, because the polynomial in the RootOf had degree 16. Most
# of the solutions are complex, though. We'd like
# to weed those out. Here's a function that decides whether a number is real or not.
#
> isreal:= t -> not(has(t,I));
> vreal:= select(isreal,[v]);
# Alternatively, we could use "fsolve" on the polynomial in the RootOf. This will give us
# the real solutions only, so we won't have to weed out the complex ones. First we have to
# dig out that polynomial.
> P:= op(subs(sol,x));
> xs:= fsolve(P,_Z);
# These are, perhaps not surprisingly, exactly the same x values as in "vreal". Now to get
# the solutions. We substitute each of these x values for the RootOf in "sol".
> subs(RootOf(P)=xs[1], sol);
> subs(RootOf(P)=xs[2],sol);
# Again, exactly the same solutions as we got from "allvalues". We could get both at once
# using the "map" command, which performs some function on all members of a set or list.
# In this case the function to perform is what we did above with xs[1] and xs[2].
> vs:=map(t -> subs(RootOf(P)=t,sol), [xs]);
# Now let's check the accuracy of these solutions.
> map(t -> subs(t, [p1(x,y), p2(x,y)]), vs);
> subs(vs[1], [p1(x,y), p2(x,y)]);
# Not really very good accuracy, especially for the second one. Look at the terms in the
# expression for y in the solution.
> convert(subs(sol,y), list);
> subs(RootOf(P)=xs[2],");
# Catastrophic cancellation is the culprit. We could do better by using a higher setting of
# Digits. The largest terms in the sum are about 5 orders of magnitude higher than the
# final answer, so we should use at least 5 extra digits. To be safer, use 6.
#
> Digits:= 16;
> xs:= fsolve(P,_Z);
> vs:=map(t -> subs(RootOf(P)=t,sol), [xs]);
# Note that in our solutions with Digits = 10, the x values were correct to 10 digits, but the
# y values weren't (expecially in the second solution where the 5th digit was already wrong).
> map(t -> subs(t, [p1(x,y), p2(x,y)]), vs);
# Of course, we could have used "fsolve" on the original polynomials. But we'd only get
# one solution (it's only for a single polynomial, not a set of them, that "fsolve" produces all
# real solutions).
> Digits:= 10: fsolve({p1(x,y),p2(x,y)},{x,y});
# A picture would be useful. The "implicitplot" command in the "plots" package plots a
# curve defined by an implicit formula such as p1(x,y)=0. Or we could do several at once
# (putting the formulas in a set).
> with(plots):
> implicitplot({p1(x,y)=0, p2(x,y)=0}, x=-4 .. 4, y = -4 .. 4);
# Note our two solutions, at the points where the curves cross.
#
# Which curve is which? On a colour monitor, they are in different colours, but it's hard to
# know which curve is which colour (because Maple can change the order of a set). We
# could produce two plots separately, assigning a colour to each, and combine them with
# "display".
> plot1:= implicitplot(p1(x,y)=0, x=-4 .. 4, y=-4 .. 4, colour=red):
> plot2:= implicitplot(p2(x,y)=0, x=-4 .. 4, y= -4 .. 4, colour=blue):
> display({plot1,plot2});
> p1(x,y), p2(x,y);
#
# Non-Polynomial Systems
# ===================
#
# Let's look at the critical points of the following function.
> f:= (x,y) -> sin(y-x^2-1) + cos(2*y^2-x);
# First we try to visualize the graph z = f(x,y). This
# is a three-dimensional graph.
> plot3d(f(x,y), x=-2 .. 2, y=-2 .. 2, axes=BOXED);
> contourplot(f(x,y), x=-3 .. 3, y = -3 .. 3, grid=[50, 50], axes=BOXED, shading=ZHUE);
# The critical points are the points where the partial
# derivatives of f with respect to x and y are 0.
#
> eq1:= D[1](f)(x,y)=0; eq2:= D[2](f)(x,y)=0;
# Do you suppose "solve" will be able to make
# anything out of this?
#
> solve({eq1, eq2},{x,y});
> s1:= fsolve({eq1,eq2},{x,y});
# There should be lots of others, of course.
> s2:= fsolve({eq1,eq2},{x,y},x=2.6 .. 3, y=-2 .. -1.6);
