# 3.1 Assignment 1 Solutions(Israel's)

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

```