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