# 3.11 Solutions to Homework 2 (Sjerve)

```#
#
#                       HOMEWORK ASSIGNMENT # 2, MATH 210
#
#
> restart;
--------------------------------------------------------------------------------
# QUESTION 1
#
#      Consider the equation x=tan(x) or xcos(x)=sin(x).
#      (a)  Show that there is exactly one root in each interval [mPi,(m+1)Pi] for each
#             integer m>0.
#      (b)  Show that the root in [mPi,(m+1)Pi] is approximately (m+1/2)Pi for m>>0.
#      (c)  For fixed m, let g(x)=mPi+arctan(x).  Show that the iterates x[n]=g(x[n-1])
#             converge to the m^{th} root of x=tan(x) for any choice of starting value x[0].
#      (d)  Find the 100^{th} and 1000^{th} root using (c).
#      (e)  Compare (d) with the use of Maple's fsolve routine for the two equations
#             the 2 approaches?
#
--------------------------------------------------------------------------------
# 1(a)
--------------------------------------------------------------------------------
> plot({x,tan(x)},x=0..5*Pi/2,y=-10..10,discont=true,title=`y=x and y=tan(x)`);
--------------------------------------------------------------------------------
#

** Maple V Graphics **

--------------------------------------------------------------------------------
# The function tan(x) is periodic of period 2*Pi,  has infinite limits at odd multiples of Pi/2
# and is strictly increasing on each interval [(m-1/2)Pi,(m+1/2)Pi].  Thus it is clear from
# the graph that y=x will cross y=tan(x) once in every interval [(m-1/2)Pi,(m+1/2)Pi], m>0,
# and that this crossing point will be just to the left of (m+1/2)Pi.
--------------------------------------------------------------------------------
> limit(tan(x),x=Pi/2,left);limit(tan(x),x=Pi/2,right);

infinity

- infinity
--------------------------------------------------------------------------------
> D(tan);

2
1 + tan
--------------------------------------------------------------------------------
# The derivative of tan(x) is strictly positive.  Another approach is given in the following
# calculations.
--------------------------------------------------------------------------------
> f:=x->tan(x)-x;

f := x -> tan(x) - x
--------------------------------------------------------------------------------
> Diff(f,x)=D(f)(x);;

d            2
---- f = tan(x)
dx
--------------------------------------------------------------------------------
> Limit(f(x),x=(m-1/2)*Pi,right)=-infinity;

Limit                tan(x) - x = - infinity
x -> ((m - 1/2) Pi)+
> Limit(f(x),x=(m+1/2)*Pi,left)=infinity;

Limit                tan(x) - x = infinity
x -> ((m + 1/2) Pi)-
--------------------------------------------------------------------------------
# The function f(x) is increasing from -infinity to +infinity on the interval
# [(m-1/2)Pi,(m+1/2)Pi] and therefore has an unique root somewhere in between.  This
# root must lie in the interval [m*Pi,(m+1/2)Pi] since tan(x) is negative in
# [(m-1/2)Pi,m*Pi].
--------------------------------------------------------------------------------
# 1(b)
--------------------------------------------------------------------------------
# This has already been answered in part (a).  Here is another approach.  Rewrite the
# equation in the form cos(x)=sin(x)/x.
--------------------------------------------------------------------------------
> plot({cos(x),sin(x)/x},x=0..20,title=`y=cos(x) and y=sin(x)/x`);
--------------------------------------------------------------------------------
#

** Maple V Graphics **

--------------------------------------------------------------------------------
# As x tends to infinity the graph of y=cos(x) oscillates infinitely often between 1 and -1
# with period 2*Pi, and the graph of y=sin(x)/x also oscillates infinitely often with the same
# period but with decaying amplitude.  It follows that the solutions of cos(x)=sin(x)/x will
# approximate those of cos(x)=0 as x tends to infinity.
--------------------------------------------------------------------------------
# 1(c)
#
--------------------------------------------------------------------------------
> g:=x->m*Pi+arctan(x);

g := x -> m Pi + arctan(x)
#
--------------------------------------------------------------------------------
# From part (a) we know that there is exactly one solution of x=tan(x) in the interval
# [(m-1/2)Pi,(m+1/2)Pi] for every m>0.  Let's call it a[m].  Then the equation g(x)=x is
# equivalent to m*Pi+arctan(x)=x, i.e. tan(x)=x and x lies in the interval
# [(m-1/2)Pi,(m+1/2)Pi].  Thus if the iteration x[n]=g(x[n-1]) converges it must converge
# to a[m].
--------------------------------------------------------------------------------
> Diff(g,x)=D(g)(x);

d         1
---- g = ------
dx           2
1 + x
--------------------------------------------------------------------------------
> diff(x-g(x),x);

1
1 - ------
2
1 + x
--------------------------------------------------------------------------------
> simplify(");

2
x
------
2
1 + x
--------------------------------------------------------------------------------
# Thus the function h(x)=x-g(x) is strictly increasing for all x.  We already know that
# h(a[m])=0, and therefore h(x)<0 for x<a[m] and h(x)>0 for x>a[m].  Thus if we start
# with x[0]<a[m] then x[0]<g(x[0])=x[1]<g(a[m])=a[m].  By induction we see that if
# x[0]<a[m] then  x[0]<x[1]<x[2]<...<a[m].  That is, the sequence x[n] is an increasing
# sequence bounded above by a[m].  Therefore the limit of x[n] as n tends to infinity exists
# and must equal a[m].  The same argument works for starting values x[0]>a[m].
#
--------------------------------------------------------------------------------
# 1(d)
--------------------------------------------------------------------------------
> g:=x->100*Pi+arctan(x);

g := x -> 100 Pi + arctan(x)
--------------------------------------------------------------------------------
> starttime:=time():
--------------------------------------------------------------------------------
> x[0]:=0.0;

x[0] := 0
--------------------------------------------------------------------------------
> for j from 1 to 5 do x[j]:=evalf(g(x[j-1]))od;j:='j':

x[1] := 314.1592654

x[2] := 315.7268786

x[3] := 315.7268944

x[4] := 315.7268944

x[5] := 315.7268944
--------------------------------------------------------------------------------
> elapsedtime:=time()-starttime;

elapsedtime := .017
--------------------------------------------------------------------------------
# Thus the 100^{th} root is 315.7268944.
--------------------------------------------------------------------------------
> g:=x->1000*Pi+arctan(x);

g := x -> 1000 Pi + arctan(x)
--------------------------------------------------------------------------------
> starttime:=time():
--------------------------------------------------------------------------------
> x[0]:=0.0;

x[0] := 0
--------------------------------------------------------------------------------
> for j from 1 to 5 do x[j]:=evalf(g(x[j-1])) od;j:='j':

x[1] := 3141.592654

x[2] := 3143.163132

x[3] := 3143.163132

x[4] := 3143.163132

x[5] := 3143.163132
--------------------------------------------------------------------------------
> elapsedtime:=time()-starttime;

elapsedtime := .033
--------------------------------------------------------------------------------
# Thus the 1000^{th} root is 3143.163132.
--------------------------------------------------------------------------------
# 1(e)
--------------------------------------------------------------------------------
# To answer this question we time the various approaches.  First x=tan(x).
--------------------------------------------------------------------------------
> starttime:=time():
--------------------------------------------------------------------------------
> fsolve(x=tan(x),x=100*Pi..101*Pi);

fsolve(x = tan(x), x, 100 Pi .. 101 Pi)
--------------------------------------------------------------------------------
# What's this?  We know there is a root in this interval, and yet fsolve did not find it.  The
# reason is that tan(x) is very large near the root.  Perhaps we should narrow the interval of
# search.
--------------------------------------------------------------------------------
> starttime:=time():
--------------------------------------------------------------------------------
> fsolve(x=tan(x),x=100.5*Pi-0.01..100.5*Pi);

315.7268944
--------------------------------------------------------------------------------
> elapsedtime:=time()-starttime;

elapsedtime := .066
--------------------------------------------------------------------------------
# This is quite a bit slower than the iteration in (d).
--------------------------------------------------------------------------------
> starttime:=time():
#
--------------------------------------------------------------------------------
> fsolve(x=tan(x),x=1000.5*Pi-0.001..1000.5*Pi);

3143.163132
--------------------------------------------------------------------------------
> elapsedtime:=time()-starttime;

elapsedtime := .067
--------------------------------------------------------------------------------
# Again quite a bit slower.
--------------------------------------------------------------------------------
# Now we will time the equation xcos(x)=sin(x).
--------------------------------------------------------------------------------
> starttime:=time():
--------------------------------------------------------------------------------
> fsolve(x*cos(x)=sin(x),x=100*Pi..101*Pi);

315.7268944
--------------------------------------------------------------------------------
> elapsedtime:=time()-starttime;

elapsedtime := .083
--------------------------------------------------------------------------------
> starttime:=time():
--------------------------------------------------------------------------------
> fsolve(x*cos(x)=sin(x),x=1000*Pi..1001*Pi);

3143.163132
--------------------------------------------------------------------------------
> elapsedtime:=time()-starttime;

elapsedtime := .050
--------------------------------------------------------------------------------
# QUESTION 2
#
#      You will not need Maple for this exercise.  It suggests a way of speeding up fixed point
# iterations which you may find useful in problems 3 and 4.  We suppose g(x) is a
# differentiable function mapping the interval [a,b] to itself.
#      (a)  For any choice of alpha different from 0 show that h(x)=(1-alpha)x+alphag(x)
#             has the same fixed points as g.
#      (b)  Suppose g(c)=c and g'(c) is not equal to 1.  Show that there is a choice of alpha
#             so that c is an attractive fixed point of h.  What is the best choice of alpha?
#      (c)  What is the range of values of g'(c) for which c is an attractive fixed point of
#             h(x)=(x+g(x))/2?
#
#
#
--------------------------------------------------------------------------------
# 2(a)
--------------------------------------------------------------------------------
> g:='g':
--------------------------------------------------------------------------------
> h(x):=(1-alpha)*x+alpha*g(x);

h(x) := (1 - alpha) x + alpha g(x)
--------------------------------------------------------------------------------
> h(x)-x;

(1 - alpha) x + alpha g(x) - x
--------------------------------------------------------------------------------
> simplify(");

- x alpha + alpha g(x)
--------------------------------------------------------------------------------
# From this it is clear that h(x)=x if, and only if, g(x)=x, since alpha is different from 0.
--------------------------------------------------------------------------------
# 2(b)
--------------------------------------------------------------------------------
> h:=x->(1-alpha)*x+alpha*g(x);

h := x -> (1 - alpha) x + alpha g(x)
--------------------------------------------------------------------------------
> D(h)(c);

1 - alpha + alpha D(g)(c)
--------------------------------------------------------------------------------
# By choosing alpha so that -1<1-alpha+alphaD(g)(c)<1 we see that c is an attractive
# fixed point for h.  The solution is given by |alpha|<2/|D(g)(c)-1|.
--------------------------------------------------------------------------------
> abs(alpha)<2/abs(D(g)(c)-1);

2
abs(alpha) < ----------------
abs(D(g)(c) - 1)
--------------------------------------------------------------------------------
# The best choice of alpha is such that D(h)(c)=0.
--------------------------------------------------------------------------------
> solve(1-alpha+(alpha)*D(g)(c)=0,alpha);

1
- -----------
D(g)(c) - 1
--------------------------------------------------------------------------------
# 2(c)
--------------------------------------------------------------------------------
# For alpha=1/2 we see that c is an attractive fixed point of h if, and only if,
# -1<1/2+1/2D(g)(c)<1.  Solving for D(g)(c) we get -3<D(g)(c)<1.
#
--------------------------------------------------------------------------------
# QUESTION 3
#
#      (a)  Using Maple show that the function g(x)=1.5cos(x) has a fixed point, but that
#             this fixed point is not attractive.
#      (b)  Experiment with the iteration x[n]=g(x[n-1]) for various starting values and
#             then guess the typical behaviour of these iterates.
#      (c)  Prove that your guess is correct by using Maple.
#      (d)  Now put g(x)=2cos(x), try a few hundred iterates of x[n]=g(x[n-1]) from
#             various starting values.  You can use the cobweb program.  What seems to be
#             happening?  You do not have to give a proof.
#
--------------------------------------------------------------------------------
# 3(a)
--------------------------------------------------------------------------------
> g:='g':
--------------------------------------------------------------------------------
> g:=x->1.5*cos(x);

g := x -> 1.5 cos(x)
--------------------------------------------------------------------------------
# All fixed points of g(x) must be in the interval [-1.5,1.5]. so we don't need to look outside
# this interval.
--------------------------------------------------------------------------------
> plot({x,g(x)},x=-1.5..1.5, title=`y=x and y=cos(x)`);
--------------------------------------------------------------------------------
#

** Maple V Graphics **

--------------------------------------------------------------------------------
# Clearly there is exactly one fixed point, and by pointing with the mouse in the plot
# window we find that the x coordinate at the point of intersection is approximately
# x=0.91.
--------------------------------------------------------------------------------
> c:=fsolve(1.5*cos(x)=x,x=0.9..1.0);

c := .9148564784
--------------------------------------------------------------------------------
> D(g)(c);

-1.188712591
--------------------------------------------------------------------------------
# The fixed point is c=0.9148564784, and it is repellant since |g'(c)|>1.
--------------------------------------------------------------------------------
# 3(b)
--------------------------------------------------------------------------------
> stair:=x->plot([[x,x],[x,g(x)],[g(x),g(x)]]);

stair := x -> plot([[x, x], [x, g(x)], [g(x), g(x)]])
> P:=plot({x,g(x)},x=0..2):
> with(plots):
--------------------------------------------------------------------------------
> x[0]:=1;

x[0] := 1
--------------------------------------------------------------------------------
> for j from 1 to 50 do x[j]:=g(x[j-1]) od:j:='j':
--------------------------------------------------------------------------------
> display({P,seq(stair(x[n]),n=1..25)},title=`a cobweb diagram`);
#
--------------------------------------------------------------------------------

** Maple V Graphics **

#
--------------------------------------------------------------------------------
# On the basis of this single experiment my guess is that there is an attractive 2-cycle
# {a,b}, where approximately a=0.13 and b=1.48.
#
--------------------------------------------------------------------------------
# 3(c)
--------------------------------------------------------------------------------
> g2:=g@g;

(2)
g2 := g
--------------------------------------------------------------------------------
> g2(x);

1.5 cos(1.5 cos(x))
--------------------------------------------------------------------------------
> a:=fsolve(g2(x)-x,x=0.12..0.14);

a := .1230755003
--------------------------------------------------------------------------------
> b:=fsolve(g2(x)-x,x=1.4..1.6);

b := 1.488653649
--------------------------------------------------------------------------------
> g(a)-b;

0
--------------------------------------------------------------------------------
> g(b)-a;

-9
.4*10
--------------------------------------------------------------------------------
> D(g2)(a);D(g2)(b);

.2752899275

.2752899284
--------------------------------------------------------------------------------
# Sure enough {a,b} is an attractive 2-cycle.
--------------------------------------------------------------------------------
# 3(d)
--------------------------------------------------------------------------------
> g:='g':
--------------------------------------------------------------------------------
> g:=x->2*cos(x);

g := x -> 2 cos(x)
--------------------------------------------------------------------------------
> stair:=x->plot([[x,x],[x,g(x)],[g(x),g(x)]]);

stair := x -> plot([[x, x], [x, g(x)], [g(x), g(x)]])
--------------------------------------------------------------------------------
> Q:=plot({x,g(x)},x=0..2):
> with(plots):
--------------------------------------------------------------------------------
> x[0]:=1;

x[0] := 1
--------------------------------------------------------------------------------
> for j from 1 to 50 do x[j]:=g(x[j-1]) od:j:='j':
--------------------------------------------------------------------------------
> display({P,seq(stair(x[n]),n=1..25)},title=`a cobweb diagram`);
--------------------------------------------------------------------------------
#

** Maple V Graphics **

--------------------------------------------------------------------------------
# This seems to be an example of chaos.
--------------------------------------------------------------------------------
# QUESTION 4
#
#      Let f(x)=x^4-7*x^2-x-11.  The iterates of of Newton's method, starting at x=1.0
# seem to settle into a cycle of order 4.
#      (a)  Using Maple find the points in this cycle to 10 decimal places of accuracy.  That is,
#             if g(x)=x-f(x)/f'(x), you are to find c[0], c[1], c[2], c[3] to 10 decimal places such
#             that g(c[n])=c[n+1], for n=0,...,3, where c[4]=c[0] and c[0] is near 1.
#      (b)  By computing an appropriate derivative show that this cycle is stable.
#      (c)  Using a cobweb program plot the iterates starting from various initial points for
#             g(x).
#
--------------------------------------------------------------------------------
# 4(a)
--------------------------------------------------------------------------------
> f:='f':g:='g':
--------------------------------------------------------------------------------
> f:=x->x^4-7*x^2-x-11;

4      2
f := x -> x  - 7 x  - x - 11
--------------------------------------------------------------------------------
> g:=x->x-f(x)/D(f)(x);

f(x)
g := x -> x - -------
D(f)(x)
--------------------------------------------------------------------------------
> x[0]:=1.0;

x[0] := 1.0
--------------------------------------------------------------------------------
> for j from 1 to 12 do x[j]:=evalf(g(x[j-1])) od;j:='j':

x[1] := -.636363636

x[2] := 1.258636612

x[3] := -.698857323

x[4] := 1.118361114

x[5] := -.627181486

x[6] := 1.282168826

x[7] := -.722507264

x[8] := 1.073217926

x[9] := -.624277596

x[10] := 1.289784578

x[11] := -.731042757

x[12] := 1.057874941
--------------------------------------------------------------------------------
# This is not too convincing.
--------------------------------------------------------------------------------
> u[0]:=1.0;

u[0] := 1.0
--------------------------------------------------------------------------------
> for j from 1 to 300 do u[j]:=evalf(g(u[j-1])) od:j:='j':
--------------------------------------------------------------------------------
> v[0]:=u[300];

v[0] := 1.061156251
--------------------------------------------------------------------------------
> for j from 1 to 8 do v[j]:=evalf(g(v[j-1])) od;j:='j':

v[1] := -.624890972

v[2] := 1.288168861

v[3] := -.729194361

v[4] := 1.061156264

v[5] := -.624890970

v[6] := 1.288168867

v[7] := -.729194367

v[8] := 1.061156251
--------------------------------------------------------------------------------
> for j from 1 to 16 do v[j]:=evalf(g(v[j-1])) od;j:='j':

v[1] := -.624890972

v[2] := 1.288168861

v[3] := -.729194361

v[4] := 1.061156264

v[5] := -.624890970

v[6] := 1.288168867

v[7] := -.729194367

v[8] := 1.061156251

v[9] := -.624890972

v[10] := 1.288168861

v[11] := -.729194361

v[12] := 1.061156264

v[13] := -.624890970

v[14] := 1.288168867

v[15] := -.729194367

v[16] := 1.061156251
--------------------------------------------------------------------------------
# Actually it is not clear that there is a 4-cycle, perhaps it is an 8-cycle.
--------------------------------------------------------------------------------
> w[0]:=1.061156251;

w[0] := 1.061156251
--------------------------------------------------------------------------------
> for j from 1 to 1000 do w[j]:=evalf(g(w[j-1])) od:j:='j':
--------------------------------------------------------------------------------
> z[0]:=w[1000];

z[0] := 1.061156251
--------------------------------------------------------------------------------
> for j from 1 to 16 do z[j]:=evalf(g(z[j-1])) od;j:='j':

z[1] := -.624890972

z[2] := 1.288168861

z[3] := -.729194361

z[4] := 1.061156264

z[5] := -.624890970

z[6] := 1.288168867

z[7] := -.729194367

z[8] := 1.061156251

z[9] := -.624890972

z[10] := 1.288168861

z[11] := -.729194361

z[12] := 1.061156264

z[13] := -.624890970

z[14] := 1.288168867

z[15] := -.729194367

z[16] := 1.061156251
--------------------------------------------------------------------------------
# Sure enough it is actually an 8-cycle.  Sorry about that.
--------------------------------------------------------------------------------
# 4(b)
--------------------------------------------------------------------------------
> g8:=g@@8;

(8)
g8 := g
--------------------------------------------------------------------------------
>  D(g8)(1.061156251);

D(g)(-.729194367) D(g)(1.288168867) D(g)(-.624890970) D(g)(1.061156264)

D(g)(-.729194361) D(g)(1.288168861) D(g)(-.624890972) D(g)(1.061156251)
--------------------------------------------------------------------------------
# Maple knows the chain rule.
--------------------------------------------------------------------------------
> h(x):=diff(g(x),x);

4      2                2
(x  - 7 x  - x - 11) (12 x  - 14)
h(x) := ---------------------------------
3            2
(4 x  - 14 x - 1)
--------------------------------------------------------------------------------
> product(subs(x=z[j],h(x)),j=1..8);

.1564786891
--------------------------------------------------------------------------------
# Yes this is a stable 8-cycle.
--------------------------------------------------------------------------------

```