# 3.12 Solutions to Homework 3 (Sjerve)

```# 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.
--------------------------------------------------------------------------------

```