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