# Lesson 13. Three Equations in Three Unknowns
# =====================================
> restart: with(plots):
# We used the following example in looking at Newton's method in several
# variables. The question was to find all solutions of this system of equations,
# for which "solve" didn't seem to work even though the equations are polynomials.
> f1:= (x,y,z) -> x^2 + x*y*z + y^2*z^3 - 3;
> f2:= (x,y,z) -> x*y^2 - y*z^2 - 2*x^2;
> f3:= (x,y,z) -> x^2*y + y^2*z+ z^4-4;
# The resultant of f1(x,y,z)=0 and f2(x,y,z)=0 with respect to x is a polynomial
# in y and z that is 0 whenever f1(x,y,z) and f2(x,y,z) are both 0 for the same x.
> r1:= unapply(resultant(f1(x,y,z), f2(x,y,z), x), (y,z));
> r2:=unapply(resultant(f1(x,y,z), f3(x,y,z),x), (y,z));
# Since r1 and r2 are functions of two variables, we can plot the curves r1(y,z)=0 and
# r2(y,z)=0 in the plane.
> implicitplot({r1(y,z), r2(y,z)}, y=-10 .. 10, z = -10 .. 10, grid = [200,200]);
# There seem to be several intersections in the region -3 < y < 0,
# 0 < z < 2, a couple with 1.5 < y < 2.5, 0 < z < 1, one near y=2, z=-.3,
# and one near y=3, z=-.9. Let's zoom in on the first region.
> implicitplot({r1(y,z),r2(y,z)}, y=-3 .. 0, z = 0 .. 2, grid=[200,200]);
# There are apparently five intersections here, near
# [-2.14, .80], [-1.44, 1.13], [-1.04, 1.43], [-.81, 1.37] and
# [.67, 1.53]. The "intersection" near [-1.79, .89] is a self-intersection
# of one of the curves.
#
# Let's pick the intersection near [-2.14, .80]. Can we find it more exactly?
> yz:= fsolve({r1(y,z), r2(y,z)}, {y, z}, y=-2.16 .. -2.12, z = .78 .. 0.82);
> f1p:= subs(yz, f1(x,y,z)); f2p:= subs(yz, f2(x,y,z)); f3p:= subs(yz, f3(x,y,z));
> solve(f1p,x);solve(f2p,x);solve(f3p,x);
# So the solution is at x=-.273643496 which (up to roundoff error) is common to
# all three.
# However, not all the intersections of r1=0 and r2=0 correspond to
# solutions of our equations. For example, what about [-1.44, 1.13]?
> yz:= fsolve({r1(y,z), r2(y,z)}, {y, z}, y=-1.46 .. -1.42, z = 1.11 .. 1.15);
> f1p:= subs(yz, f1(x,y,z)); f2p:= subs(yz, f2(x,y,z)); f3p:= subs(yz, f3(x,y,z));
> solve(f1p,x); solve(f2p,x); solve(f3p,x);
# How did this happen? f1(x, -1.448948417, 1.129687620) and
# f2(x, -1.448948417, 1.129687620) have a common root, and so do
# f1(x, -1.448948417, 1.129687620) and f3(x, -1.448948417, 1.129687620),
# it's just not the same root!
#
# Instead of plotting the curves, a more sophisticated method of finding
# the intersections of r1 = 0 and r2 = 0 is to use the resultant of r1 and r2 with
# respect to one variable, say y.
> r3:= unapply(resultant(r1(y,z),r2(y,z),y), z);
# Among the roots of this polynomial are the z values for our intersections.
> s3:= [fsolve(r3(z), z)];
# We have reason to doubt the accuracy of this calculation: catastrophic cancellation is
# almost certain to occur in trying to calculate r3. Let's check
# the accuracy.
> map(r3, s3);
# In the worst case, the result is about 20 orders of magnitude worse than what we
# would like! So let's use 20 more digits.
> Digits:= 30;
> s3:=[fsolve(r3(z), z)];
# Surprise! The solutions for Digits = 10 are actually correct to 10 digits
# (assuming these are correct). It seems that "fsolve" actually uses a larger
# setting of Digits, and then rounds the answers off to obtain an accurate result.
> map( r3, s3);
# Let's just keep 15 digits from now on.
> Digits:= 15;
# Here is a command to take a z value and find the list of y's that satisfy r1=0 with that z.
> gety:= z -> [fsolve(r1(y,z),y)];
> gety(0);
# It would be convenient to put the z values in too, getting a sequence of [y,z] pairs.
> makepairs:= (z,L) -> seq([L[jj],z], jj=1..nops(L));
> makepairs(0, "");
# So here is the list of all pairs [y,z] where r3(z) = 0 and r1(y,z) = 0.
> yzs:= map(z -> makepairs(z,gety(z)), s3);
# We'd like to select those [y,z] pairs for which r2(y,z) = 0. First, here's a function that,
# when applied to a pair [y,z], will return true when r2(y,z) is (almost) 0 and false
# otherwise.
> r2is0:= L -> (evalf(abs(r2(op(L)))) < 10^(-10));
# The "select" command finds those members of a list or set for which a certain
# function returns true (assuming the function always returns true or false).
> yzs:= select(r2is0, yzs);
# Now the corresponding x values, for which f1(x,y,z) = 0.
> getx:= L -> [fsolve(f1(x,op(L)),x)];
> map(L -> makepairs(L, getx(L)), yzs);
# This would be fine except that each member of the list is of the form [x, [y,z]],
# not [x,y,z].
> xyzs:= map(t -> [t[1], op(t[2])], ");
# Now each member of xyzs is a triple [x,y,z] with f1(x,y,z) = 0. We'd like to
# select those for which f2 and f3 are also (almost) 0.
> f2f3are0:= L -> (evalf(abs(f2(op(L))) + abs(f3(op(L)))) < 10^(-10));
> xyzs:= select(f2f3are0, xyzs);
# Here's a less sophisticated method. While resultants are only useful for
# polynomials, the idea of this method can often be used for non-polynomials
# as well.
# Think of f2(x,y,z)=0 as an equation in x and z, for fixed y.
# Complete the square in the x variable.
> 2*(x- y^2/4)^2 + y*z^2 = expand(f2(x,y,z) + 2*(x- y^2/4)^2 + y*z^2) ;
# If y > 0 this is the equation of an ellipse. If y < 0 it is the equation of a hyperbola.
# We can discard y=0 pretty easily:
#
> f1(x,0,z), f2(x,0,z), f3(x,0,z);
# For y > 0 we can parametrize the ellipse:
# x - y^2/4 = y^2/4 cos(t), z = (y/2)^(3/2) sin(t)
> par1:= { x = y^2/4 * (1 + cos(t)), z = (y/2)^(3/2)*sin(t)};
> expand(subs(par1, f2(x,y,z)));
# This makes f1 and f3 into the following functions of y and t.
> f1p:= unapply(subs(par1, f1(x,y,z)), (y,t));
> f3p:= unapply(subs(par1, f3(x,y,z)), (y,t));
> implicitplot({f1p(y,t), f3p(y,t)}, y=0..10, t=0.. 2*Pi);
#
#
> implicitplot({f1p(y,t), f3p(y,t)}, y= 1.4 .. 2.6, t= 0 .. 2.4);
# This indicates one solution around y = 1.68, t = 0.82 and another around y=1.99, t=2.17.
> yt1:= fsolve({f1p(y,t), f3p(y,t)},{y,t}, y=1.66 .. 1.7, t = 0.8 .. 0.84);
> xyz1:= evalf(subs(par1, yt1, [x,y,z]));
# Compare to one of the solutions we found the other way.
> evalm("-xyzs[4]);
> yt2:= fsolve({f1p(y,t), f3p(y,t)},{y,t}, y=1.98 .. 2, t=2.16 .. 2.18);
> yt2:= fsolve({f1p(y,t), f3p(y,t)}, {y,t}, y = 1.95 .. 2.05, t = 2.1 .. 2.2);
> xyz2:= evalf(subs(par1, yt2, [x,y,z]));
> evalm("-xyzs[5]);
> implicitplot({f1p(y,t), f3p(y,t)}, y= 1.9 .. 10, t = 2.4 .. 3.3);
# This doesn't show any solutions (it doesn't eliminate the possibility of some
# solutions for y > 10 though).
#
> implicitplot({f1p(y,t), f3p(y,t)}, y= 10 .. 100, t = 3.1 .. 3.2);
# For y < 0 there are two branches of the hyperbola. In one we can take
# x - y^2/4 = y^2/4 cosh(t), z = (-y/2)^(3/2) sinh(t). Actually we may as well
# take y = -s so s > 0.
> par2:= { x= s^2/4 * (cosh(t)+1), y=-s, z=(s/2)^(3/2)*sinh(t) };
> expand(subs(par2, f2(x,y,z)));
> simplify(");
> f1p2:= unapply(subs(par2, f1(x,y,z)), (s,t));
> f3p2:= unapply(subs(par2, f3(x,y,z)), (s,t));
> implicitplot({f1p2(s,t), f3p2(s,t)}, s=0..10, t=-5 .. 5);
> implicitplot({f1p2(s,t), f3p2(s,t)}, s= 0.3 .. 1.4, t = 1.7 .. 3.6);
> fsolve({f1p2(s,t), f3p2(s,t)}, {s,t}, s= 1 .. 1.2, t = 2 .. 2.2);
> xyz3:= evalf(subs(par2, ", [x,y,z]));
> evalm("-truexyz[3]);
# The other branch of the hyperbola has x - y^2/4 = - y^2/4 cosh(t) and
# z = (-y/2)^(3/2) sinh(t).
> par3:= { x= s^2/4 * (-cosh(t)+1), y=-s, z=(s/2)^(3/2)*sinh(t) };
> f1p3:= unapply(subs(par3, f1(x,y,z)), (s,t));
> f3p3:= unapply(subs(par3, f3(x,y,z)), (s,t));
> implicitplot({f1p3(s,t), f3p3(s,t)}, s=0..10, t=-5 .. 5);
> implicitplot({f1p3(s,t), f3p3(s,t)}, s=0..2, t=1 .. 5);
> fsolve({f1p3(s,t), f3p3(s,t)}, {s,t}, s= 0.7 .. 1, t = 2.1 .. 2.6);
> xyz4:= evalf(subs(par3, ", [x,y,z]));
> evalm("-xyzs[2]);
> implicitplot({f1p3(s,t), f3p3(s,t)}, s=2 .. 10, t=0 .. 1);
> fsolve({f1p3(s,t), f3p3(s,t)}, {s,t}, s= 2 .. 2.4, t = 0.6 .. 0.8);
> xyz5:= evalf(subs(par3, ", [x,y,z]));
> evalm("-xyzs[1]);
