# ASSIGNMENT # 3, SOLUTIONS
#
#
#
# Let f(x)=x^3-2*x^2*y^2+y^3-2 and g(x)=x^2+3*x*y_y^4-2.
--------------------------------------------------------------------------------
# (a) Plot the curves f(x,y)=0 and g(x,y)=0 individually and together for x=-4..4,
# y=-4..4.
--------------------------------------------------------------------------------
> restart;
--------------------------------------------------------------------------------
> f:=x^3-2*x^2*y^2+y^3-2;g:=x^2+3*x*y+y^4-2;
3 2 2 3
f := x - 2 x y + y - 2
2 4
g := x + 3 x y + y - 2
--------------------------------------------------------------------------------
> with(plots):
--------------------------------------------------------------------------------
> F:=implicitplot(f,x=-4..4,y=-4..4):
--------------------------------------------------------------------------------
> G:=implicitplot(g,x=-4..4,y=-4..4):
--------------------------------------------------------------------------------
> display(F,title=`f(x,y)=0`);
--------------------------------------------------------------------------------
#
** Maple V Graphics **
--------------------------------------------------------------------------------
> display(G,title=`g(x,y)=0`);
--------------------------------------------------------------------------------
#
** Maple V Graphics **
--------------------------------------------------------------------------------
> display([F,G],title=`f(x,y)=0 and g(x,y)=0`);
--------------------------------------------------------------------------------
#
** Maple V Graphics **
--------------------------------------------------------------------------------
# (b) How many intersection points do these curves seem to have? You may want to
# "zoom-in" on part of the second quadrant.
--------------------------------------------------------------------------------
# Clicking with the mouse we find the following approximations for the intersection points:
# (3.6, -1.3), (1.3, 0.08) and (-0.23, 1.3). There might be another point of intersection near
# this last one. We "zoom-in" to see what the graph looks like near this point.
--------------------------------------------------------------------------------
> Fzoom:=implicitplot(f,x=-0.3..-0.2,y=1.29..1.33):
--------------------------------------------------------------------------------
> Gzoom:=implicitplot(g,x=-0.3..-0.2,y=1.29..1.33):
--------------------------------------------------------------------------------
> display([Fzoom,Gzoom],title=`a zoom-in`);
--------------------------------------------------------------------------------
#
** Maple V Graphics **
--------------------------------------------------------------------------------
# Yes there are 2 points of intersection here. The coordinates are approximately
# (-0.28, 1.32) and (-0.23, 1.3). Thus there are 4 points of intersection to these graphs.
--------------------------------------------------------------------------------
# (c) Find numerical approximations to the intersection points using fsolve.
--------------------------------------------------------------------------------
> fsolve({f,g},{x,y},x=3.5..3.7,y=-1.4..-1.2);
{x = 3.662985358, y = -1.294642806}
--------------------------------------------------------------------------------
> assign(");
#
--------------------------------------------------------------------------------
> X[1]:=x;Y[1]:=y;x:='x':y:='y':
X[1] := 3.662985358
Y[1] := -1.294642806
--------------------------------------------------------------------------------
> fsolve({f,g},{x,y},x=1.2..1.4,y=0.07..0.09);
Error, (in fsolve/genroot) cannot converge to a solution
--------------------------------------------------------------------------------
# Remember, clicking with the mouse is not too accurate.
--------------------------------------------------------------------------------
> fsolve({f,g},{x,y},x=1.1..1.5,y=0.06..1.1);
{y = .1038732549, x = 1.266919995}
--------------------------------------------------------------------------------
> assign(");
--------------------------------------------------------------------------------
> X[2]:=x;Y[2]:=y;x:='x':y:='y':
X[2] := 1.266919995
Y[2] := .1038732549
--------------------------------------------------------------------------------
> fsolve({f,g},{x,y},x=-0.3..-0.26,y=1.31..1.33);
{y = 1.321820092, x = -.2861212379}
--------------------------------------------------------------------------------
> assign(");
--------------------------------------------------------------------------------
> X[3]:=x;Y[3]:=y;x:='x':y:='y':
X[3] := -.2861212379
Y[3] := 1.321820092
--------------------------------------------------------------------------------
> fsolve({f,g},{x,y},x=-0.26..-0.21,y=1.28..1.34);
{y = 1.300003031, x = -.2334990575}
--------------------------------------------------------------------------------
> assign(");
--------------------------------------------------------------------------------
> X[4]:=x;Y[4]:=y;x:='x':y:='y':
X[4] := -.2334990575
Y[4] := 1.300003031
--------------------------------------------------------------------------------
# We have now found the 4 solutions, assigned them to the names (X[k],Y[k]), and
# preserved x and y as variables.
--------------------------------------------------------------------------------
> for k from 1 to 4 do print(X[k],Y[k]) od;k:='k':
3.662985358, -1.294642806
1.266919995, .1038732549
-.2861212379, 1.321820092
-.2334990575, 1.300003031
--------------------------------------------------------------------------------
# (d) Compute the residuals for the solutions.
--------------------------------------------------------------------------------
> sub:=w->subs(x=X[k],y=Y[k],w);
sub := w -> subs(x = X[k], y = Y[k], w)
--------------------------------------------------------------------------------
# This constructs a substitution operator which we can now use in a loop.
--------------------------------------------------------------------------------
> for k from 1 to 4 do print(sub(f),sub(g)) od;k:='k':
-7 -8
-.11*10 , .4*10
-8 -8
.2265*10 , .13137*10
-8 -8
-.2*10 , -.2*10
-8 -8
-.1*10 , -.1*10
--------------------------------------------------------------------------------
# The residuals are small, but they do not have very many significant digits.
--------------------------------------------------------------------------------
> Digits:=20:
--------------------------------------------------------------------------------
> for k from 1 to 4 do print(sub(f),sub(g)) od;k:='k':
-7 -10
-.219592123062*10 , -.651301577*10
-8 -8
.21492693746339*10 , .108008792221497*10
-8 -8
-.13373037425*10 , -.23044847152*10
-9 -8
-.6427349228*10 , -.10823856120*10
--------------------------------------------------------------------------------
# The residuals are indeed small.
--------------------------------------------------------------------------------
# (e) Show that inverse(J(X,Y))&(f(X,Y),g(X,Y))^{tr} is a reasonable measure of the
# accuracy of the solution.
--------------------------------------------------------------------------------
# This was done in class.
--------------------------------------------------------------------------------
# (f) Using (e) test the accuracy of your solutions.
--------------------------------------------------------------------------------
> with(linalg):
Warning: new definition for norm
Warning: new definition for trace
--------------------------------------------------------------------------------
> J:=jacobian([f,g],[x,y]);
[ 2 2 2 2 ]
[ 3 x - 4 x y - 4 x y + 3 y ]
J := [ ]
[ 3 ]
[ 2 x + 3 y 3 x + 4 y ]
--------------------------------------------------------------------------------
> for k from 1 to 4 do j[k]:=map(sub,J) od;k:='k':
[ 15.694266235982763248 74.511581214615713679 ]
j[1] := [ ]
[ 3.442042298 2.3091528706859105312 ]
[ 4.7605803122698986077 -.63453318334036113788 ]
j[2] := [ ]
[ 2.8454597547 3.8052430105401083312 ]
[ 2.2452497586433064012 4.8087799414123702781 ]
j[3] := [ ]
[ 3.3932178002 8.3796167237481623444 ]
[ 1.7420264187194356437 4.7865095695675190280 ]
j[4] := [ ]
[ 3.4330109780 8.0875642963233167028 ]
--------------------------------------------------------------------------------
# These 4 matrices are the Jacobians at the 4 aprroximations (X[k],Y[k]), k=1..4.
--------------------------------------------------------------------------------
> ?map
--------------------------------------------------------------------------------
> for k from 1 to 4 do r[k]:=vector([sub(f),sub(g)]) od;k:='k':
-7 -10
r[1] := [ -.219592123062*10 , -.651301577*10 ]
-8 -8
r[2] := [ .21492693746339*10 , .108008792221497*10 ]
-8 -8
r[3] := [ -.13373037425*10 , -.23044847152*10 ]
-9 -8
r[4] := [ -.6427349228*10 , -.10823856120*10 ]
--------------------------------------------------------------------------------
# These 4 vectors are the result of substituting the solutions (X[k],Y[k]) into f(x,y) and
# g(x,y), that is the residual vectors.
--------------------------------------------------------------------------------
# Now we multiply the inverse of the Jacobian against these vectors to estimate the errors.
--------------------------------------------------------------------------------
> for k from 1 to 4 do error[k]:=evalm(inverse(j[k])&*r[k]) od;k:='k':
-9 -9
error[1] := [ .20820916078870549575*10 , -.33856350780863766785*10 ]
-9 -10
error[2] := [ .44495636665782432623*10 , -.4888452884406011296*10 ]
-10 -9
error[3] := [ -.497910351745407365*10 , -.2548485161980237262*10 ]
-11 -9
error[4] := [ .73871312565981075*10 , -.13696901491137695527*10 ]
--------------------------------------------------------------------------------
# All solutions appear to be accurate to 10 decimal places. Error[k] is supposed to measure
# the size of the column vector (X[k]-a[k],Y[k]-b[k])^{tr}, k=1..4, where (a[k],b[k]) is an
# exact solution and (X[k],Y[k]) is an approximation.
--------------------------------------------------------------------------------
# (g) Find the exact solutions of the 2 equations using solve.
--------------------------------------------------------------------------------
> Digits:=10:
--------------------------------------------------------------------------------
> S:=solve({f,g},{x,y});
S := {y = %1,
5746853 6 2072275 5 2875722 3 5795527 2 960641 4
x = - ------- %1 - ------- %1 + ------- %1 + ------- %1 - ------- %1
2458382 1229191 1229191 1229191 1229191
118418 9 310480 11 338885 8 773806 7 436411 10
- ------- %1 + ------- %1 + ------- %1 - ------- %1 + ------- %1
1229191 1229191 2458382 1229191 1229191
1529466 363842
+ ------- - ------- %1
1229191 1229191
}
3 4 2 9 11 12 8
%1 := RootOf(74 _Z + 46 _Z + 16 _Z + 4 _Z + 6 _Z + 5 _Z - 13 _Z
7 5 6
- 42 _Z - 26 _Z - 4 + 36 _Z - 34 _Z )
--------------------------------------------------------------------------------
# (h) Use (g) to find numerical approximations to the solutions and then compare these
# with those found above.
--------------------------------------------------------------------------------
> assign(");
--------------------------------------------------------------------------------
> xsol:=x;ysol:=y;x:='x':y:='y':
5746853 6 2072275 5 2875722 3 5795527 2 960641 4
xsol := - ------- %1 - ------- %1 + ------- %1 + ------- %1 - ------- %1
2458382 1229191 1229191 1229191 1229191
118418 9 310480 11 338885 8 773806 7 436411 10
- ------- %1 + ------- %1 + ------- %1 - ------- %1 + ------- %1
1229191 1229191 2458382 1229191 1229191
1529466 363842
+ ------- - ------- %1
1229191 1229191
3 4 2 9 11 12 8
%1 := RootOf(74 _Z + 46 _Z + 16 _Z + 4 _Z + 6 _Z + 5 _Z - 13 _Z
7 5 6
- 42 _Z - 26 _Z - 4 + 36 _Z - 34 _Z )
3 4 2 9 11 12 8
ysol := RootOf(74 _Z + 46 _Z + 16 _Z + 4 _Z + 6 _Z + 5 _Z - 13 _Z
7 5 6
- 42 _Z - 26 _Z - 4 + 36 _Z - 34 _Z )
--------------------------------------------------------------------------------
> p:=op(%1);
3 4 2 9 11 12 8 7
p := 74 _Z + 46 _Z + 16 _Z + 4 _Z + 6 _Z + 5 _Z - 13 _Z - 42 _Z
5 6
- 26 _Z - 4 + 36 _Z - 34 _Z
--------------------------------------------------------------------------------
> ?op;
--------------------------------------------------------------------------------
> Y2:=fsolve(p,_Z);
Y2 := -1.294642806, .1038732549, 1.300003031, 1.321820092
--------------------------------------------------------------------------------
> for k from 1 to 4 do X2[k]:=subs(%1=Y2[k],xsol) od;k:='k':
X2[1] := 3.662985356
X2[2] := 1.266919994
X2[3] := -.2334990535
X2[4] := -.2861212355
--------------------------------------------------------------------------------
> for k from 1 to 4 do print(X2[k],Y2[k]) od;k:='k':
3.662985356, -1.294642806
1.266919994, .1038732549
-.2334990535, 1.300003031
-.2861212355, 1.321820092
--------------------------------------------------------------------------------
# Now we compare with the earlier solutions.
--------------------------------------------------------------------------------
> for k from 1 to 4 do print(X[k],Y[k]) od;k:='k':
3.662985358, -1.294642806
1.266919995, .1038732549
-.2861212379, 1.321820092
-.2334990575, 1.300003031
--------------------------------------------------------------------------------
> for k from 1 to 4 do print(X[k]-X2[k],Y[k]-Y2[k]) od;k:='k':
-8
.2*10 , 0
-8
.1*10 , 0
-.0526221844, .021817061
.0526221780, -.021817061
--------------------------------------------------------------------------------
# I did not notice that the solutions were printed out in a different order.
--------------------------------------------------------------------------------
> for k from 1 to 2 do print(X[k]-X2[k],Y[k]-Y2[k]) od;k:='k':
-8
.2*10 , 0
-8
.1*10 , 0
--------------------------------------------------------------------------------
> print(X[3]-X2[4],Y[3]-Y2[4]); print(X[4]-X2[3],Y[4]-Y2[3]);
-8
-.24*10 , 0
-8
-.40*10 , 0
--------------------------------------------------------------------------------
# Thus the y values are exactly the same and the x values are within 10^{-8}.
--------------------------------------------------------------------------------
# (i) Compute resultant(f,g,x) and resultant(f,g,y) and explain how they can be used with
# part (a) above to find numerical approximations to the solutions as an alternative to (g)
# and (h).
--------------------------------------------------------------------------------
> r:=resultant(f,g,x);
9 5 8 12 11 3 4 7 2
r := 4 y - 26 y - 13 y + 5 y + 6 y + 74 y + 46 y - 42 y + 16 y - 4
6
- 34 y + 36 y
--------------------------------------------------------------------------------
> p;
3 4 2 9 11 12 8 7 5
74 _Z + 46 _Z + 16 _Z + 4 _Z + 6 _Z + 5 _Z - 13 _Z - 42 _Z - 26 _Z
6
- 4 + 36 _Z - 34 _Z
--------------------------------------------------------------------------------
# I thought the polnomial r looked familiar.
--------------------------------------------------------------------------------
> s:=resultant(f,g,y);
12 11 10 9 8 7 6 5
s := 25 x - 72 x - 77 x - 48 x + 211 x + 132 x + 46 x - 116 x
4 3 2
- 242 x + 46 x + 184 x + 72 x + 8
--------------------------------------------------------------------------------
> Y3:=fsolve(r,y);
Y3 := -1.294642806, .1038732549, 1.300003031, 1.321820092
--------------------------------------------------------------------------------
> X3:=fsolve(s,x);
X3 := -.2861212379, -.2334990575, 1.266919995, 3.662985358
--------------------------------------------------------------------------------
# The solutions of r(y)=0 give the y-coordinates of the points of intersection and the
# solutions of s(x)=0 give the x-coordinates. We can not take all possible pairs (X[k],Y[l])
# since this will produce too many pairs, so we go back to the graphs to help us decide
# which x,y coordinates to pair. We see that the correct paiting is (X3[1],Y3[4]),
# (X3[2],Y3[3]), (X3[3],Y3[2]) and (X3[4],Y3[1]).
--------------------------------------------------------------------------------
> for k from 1 to 4 do X3new[k]:=X3[5-k] od;k:='k':
X3new[1] := 3.662985358
X3new[2] := 1.266919995
X3new[3] := -.2334990575
X3new[4] := -.2861212379
--------------------------------------------------------------------------------
> for k from 1 to 4 do print(X3new[k],Y3[k]) od;k:='k':
3.662985358, -1.294642806
1.266919995, .1038732549
-.2334990575, 1.300003031
-.2861212379, 1.321820092
--------------------------------------------------------------------------------
> for k from 1 to 4 do print(X[k],Y[k]) od;k:='k':
3.662985358, -1.294642806
1.266919995, .1038732549
-.2861212379, 1.321820092
-.2334990575, 1.300003031
--------------------------------------------------------------------------------
# These are exactly the same solutions we found earlier, just in a slightly different order.
--------------------------------------------------------------------------------
# (j) Which approach do you prefer?
--------------------------------------------------------------------------------
# This is a matter of personal preference, but we should note that the approach using
# resultants will work only for polynomials, whereas the approach using graphs together
# with fsolve (i.e., (a) and (c)) works more generally.
--------------------------------------------------------------------------------
