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