# Math 210 Sec. 201
# Solutions to Assignment 1
#
# 1(a)
> eqn:= x^4 - 11*x^3 + 54*x^2 -108*x + 55 = 0;
4 3 2
eqn := x - 11 x + 54 x - 108 x + 55 = 0
> _EnvExplicit:= true;
_EnvExplicit := true
> exactsoln:= solve(eqn, x);
exactsoln :=
1/2 1/2
11/4 + 1/12 3 %2
/ 1/3 1/2 1/2 2/3 1/2 1/2 1/3\1/2
|414 %1 %2 + 36 %2 %1 + 48 %2 + 3258 3 %1 |
+ 1/12 I |-------------------------------------------------------------|
| 1/3 1/2 |
\ %1 %2 /
,
1/2 1/2
11/4 + 1/12 3 %2
/ 1/3 1/2 1/2 2/3 1/2 1/2 1/3\1/2
|414 %1 %2 + 36 %2 %1 + 48 %2 + 3258 3 %1 |
- 1/12 I |-------------------------------------------------------------|
| 1/3 1/2 |
\ %1 %2 /
,
1/2 1/2
11/4 - 1/12 3 %2
/ 1/3 1/2 1/2 2/3 1/2 1/2 1/3\1/2
|- 414 %1 %2 - 36 %2 %1 - 48 %2 + 3258 3 %1 |
+ 1/12 |---------------------------------------------------------------|
| 1/3 1/2 |
\ %1 %2 /
,
1/2 1/2
11/4 - 1/12 3 %2
/ 1/3 1/2 1/2 2/3 1/2 1/2 1/3\1/2
|- 414 %1 %2 - 36 %2 %1 - 48 %2 + 3258 3 %1 |
- 1/12 |---------------------------------------------------------------|
| 1/3 1/2 |
\ %1 %2 /
1/2
%1 := 679/2 + 1/18 37343553
1/3 2/3
- 69 %1 + 12 %1 + 16
%2 := --------------------------
1/3
%1
# Note that "exactsoln" is a sequence. To form a list, we put brackets [] around it.
#
# (b)
> Digits:= 14;
Digits := 14
> evalfsoln:= evalf([exactsoln]);
evalfsoln := [3.6431728008923 + 3.3827531371481 I,
3.6431728008923 - 3.3827531371481 I, 2.9624959932669, .7511584049487]
# (c)
> fsolvesoln:= [fsolve(eqn, x)];
fsolvesoln := [.75115840494860, 2.9624959932669]
# (d) The second answer from (c) was the same as one of the answers from (b), but the
# first was slightly different. To see which was really more accurate (without more
# roundoff complicating matters), it's best to increase Digits some more.
> Digits:= 20;
Digits := 20
> evalfacc:= subs(x=evalfsoln[4], eqn);
-11
evalfacc := -.4568114*10 = 0
> fsolveacc:= subs(x=fsolvesoln[1], eqn);
-12
fsolveacc := -.188170*10 = 0
# It appears that the solution from (c) was more accurate (although the solution from (b)
# was not too far off). The error in (b) is roundoff error from calculating the rather
# complicated exact formula for the roots. The reason for the high accuracy of the "fsolve"
# result turns out to be that "fsolve" internally uses a higher value of Digits, and then
# rounds its result to the original value of Digits.
#
# I'll now return Digits to its default value.
> Digits:= 10:
#
# 2. (a)
> f:= x -> x^3 + exp(1/x);
3
f := x -> x + exp(1/x)
> D(f)(x), (D@@2)(f)(x);
2 exp(1/x) exp(1/x) exp(1/x)
3 x - --------, 6 x + 2 -------- + --------
2 3 4
x x x
# (b) Note that f''(x) > 0 when x > 0. This means that f' is an increasing function.
# So for x > 0 there is at most one interval (say 0 < x < a) where f'(x) < 0 (and f is
# decreasing) and one (a < x < infinity) where f'(x) > 0 and f is increasing. Therefore
# f(x) = c can have at most one solution in the first interval and at most one solution in
# the second interval.
#
# (c) The smallest c is the minimum value of f(x) for x > 0. This will occur at a critical
# point of f, i.e. a solution of f'(x) = 0. The fact that f''(x) > 0 for all x > 0 means that if
# such a critical point exists, the minimum will occur there.
> solve(D(f)(x)=0, x);
# That was worth a try, although I didn't really expect it to work.
> a:= fsolve(D(f)(x)=0, x);
a := .9805089961
> c:= f(a);
c := 3.715516997
# (d) Let's first see what "fsolve" by itself gives us.
> soln1:= fsolve(f(x)=4, x);
soln1 := 1.189115994
# Since this is greater than "a", the other solution should be in the interval 0 .. a.
> soln2:= fsolve(f(x)=4, x, x = 0 .. a);
soln2 := .8006647233
# (e) To check these solutions, I'll evaluate f at each of them, using a higher setting of
# Digits to reduce roundoff error, and compare to 4.
> Digits:= 20:
> f(soln1) - 4;
-9
-.2116296597*10
> f(soln2) - 4;
-11
.33701849*10
# Both of the solutions are about as accurate as you could expect with Digits = 10.
# In fact, each is a "true" solution rounded to 10 significant digits.
> fsolve(f(x)=4, x);
1.1891159940813257727
> fsolve(f(x)=4, x, x = 0 .. a);
.80066472330095859085
> Digits:= 10:
#
# 3 (a) A fixed point is a solution of g(x) = x.
> g:= x -> 1.5 * cos(x);
g := x -> 1.5 cos(x)
> solve(g(x) = x, x);
> fp:= fsolve(g(x) = x, x);
fp := .9148564784
# So there is a fixed point fp at (approximately) .9148564784. To check whether it's an
# attractor or repeller, we check whether |g'(fp)| < 1 or > 1.
> D(g)(fp);
-1.188712591
# Thus it is a repeller.
#
# (b) Several starting points could be used. I'll just try one, starting at x = 1. To save space,
# I won't show the first 20 iterates, just the 21st to 24th.
> X:= 1: for jj from 1 to 20 do X:= evalf(g(X)) od:
> for jj from 21 to 24 do X:= evalf(g(X)) od;
X := .1252851681
X := 1.488243110
X := .1236892224
X := 1.488540354
# It appears to be settling down to a 2-cycle.
#
# (c) We can't really prove that most points will approach the 2-cycle, but we can show
# that the 2-cycle exists and is an attractor. The points on the 2-cycle will
# be solutions of g(g(x)) = x.
> x1:= fsolve(g(g(x))=x, x);
x1 := 1.488653649
# This was lucky: "fsolve" might have given us the fixed point instead. If so, I could have
# asked "fsolve" for a solution in an interval that contained a point of the 2-cycle but not
# the fixed point. Now find the other point of the 2-cycle.
> x2:= g(x1);
x2 := .1230755007
# Now, is the 2-cycle an attractor or repeller? The test is whether |g'(x1) g'(x2)| < 1 or > 1.
> D(g)(x1) * D(g)(x2);
.2752899284
# So it is an attractor. This means that if we start close to x1 or x2, the iterations will
# approach the 2-cycle.
#
# (d) I'll change the definition of g, and again try an experiment starting at x = 1. This
# time no pattern is apparent, even after 100 iterations. I'll just print the 90th to 100th
# iterations. Since I want to draw a cobweb diagram, I'll save the results in a table.
> g:= x -> 2*cos(x);
g := x -> 2 cos(x)
> X[0]:= 1: for jj from 1 to 100 do X[jj]:=evalf(g(X[jj-1])) od:
> seq(X[jj], jj = 90 .. 100);
-.1735617382, 1.969951867, -.7772808092, 1.425646517, .2892813316, 1.916898266,
-.6784669592, 1.557071540, .02744871182, 1.999246616, -.8309233368
# The cobweb diagram is made using the "stair" command that was defined in Lesson 4.
> stair:= x -> plot([[x,x], [x,g(x)], [g(x), g(x)]]);
stair := x -> plot([[x, x], [x, g(x)], [g(x), g(x)]])
> with(plots):
# Note that -2 <= g(x) <= 2 for every x, so x = -2 .. 2 is a reasonable interval to use
# for the diagram.
> p:= plot({x, g(x)}, x = - 2 .. 2):
# The first few steps spiral away from the initial point, which was near the fixed point.
# But up to the 10th step, there's no evidence that the iterations are approaching anything.
> display({p, seq(stair(X[jj]), jj = 0 .. 10) });
# Here are the 80th to 100th steps, still apparently without approaching an attractor. It
# seems that this is an instance of chaos.
> display({p, seq(stair(X[jj]), jj = 80 .. 100) });
# (e) A way to ensure that a cycle is an attractor is to have g'(x) = 0 for some point on the
# cycle. In this case we have g'(0) = 0 (other possibilities are multiples of Pi). To have 0
# on a cycle of period 3, we need (g@@3)(0) = 0. This is an equation which we must
# solve for c. I gave "c" a value back in problem 2, so I must "unassign" it to use it as
# a symbolic variable.
> c:= 'c';
> g:= x -> c * cos(x);
c := c
g := x -> c cos(x)
> eqn:= (g@@3)(0) = 0;
eqn := c cos(c cos(c)) = 0
> c1:= fsolve(eqn, c);
c1 := 0
# Wrong solution: this makes x=0 a fixed point, not a point on a 3-cycle. We need to give
# "fsolve" an interval in which to look. Let's make a graph to see
# what interval would be reasonable.
> plot(lhs(eqn), c = 0 .. 5);
# It looks like the first positive c that works is between 2 and 3.
#
> c2:= fsolve(eqn, c, c = 2 .. 3);
c2 := 2.316112131
>
