1.9 Lesson 9

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