# 3.2 Assignment 2 Solutions(Israel's)

```# Math 210 Sec. 201
# Solutions to Assignment 2
#
# 1.(a) Because typically x[n+1] - r is approximately c (x[n]-r)^2 for some c.
#
#    (b)
> Digits:= 80: x[0]:= 2: f:= x -> 1 + sin(2*x)-x:
> newt:= unapply(x - f(x)/D(f)(x), x);

1 + sin(2 x) - x
newt := x -> x - ----------------
2 cos(2 x) - 1
> for nn from 1 to 10 do x[nn]:= evalf(newt(x[nn-1])) od:
# An indication of the number of correct digits is log[10]|x[n]-x[10]|.
# If this is between -k and -k+1, i.e. 10^(-k) < |x[n]-x[10]| < 10^(-k+1),
# then at most k digits (including the one before the decimal point) are correct.
# It could be less.  We can tell by looking at the first k digits, which are the
# integer part of x[n]*10^(k-1).  The "floor" function takes the integer part
# of a number.
#
# Of course you could do all this just by counting, but it wouldn't be as much fun!
> floor(log[10](abs(x[1]-x[10])));

-1
> floor(x[1])=floor(x[10]);

1 = 1
> floor(log[10](abs(x[2]-x[10])));

-3
> trunc(x[2]*10^2)=trunc(x[10]*10^2);

138 = 137
> floor(log[10](abs(x[3]-x[10])));

-5
> trunc(x[3]*10^4)=trunc(x[10]*10^4);

13773 = 13773
> floor(log[10](abs(x[4]-x[10])));

-11
> trunc(x[4]*10^10)=trunc(x[10]*10^10);

13773368771 = 13773368771
> floor(log[10](abs(x[5]-x[10])));

-21
> trunc(x[5]*10^20)=trunc(x[10]*10^20);

137733687712314295169 = 137733687712314295169
> floor(log[10](abs(x[6]-x[10])));

-43
> trunc(x[6]*10^42)=trunc(x[10]*10^42);

1377336877123142951690481942218603304051598 =

1377336877123142951690481942218603304051597
> floor(log[10](abs(x[7]-x[10])));
Error, (in ln) singularity encountered

> x[7]-x[10];

0
# Only the first digit is correct in x[1], the first two in x[2], five in x[3], 11 in x[4], 21 in
# x[5], 42 in x[6].  From x[7] on we have x[n+1] = x[n] to the given accuracy.
#
# (c) Because the proof of the error estimate depended
# on f' being nonzero at the root.  Or because when f' is
# near 0, a small change in f(x) can correspond to a much larger change in x.
> Digits:= 10:
> f:= x -> (x-1)^3: newt:= unapply(x - f(x)/D(f)(x), x);

newt := x -> 2/3 x + 1/3
> x[0]:= 2:
> normal((newt(x)-1)/(x-1));

2/3
# So we should have x[n]-1 = A (2/3)^n for some constant A.  For n=0 we want A =
# x[0]-1 = 1, so
# the formula is x[n] = 1 + (2/3)^n.
> for nn from 1 to 10 do x[nn]:= newt(x[nn-1]) od:
> seq(x[nn]=1+(2/3)^nn, nn=1..5);

35     35    97     97   275   275
5/3 = 5/3, 13/9 = 13/9, ---- = ----, ---- = ----, --- = ---
27     27    81     81   243   243
# Multiplying the error by 2/3 at each stage, the number of correct digits will usually stay
# the same, but sometimes increase by 1.
> solve((2/3)^n = 1/10, n);

ln(10)
- -------
ln(2/3)
> evalf(");

5.678873587
# So we can expect one more digit every 5 or 6 iterations.  The root is correct to 10 decimal
# digits,
# i.e. the error is less than .5 * 10^(-10), if n is at least:
> evalf(ln(.5*10^(-10))/ln(2/3));

58.49824715
> evalf(1+(2/3)^59, 15);

1.00000000004080
# So we need n to be at least 59.
#
# (d)
> f:= x -> x^(1/x): x[0]:= 0.5:
> newt:= unapply(x - f(x)/D(f)(x), x);

1
newt := x -> x - --------------
ln(x)     1
- ----- + ----
2      2
x      x
> for nn from 1 to 30 do x[nn]:= evalf(newt(x[nn-1])) od:
> x[19], x[20];

.1014982263, .09836477550
# The first one less than 1/10 is x[20].
#
# 2(a).  We might try to get newt(0)=1, newt(1)=2, and newt(2)=0 for the 3-cycle.
# To ensure that the cycle is stable, we can require D(newt)(0)=0.  This will be
# true if (D@@2)(f)(0)=0.
> f:= 'f':
> newt:= x -> x - f(x)/D(f)(x);

f(x)
newt := x -> x - -------
D(f)(x)
> D(newt)(0);

(2)
f(0) D   (f)(0)
---------------
2
D(f)(0)
# We might try a cubic polynomial, but that turns out not to work.  How about
# a quartic?
> f:= x -> a*x^4 + b*x^3 + c*x + d;

4      3
f := x -> a x  + b x  + c x + d
# There's no x^2 term so that (D@@2)(f)(0)=0.
> equations:= {newt(0)=1, newt(1)=2, newt(2)=0};

a + b + c + d          16 a + 8 b + 2 c + d
equations := {- d/c = 1, 1 - ------------- = 2, 2 - -------------------- = 0}
4 a + 3 b + c             32 a + 12 b + c
> solve(equations);

43
{d = d, a = - 3/28 d, b = --- d, c = - d}
112
# We could take any nonzero value for d.  To make all coefficients integers,
# try d=112.
> subs(", f(x));

4    43    3
- 3/28 d x  + --- d x  - d x + d
112
> subs(d=112,");

4       3
- 12 x  + 43 x  - 112 x + 112
> f:= unapply(",x);

4       3
f := x -> - 12 x  + 43 x  - 112 x + 112
# (b) We can be sure that all points in the interval (-a,a) are attracted to the
# 3-cycle if |g(x)/x| < 1 for 0 < |x| < a, where g(x) = newt(newt(newt(x))).
> g:= unapply(newt(newt(newt(x))), x):
# This is an extremely complicated function.  Its derivative is even more complicated.  I
# crashed Maple on my PC several times trying to use "fsolve"
# with it.  So I'll have to be careful.
> plot({g(x)/x, -1, 1}, x = -.1 .. .1);
# The plot shows g(x)/x = 1 at some point near -.037, g(x)/x = -1 at some point
# near +.04, and g(x)/x between 1 and -1 for x between these points.
> g(-.037)/(-.037);

1.019476622
> g(-.036)/(-.036);

.9885820001
# So the interval (-.036, .036) will do.  More detailed investigation shows that
# (-.03637, .03779) is a suitable interval.
#
# 3(a)
#
> eq1:= x^3 + 3 *x *y^3 - y = 0 ;

3        3
eq1 := x  + 3 x y  - y = 0
> eq2:= x^2 *y^2 - 2 *x^3 - 3 *y^3 = 4;

2  2      3      3
eq2 := x  y  - 2 x  - 3 y  = 4
> with(plots, implicitplot):
> implicitplot({eq1, eq2}, x=-5 .. 5, y = -5 .. 5);
# My mouse locates the two intersections at about
# [-1.137, -.6161] and [.2243, -1.111].
#
# (b)  Note the use of ranges for x and y to help "fsolve"
# find the solutions we want.  However, we want to be careful not to make the ranges too
# narrow, because the
# mouse's locations for the intersections are not very
# accurate.
> sol1:= fsolve({eq1,eq2},{x,y},{x=-1.3 .. -1, y=-.8 .. -.5});

sol1 := {x = -1.120804224, y = -.6172428838}
> sol2:= fsolve({eq1,eq2},{x,y},{x=0.1 .. 0.4, y = -1.3 .. -1});

sol2 := {x = .2832771202, y = -1.095952011}
# (c)
> res:= resultant(lhs(eq1), lhs(eq2)-4, x);

13        12       11       10        9        8       7        6
res := - 27 y   - 324 y   - 54 y   - 36 y   - 477 y  - 143 y  - 42 y  - 132 y

5        4        3       2
- 132 y  - 144 y  - 152 y  - 48 y  - 96 y - 64
> ys:= [fsolve(res,y)];

ys := [-11.83006098, -1.095952011, -.6172428838]
> subs(y=ys[1], eq1);

3
x  - 4966.863267 x + 11.83006098 = 0
> fsolve(");

-70.47716739, .002381797194, 70.47478560
> subs(y=ys[1], eq2);

2      3
139.9503428 x  - 2 x  + 4966.863267 = 4
> fsolve(");

70.47478562
# Thus there is a solution of the equations at
# [70.4747856, -11.83006098].   This is a surprise because it is far outside the range of our
# plot, which gave little hint that there was a solution in this direction.
# The other two solutions are the ones we know.
> fsolve(subs(y=ys[2], eq1),x);

-2.113666698, .2832771200, 1.830389578
> fsolve(subs(y=ys[2], eq2), x);

-.1805423016, .2832770998, .4978206068
# This finds the solution at approximately [.2832771, -1.095952011].
> fsolve(subs(y=ys[3],eq1),x);

-1.120804224
> fsolve(subs(y=ys[3],eq2),x);

-1.120804224
# This finds the solution at [-1.120804224, -.6172428838].
# (d)  A higher setting of Digits may be needed because there is likely to be some
# catastrophic cancellation in the calculations.
> Digits:= 20:
> soln:=solve({eq1,eq2},{x,y});

soln := {x = %1,

30739345564      4367504387   8   40237436898   4   32679874755   3
y = - ----------- %1 + ---------- %1  - ----------- %1  + ----------- %1
7794397861      7794397861        7794397861        7794397861

40304307    12   2919034588   11   16973544208   10   8309488874   9
- ---------- %1   + ---------- %1   - ----------- %1   + ---------- %1
7794397861        7794397861        23383193583        7794397861

21068161718   7   30183552660   6    2750028508   5    6931063024   2
- ----------- %1  + ----------- %1  - ----------- %1  + ----------- %1
7794397861        7794397861       23383193583       23383193583

514800932
- ----------
7794397861

}

13        12         11        10         9         8        7
%1 := RootOf(_Z   - 72 _Z   + 108 _Z   - 30 _Z   - 429 _Z  + 427 _Z  - 12 _Z

6         5        4         3
- 876 _Z  + 428 _Z  + 96 _Z  - 570 _Z  + 12)
> ro:= subs(soln, x);

13        12         11        10         9         8        7
ro := RootOf(_Z   - 72 _Z   + 108 _Z   - 30 _Z   - 429 _Z  + 427 _Z  - 12 _Z

6         5        4         3
- 876 _Z  + 428 _Z  + 96 _Z  - 570 _Z  + 12)
> xs:= [fsolve(op(ro), _Z)];

xs := [-1.1208042240180775282, .28327712023380441512, 70.474785577010416057]
> sol1:= subs(ro=xs[1], soln);

sol1 := {x = -1.1208042240180775282, y = -.61724288377170431638}
> sol2:= subs(ro=xs[2], soln);

sol2 := {x = .28327712023380441512, y = -1.0959520105979225755}
> sol3:= subs(ro=xs[3], soln);

sol3 := {x = 70.474785577010416057, y = -12.254689243523496846}
# Note that the y value is quite different from what I found in (c).  That would have been
# improved by using the higher setting of Digits.
> Digits:= 10:
# (e) We use the several-variables Newton's method
# from Lesson 12.
> with(linalg):
Warning: new definition for   norm
Warning: new definition for   trace

> f1:= unapply(lhs(eq1),(x,y));

3        3
f1 := (x,y) -> x  + 3 x y  - y
> f2:= unapply(lhs(eq2)-4,(x,y));

2  2      3      3
f2 := (x,y) -> x  y  - 2 x  - 3 y  - 4
> F:= V -> [f1(op(V)), f2(op(V))];

F := V -> [f1(op(V)), f2(op(V))]
> xy:= [x,y]:
> jacobian(F(xy),xy);

[     2      3          2      ]
[  3 x  + 3 y      9 x y  - 1  ]
[                              ]
[      2      2     2        2 ]
[ 2 x y  - 6 x   2 x  y - 9 y  ]
> jlist:= unapply(convert(",listlist),(x,y));

2      3       2             2      2     2        2
jlist := (x,y) -> [[3 x  + 3 y , 9 x y  - 1], [2 x y  - 6 x , 2 x  y - 9 y ]]
> j:= V -> matrix(jlist(op(V)));

j := V -> matrix(jlist(op(V)))
> evl:= V -> convert(evalm(V),list);

evl := V -> convert(evalm(V), list)
> evl(F([x0,y0]) + j([x0,y0]) &* (xy-[x0,y0]));

3          3       2         3            2
[- 2 x0  - 9 x0 y0  + 3 x0  x + 3 y0  x + 9 x0 y0  y - y,

2   2       3       3              2         2         2            2
- 3 x0  y0  + 4 x0  + 6 y0  - 4 + 2 x0 y0  x - 6 x0  x + 2 x0  y0 y - 9 y0  y

]
> linapp:= unapply(",(x0,y0,z0));

linapp := (x0,y0,z0) -> [

3          3       2         3            2
- 2 x0  - 9 x0 y0  + 3 x0  x + 3 y0  x + 9 x0 y0  y - y,

2   2       3       3              2         2         2            2
- 3 x0  y0  + 4 x0  + 6 y0  - 4 + 2 x0 y0  x - 6 x0  x + 2 x0  y0 y - 9 y0  y

]
> newt:= V0 -> evalf(subs(solve({op(linapp(op(V0)))},{x,y}),[x,y]));

newt := V0 -> evalf(subs(solve({op(linapp(op(V0)))}, {x, y}), [x, y]))
> V:= [1,1];

V := [1, 1]
> for jj from 1 to 12 do V:= newt(V) od;

V := [5.300000000, -2.600000000]

V := [3.452098609, -2.020542085]

V := [2.298710571, -1.583785834]

V := [1.500462876, -1.253562819]

V := [.8158876771, -1.015508260]

V := [.1232960249, -.9832738501]

V := [.3457538237, -1.109753445]

V := [.2849522878, -1.095959586]

V := [.2832765774, -1.095951893]

V := [.2832771208, -1.095952009]

V := [.2832771203, -1.095952011]

V := [.2832771200, -1.095952011]
# (f)
#
> newt([x0,y0]);

5         4   2        3        3   4         3   2        2   2
[(4. y0 x0  + 36. x0  y0  - 4. x0  - 9. x0  y0  - 18. x0  y0  + 3. x0  y0

5            2        3         /       4            2   2
- 27. x0 y0  - 36. x0 y0  - 6. y0  + 4.)  /  (6. x0  y0 - 27. x0  y0
/

4   2         5           2         3   2        2
- 12. y0  x0  - 27. y0  + 2. x0 y0  + 54. x0  y0  - 6. x0 ),

4   2         2   3         2         6         3        5   2
(5. x0  y0  - 18. x0  y0  + 12. x0  - 18. y0  + 12. y0  - 9. y0  x0

3   3    /       4            2   2         4   2         5
+ 42. y0  x0 )  /  (6. x0  y0 - 27. x0  y0  - 12. y0  x0  - 27. y0
/

2         3   2        2
+ 2. x0 y0  + 54. x0  y0  - 6. x0 )                                ]
> eqs:= { "[1] = 70.47478558, "[2] = -12.25468924 };

eqs := {

5         4   2        3        3   4         3   2        2   2
(4. y0 x0  + 36. x0  y0  - 4. x0  - 9. x0  y0  - 18. x0  y0  + 3. x0  y0

5            2        3         /       4            2   2
- 27. x0 y0  - 36. x0 y0  - 6. y0  + 4.)  /  (6. x0  y0 - 27. x0  y0
/

4   2         5           2         3   2        2
- 12. y0  x0  - 27. y0  + 2. x0 y0  + 54. x0  y0  - 6. x0 ) = 70.47478558,

4   2         2   3         2         6         3        5   2
(5. x0  y0  - 18. x0  y0  + 12. x0  - 18. y0  + 12. y0  - 9. y0  x0

3   3    /       4            2   2         4   2         5
+ 42. y0  x0 )  /  (6. x0  y0 - 27. x0  y0  - 12. y0  x0  - 27. y0
/

2         3   2        2
+ 2. x0 y0  + 54. x0  y0  - 6. x0 ) = -12.25468924                 }
> fsolve(eqs, {x0,y0}, {x0=-1..1, y0=-1..1});
Error, (in fsolve/genroot) cannot converge to a solution

# We may be running into trouble with the denominator hitting 0.
> eqs[1];

5         4   2        3        3   4         3   2        2   2
(4. y0 x0  + 36. x0  y0  - 4. x0  - 9. x0  y0  - 18. x0  y0  + 3. x0  y0

5            2        3         /       4            2   2
- 27. x0 y0  - 36. x0 y0  - 6. y0  + 4.)  /  (6. x0  y0 - 27. x0  y0
/

4   2         5           2         3   2        2
- 12. y0  x0  - 27. y0  + 2. x0 y0  + 54. x0  y0  - 6. x0 ) = 70.47478558
> numer(lhs(")-rhs("));

5         4   2        3        3   4                 3   2
4. y0 x0  + 36. x0  y0  - 4. x0  - 9. x0  y0  - 3823.638421 x0  y0

2   2            5                    2        3
+ 1905.819211 x0  y0  - 27. x0 y0  - 176.9495712 x0 y0  - 6. y0  + 4.

4                    4   2                 5
- 422.8487135 x0  y0 + 845.6974270 y0  x0  + 1902.819211 y0

2
+ 422.8487135 x0
> neweqs:= { ", numer(lhs(eqs[2])-rhs(eqs[2]))};

neweqs := {

4   2         2   3                 2         6         3        5   2
5. x0  y0  - 18. x0  y0  - 61.52813544 x0  - 18. y0  + 12. y0  - 9. y0  x0

3   3                 4                    2   2
+ 42. y0  x0  + 73.52813544 x0  y0 - 330.8766095 x0  y0

4   2                 5                    2
- 147.0562709 y0  x0  - 330.8766095 y0  + 24.50937848 x0 y0

3   2
+ 661.7532190 x0  y0 ,

5         4   2        3        3   4                 3   2
4. y0 x0  + 36. x0  y0  - 4. x0  - 9. x0  y0  - 3823.638421 x0  y0

2   2            5                    2        3
+ 1905.819211 x0  y0  - 27. x0 y0  - 176.9495712 x0 y0  - 6. y0  + 4.

4                    4   2                 5
- 422.8487135 x0  y0 + 845.6974270 y0  x0  + 1902.819211 y0

2
+ 422.8487135 x0                                                     }
> fsolve(neweqs, {x0,y0}, {x0=-1..1, y0=-1..1});

{y0 = -.4545442561, x0 = .2488423365}
> newt(subs(",[x0,y0]));

[70.47478557, -12.25468923]
# 4. (a) Critical points of f in the region:
>  f:= (x,y) -> x^2 -x*y + y^3 - x:\
g1:= (x,y) -> x^2 + x*y + y^2 - 3:\
g2:= (x,y) -> x^2 + y^2 - 4:\
\

> eqa:= { diff(f(x,y),x)=0, diff(f(x,y),y)=0 };

2
eqa := {- x + 3 y  = 0, 2 x - y - 1 = 0}
> cpa:= solve(eqa, {x,y});

cpa := {y = -1/3, x = 1/3}, {y = 1/2, x = 3/4}
> subs(cpa[1], [g1(x,y), g2(x,y)]);

[-26/9, -34/9]
> subs(cpa[2],[g1(x,y), g2(x,y)]);

29      51
[- ----, - ----]
16      16
# Both of the critical points are in the region g1(x,y) <= 0, g2(x,y) <= 0.  The values at
# these critical points:
> valuesa:= evalf([subs(cpa[1], f(x,y)),  subs(cpa[2], f(x,y))]);

valuesa := [-.1481481481, -.4375000000]
# (b) Critical points of restriction of f(x,y) to g2(x,y) = 0 with g1(x,y) <= 0.
>  L:= (x,y,t) -> f(x,y) + t*g2(x,y):\

>  eqb:= { diff(L(x,y,t),x)=0, diff(L(x,y,t),y)=0,\
diff(L(x,y,t),t)=0 };\
\

2    2                                            2
eqb := {x  + y  - 4 = 0, 2 x - y - 1 + 2 t x = 0, - x + 3 y  + 2 t y = 0}
> cpb:= solve(eqb, {x,y,t});\
\

15       15    4   135   2   193   3   279   5    89
cpb := {y = %1, t = - --- %1 + ---- %1  - --- %1  + --- %1  - --- %1  - ---,
704       88        176       176       704       176

2    63    5    67    3          4   151       31
x = 5/22 %1  - ---- %1  + ---- %1  - 3/11 %1  - --- %1 + ----}
88         22                    88       22

6        5        4        3        2
%1 :=   RootOf(9 _Z  - 12 _Z  - 28 _Z  + 52 _Z  - 31 _Z  - 8 _Z + 16)
>  ro:= subs(sols, y);

6        5        4        3        2
ro := RootOf(9 _Z  - 12 _Z  - 28 _Z  + 52 _Z  - 31 _Z  - 8 _Z + 16)
>  ys:= [fsolve(op(ro), _Z)];

ys := [-1.996168431, -.5578076195, .9336779718, 1.866125286]
> cpb:= map(r -> subs(ro=r, cpb), ys);

cpb := [{t = 3.025247044, x = -.123740028, y = -1.996168431},

{x = 1.920638086, y = -.5578076195, t = -.8848839916},

{y = .9336779718, t = -.4533570598, x = 1.768684665},

{t = -2.991947713, x = -.719427820, y = 1.866125286}]
> map(S -> subs(S, g1(x,y)), cpb);

[1.247005937, -.071346562, 2.651381910, -.342542475]
# The second and fourth solutions are in our region.
> cpb:= [cpb[2], cpb[4]];

cpb := [{x = 1.920638086, y = -.5578076195, t = -.8848839916},

{t = -2.991947713, x = -.719427820, y = 1.866125286}]
# Here are the corresponding values of f.
> valuesb:= [subs(cpb[1], f(x,y)), subs(cpb[2], f(x,y))];

valuesb := [2.665997657, 9.078185359]
# (c) The intersections of the curves g1(x,y)=0 and g2(x,y) = 0.
> sols:= solve({g1(x,y)=0, g2(x,y)=0}, {x,y});

4       2
sols := {x = RootOf(_Z  - 4 _Z  + 1),

4       2                 4       2     3
y = - 4 RootOf(_Z  - 4 _Z  + 1) + RootOf(_Z  - 4 _Z  + 1) }
> ro:= subs(sols, x);

4       2
ro := RootOf(_Z  - 4 _Z  + 1)
> xs:= [fsolve(op(ro),_Z)];

xs := [-1.931851653, -.5176380902, .5176380902, 1.931851653]
> inters:= map(t -> subs(ro=t, sols), xs);

inters := [{x = -1.931851653, y = .517638087},

{x = -.5176380902, y = 1.931851653}, {x = .5176380902, y = -1.931851653},

{x = 1.931851653, y = -.517638087}]
# There are four intersections.  Here are the corresponding f values.
> valuesc:= map(s -> subs(s, f(x,y)), inters);

valuesc := [6.802603162, 8.995355807, -6.459457423, 2.661498444]
# In Lesson 14 we found one critical point for the restriction of f(x,y) to  g1(x,y) = 0 with
# g2(x,y) <= 0.
>     cp14:=    {x = 1.213451017, y = .7701014499, t = -.2054425610}:
>     value14:=  -.2187545609;

value14 := -.2187545609
# So here are all our values:
> values:= [ op(valuesa), op(valuesb), op(valuesc), value14];

values := [-.1481481481, -.4375000000, 2.665997657, 9.078185359, 6.802603162,

8.995355807, -6.459457423, 2.661498444, -.2187545609]
# The maximum value is 9.078185359, which was from the second critical point in (b).
> maxpt:= subs(cpb[2], [x,y]);

maxpt := [-.719427820, 1.866125286]
> f(op(maxpt));

9.078185359
# The minimum value is -6.459457423, which was from the third intersection point in (c).
> minpt:= subs(inters[3], [x,y]);

minpt := [.5176380902, -1.931851653]
> f(op(minpt));

-6.459457423
>

```