#
#
# 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
# x=tan(x) and xcos(x)=sin(x). What are the advantages and disadvantages of
# 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.
--------------------------------------------------------------------------------
