-------------------------------------------------------------------
Math 1225                 Laboratory Assignment  #5             March
10, 1999
Treibergs                                                   Due March 
24, 1999
--------------------------------------------------------------------
BISECTION METHOD
This search algorithm doesn't make use of derivatives. Instead it uses
the idea that if  y=f(x)  is continuous and there are two numbers  a<b

where, say, f(a)<0 and f(b)>0 then somewhere in between there is a  c,
a<c<b
so that f(c)=0. (Intermediate value theorem.) The idea of bisection is
that
one replaces the interval [a.0,b.0]  by one that is half as long
[a.1,b.1]
 where either a.1=a.0  and  b.1=(a.0+b.0)/2  or  a.1=(a.0+b.0)/2 and 
b.1=b.0. The choice depends on the sign of  f((a.0+b.0)/2).  If it is
positive,
make the first replacement. If negative, make the second. For example
lets
compute the square root of  5.
> f:=t-> 5-t^2;`Actual square rot of five =`,sqrt(5.0);
Then guess an interval, here  a=1, b=3.  f(a)<0 and f(b)>0.
> a.0:=1.00000000;b.0:=3.00000000;f(a.0)*f(b.0);
Note that if  f(a)*f(b)<0 we may proceed in [a,b]. After n steps we
have
divided the interval in  half n times. Thus the number (a.n+b.n)/2 is
within
(b.0-a.0)*(0.5)^(n+1)  of the zero.
> n:=11; error:= (b.0-a.0)*(0.5)^(n+1);
We make a do loop that stores the endpoints of the intervals. The 
"if" asks
whether the middle point and the right point have different signs,and
if so to
set the new interval to be the second half of the old interval.
> for k from 1 to n do;
>    c.(k-1):=(a.(k-1) + b.(k-1) ) / 2;
>    if( f(c.(k-1))*f(b.(k-1)) <= 0 ) then
>       a.k:= c.(k-1);  b.k:= b.(k-1);
>    else
>       a.k:=a.(k-1) ; b.k:= c.(k-1);
>    fi;
> od:
> c.n := (a.n + b.n ) / 2:
> seq([a.k,b.k,c.k,b.k-c.k],k=0..n);
EXTRA: We provide a graphic version: It shows a sequence of shorter
and shorter nested intervals which bracket the zero. This time
plot(T,x=a.0..b.0); where  T is a set of lists of points. ([1,4] is a
single point,
[[1,4],[3,0],[5,4]] is a list of three points and
{[[1,4],[3,0],[5,4]]} is a set consisting of a list of points.
For example if T:={[[1,4],[3,0],[5,4]]} then"plot(T,x=a.0..b.0);"
produces on "vee" shaped zigzag. Then "T:=T union
{[[2,4],[3,2],[3,4]]};" will add
a second "vee" to the set  T  and the print instruction
"plot(T,x=a.0..b.0);"
will draw two vees.
> es:=0.2*max(abs(f(a.0)),abs(f(b.0))):
> ep:= es/n:
> T:={[[a.0,es],[a.0,-es],[b.0,-es],[b.0,es]]}:
> for k from 1 to n do;
>    ez := k * ep:
>    T:=T union {[[a.k,es-ez],[a.k,ez-es],[b.k,ez-es],[b.k,es-ez]]}:
> od:
> plots[display]({
>    plots[polygonplot](T),
>    plot(f(t),t=a.0..b.0),
>    plot([[c.n,2*es],[c.n,-2*es]],x=a.0..b.0,style=line),
>    plots[textplot]([c.n,2.5*es,`Root`])
> }, title=`Bisection Method`);
METHOD OF ITERATION - - SPIDER WEBS

The idea to solve the equation  x=f(x) where  f(x)  is a continuaous
function,  is to guess some solution  x.0  and try to improve it by
replacing it by  x.1 = f(x.0), then x.2= f(x.1)  and so on. In general
we let x.n=f(x.(n-1)).  This procedure will give a solution if we can
show
that  the sequence  converges to some number 
z=limit(x.n,n=1..infinity); If this is the case, then  z =
limit(x.n,n=1..infinity) =
limit(f(x.n),n=1..infinity) =f(limit(x.n,n=1..infinity))=f(z).  The
theory says that the sequence
will converge provided that we can find two numbers a < b  so that
|f'(x)|<= c < 1 for all a <= x <= b  and that  if  x  is any number in
the interval a <= x <=b  then  a <= f(x) <= b.
Let's try it in a particular case. Suppose we wanted to find the
square root of  two using iteration. That is, suppose we wanted to
solve the equation  0 = 2-x^2.  We have to convert this to  the form 
x=f(x). 
To do this we may try adding  x  to both sides.  Then  x = 2 - x^2 +
x. Now if 
f(x) = 2 - x^2 + x then f' = -2*x + 1 and |f'(x)|<1 whenever 0<x<1 
but  we
expect the root to be greater than  1.  To remedy this situation we
try instread
"f(x)=(2-x^2)/4+4;".  This time again the square root satisfies 
"z=f(z)" but  |f'(x)|<1 whenever  0<x<4  surrounding the desired root.
 Since the maximum occurs for  x=1 when f(1)=2, we can take
 a=0.1 and  b = 2  (which is OK since we expect the root between  0.1
and 2). Then
|f'(x)|<=.95 for every  x satisfying  0.1 <= x <= 2.  
Also  if  0.1 <=x <= 2
then  0.1 <= f(x) <= 2 since f is increasing in [0.1,2] and f(0.1) =
.59975>0.1.  Lets try a few
iterations. We starty with  x.0=0.3 but any guess in [a,b] will work.
 
> f:=t->(2-t^2)/4+t;
> x.0:=0.3;`Actual square root of two =`,sqrt(2.0);
> for k from 1 to 10 do;x.k:=f(x.(k-1));od;
Now we proceed to display the iterations more graphically. We
reinitialize the iteration. We will start from the point[x.0,0]. Then
draw a
vertical line until we reach [x.0,x.1]  where  x.1=f(x.0). Then make a
horizontal
line to  [x.1,x.1], in effect replascing y by  x. Then we make a
vertical line
to  [x.1,x.2]  where x.2=f(x.1). And proceed throuch several
iterations.
The  instruction  "plot(S,x=0..2,style=point;" plots a zig-zag line
through
the points in the list S. For example, if 
"S:=[[1,0],[1,2],[2,2],[2,1]];"
then the plot will be line segments from  (01,0) to (1,2) to (2,2) to
(2,1).
 The way that you can add points to the list. For example if you
execute "S:=[op(S),[1,1],[1,0]];" then  "S"  will now contain
[[1,0],[1,2],[2,2],[2,1],[1,1],[1,0]]; By ending the  "do"  loop with
colon means it doesn't print.  After accumulating the points "S" we
plot
them together with a plot of "y=x" and "y=f(x)";   Try running this
loop
with different initial x.0's.
> x.0:=0.1:S:=[[x.0,0]]:
> for k from 1 to 6 do;
> x.k:=f(x.(k-1));
> S:=[op(S),[x.(k-1),x.k],[x.k,x.k]];
> od:
> 
> plots[display]({
> plot({x,f(x)},x=0..3,scaling=constrained),
> plot(S,x=0..4,style=line,color=blue)
> },title=`Iteration Method`);
We get a different picture if  F  is decreasing. This time lets do an 
viteration to find square root of  3.
> `Actual root of  3  is  `,sqrt(3.0);
Consider ther function  "g(x)=(3-x^2)/2+x;"  Again the root satisfies
"z=g(z)". Now  g'(x)=-x+1 so that  |g'(x)|<1 if  0 < x < 2. Now if  z 
is the square root of  3, then 1.2 < z < 1.98  because  1.2^2=1.44 < 3
<
1.98^2 = 3.9204. On the other hand  g(1.2)=1.98 and g(1.98)=1.5189.
Since  g 
is decreasing in the interval [a,b]  where  a=1.2 and b=1.98 we get
that 
if 1.2 <= x <= 1.98 then  1.2 <= g(x) <= 1.98. Furthermore, if 1.2 <=
x <=
1.98 then
|g'(x)| <= .98 < 1.  Thus again the conditions for the existence of a
limit under iteration is guaranteed.  (In practice one doesn't have to
check
the conditions for  g  but one risks the possibility of chaos. After
generating the points, we plot list out all the  x's.
> g:=t->(3-t^2)/2+t;
> x.0:=0.05:S:=[[x.0,0]]:
> for k from 1 to 15 do;x.k:=g(x.(k-1));
> S:=[op(S),[x.(k-1),x.k],[x.k,x.k]];
> od:
> plots[display]({plot({x,g(x)},x=0..3,scaling=constrained),
> plot(S,x=0..4,style=line,color=blue)
> },title=`Iteration method`);
Now lets see what happens if the derivative doesn't stay bounded. We
may get chaotic behavior. For example,  if  q(x)= 4*x-2*x^2. Then  if 
 0 <= x <= 2 
we have 0 <= q(x) <= 2.Try various initial data. 
> q:=t-> 4*t-2*t^2;
> x.0:=0.1:S:=[[x.0,0]]:
> for k from 1 to 100 do;x.k:=q(x.(k-1));
> S:=[op(S),[x.(k-1),x.k],[x.k,x.k]];
> od:
> plots[display]({plot({x,q(x)},x=0..2,scaling=constrained),
> plot(S,x=0..2,style=line,color=gold)
> },title=`Chaos`);

NEWTON'S METHOD
The idea for Newton's method is to find the tangent line at each guess
point [x.k,f(x.k)]. then select the new point to be the intercept with
the 
x-axis. The iteration follows "x.k:=x.(k-1)- f(x.(k-1))/f'(x.(k-1));"
This
method converges fantastically fast! For example, to compute the
square root
of seven, we store the derivative in "fp". In MAPLE  if  "f"  is a
function then "D(f)"  is the derivative function.  The instruction
"diff(expr,x);"
differentiates expressions and leaves expressons.
> f:=t-> 7-t^2;
> fp:=D(f);
> x.0:=5.111111111:sre:=sqrt(7.0);
> for m from 1 to 6 do;
>    x.m:= x.(m-1) - f(x.(m-1))/fp(x.(m-1));
>    printf(`   % 16.10f   % 16.10f \n`,x.m,x.m-sre);
> od:
> 
Try to find the cube root of 13. This time we plot the function and
the
tangent lines.
> g:=x->x^3-13;
> gd:=D(g);
> x.0:=3.8999999998;
> cr:=(13.0)^(1/3);
> S:=[[x.0,0],[x.0,g(x.0)]]:
> for m from 1 to 6 do;
>    x.m:=x.(m-1) - g(x.(m-1))/gd(x.(m-1));
>    S:=[op(S),[x.m,0],[x.m,g(x.m)]];
>    printf(`   % 16.10f   % 16.10f \n`,x.m, x.m-cr);
> od:
EXTRA. We plot Newton's method in action. The diagonals ate tangent
lines. 
Then we zoom in.
> plots[display]({
>    plot(g(x),x=0..4,color=blue),
>    plot(S,x=0..4,style=line)
> },title=`Newton's Method`);
> 
> plots[display]({
>    plot(g(x),x=2.3..2.5,y=-1..2,color=blue),
>    plot(S,x=2.3..2.5,style=line)
> },title=`Newton's Method`);


ASSIGNED PROBLEMS FOR LAB 5.

1.  According to the Greeks, the most beautiful rectangle has the
proportions  x : 1  where  x>1  is the "golden mean" (see 520[52]). If
one removes a square from a golden rectangle, then the remaining
subrectangle is similar to the original rectangle. Find an equation
for  x. (1=x^2-x)  Solve the equation using  algebra. Solve the equation to four
digits
accuracy(error less than 5 x 10^(-5)) using Bisection. Estimate the
error.
> plot([[1,0],[1,1],[0,1],[0,0],[1.618,0],[1.618,1],[1,1]],x=-.2..2,scal
> ing=constrained, axes=none, title=`Golden Rectangle`);

2.  Find the golden mean using Newton's method.  By making one or two
test runs, show that it may take more iterations for poor choices of
the
initial guess  x.0.
3.  Find the golden mean using the method of iteration.  You may have
to tweak your iteration as we did above to get your procedure to
converge. Make a "spiderweb" picture showing your convergence. 
4.  Here is a question from statistics.  Using any one of the three
methods, find a number  z  so that  the function 
erf(z/sqrt(2.0))/2=.475. 
This means that  95% of the time, a standard normal variable will lie
between  -z
and z. This  z can be used to give 95% confidence intervals for how
well the
mean of a sample reflects the mean of a normally distributed
population in
case the population standard deviation is known.