# This worksheet contains the Maple commands from Chapter 9 of
# Introduction to Scientific Programming by Joseph L. Zachary.
# 
# (9.1)  The area of the "slice" of circle from Figure 9.1.
# 
> X := theta/(2*Pi) * Pi*R^2;
# 
# (9.2)  The length of the side R from Figures 9.2 and 9.3.
# 
> R := 2 * r * cos(theta/2);
# 
# (9.3)  The area of each of the triangles from Figure 9.2.
# 
> Y := 1/2 * R * r * sin(theta/2);
# 
# (9.4)  The area of region Z from Figure 9.2.
# 
> Z := (2*theta)/(2*Pi) * Pi*r^2;
# 
# (9.5)  The area of the region where the cow will be unable to graze.
# 
> A := Z - (X - 2*Y);
# 
# (9.6)  We simplify the expression from (9.5) by using the combine
# function.
# 
> combine(A, trig);
# 
# (9.9)  A simple example of using the solve function.
# 
> solve(2*x = 4, x);
# 
# (9.10)  Using the solve function to solve two simultaneous linear
# equations.
# 
> solve({3*x + 2*y = 7, x + y = 3}, {x,y});
# 
# (9.11)  Using the solve function to solve a quadratic equation.
# 
> solve(2*x^2 - 3*x - 9 = 0, x);
# 
# (9.12)  Using the solve function to solve a trigonometric equation.
# 
> solve(cos(alpha) = sin(alpha), alpha);
# 
# (9.13)  Using the solve function to solve a symbolic quadratic
# equation.  Notice that the result is the quadratic formula.
# 
> solve(a*x^2 + b*x + c = 0, x);
# 
# (9.14)  Unfortunately, solve can't solve just anything.  Here we try
# (and fail) to solve the equation that is at the root of Old
# MacDonald's problem.
# 
> solve(sin(theta) - theta*cos(theta) = Pi/2, theta);
# 
# (9.15)  The fsolve function uses numerical techniques to solve
# equation.  It often succeeds where solve fails.
# 
> angle := fsolve(sin(theta) - theta*cos(theta) = Pi/2, theta);
# 
# (9.16)  We compare the values of the two sides of the equation from
# (9.15) when the solution (angle) is substituted for theta.
# 
> sin(angle) - angle*cos(angle);
> evalf(Pi/2);
# 
# (9.17)  Because it operates numerically, fsolve cannot be used to find
# symbolic solutions.
# 
> fsolve(a*x^2 + b*x + c = 0, x);
# 
# (9.18)  We need to find the root of this function to solve Old
# MacDonald's problem.
# 
> cow := (theta) -> sin(theta) - theta*cos(theta) - evalf(Pi/2);
# 
# (9.19)  This demonstrates that the cow function is positive at 2.5 and
# is negative at 1.5.  Thus, a root of cow must lie somewhere in
# between.
# 
> cow(2.5);
> cow(1.5);
# 
# (9.20)  This plot verifies that a root of cow lies between 1.5 and
# 2.5.
# 
> plot(cow(theta), theta = 1.5 .. 2.5);
# 
# (9.21)  Midway between 1.5 and 2.5, cow is positive.  This means that
# a root lies somewhere between 1.5 (where cow is negative) and 2.0.
# 
> cow(2.0);
# 
# (9.22)  This plot verifies that a root of cow lies between 1.5 and
# 2.0.
# 
> plot(cow(theta), theta = 1.5 .. 2.0);
# 
# (9.23)  Midway between 1.5 and 2.0, cow is negative.  This narrows the
# interval in which a root of cow must lie to 1.75--2.0.
# 
> cow(1.75);
# 
# (9.24)  It would be nice if we had a bisection procedure that worked
# as illustrated below.  We would pass it the function whose root we
# wished to find, a value at which the function was positive, and a
# value at which the function was negative.  For the remainder of this
# worksheet we will work towards developing such a procedure.
# 
> bisection(cow, 2.5, 1.5);
# 
# (9.25)  The order in which statements are executed out can make a
# difference.  Compare the final values of "x" and "y" in this sequence
# ...
# 
> x := 1:
> x := 2;
> y := x^2;
# 
# (9.26) ... to the final values of "x" and "y" in this sequence.  The
# only difference is that we have permuted the order of the last two
# assignment statements.
# 
> x := 1:
> y := x^2;
> x := 2;
# 
# (Figure 9.5)  This is the beginning of our attempt to automate the
# bisection process.  We begin with a function f to be solved and values
# pos and neg such that f(pos)>0 and f(neg)<0.  We repeatedly average
# pos and neg, test the value of f at the average, and update either pos
# or neg appropriately.  Notice that we (as the creators of the
# statement sequence) had to determine whether to update pos or neg at
# each step.  If we were to repeat this process enough times, we would
# narrow down on a root of f.
# 
> f := cow;
> pos := 2.5;
> neg := 1.5;
> 
> avg := (pos + neg) / 2.0;
> f(avg);
> pos := avg;
> 
> avg := (pos + neg) / 2.0;
> f(avg);
> neg := avg;
# 
# (Figure 9.6)  Here we have inserted conditional statements into our
# statement sequences.  These conditionals decide whether to update pos
# or neg at each step. 
# 
> f := cow;
> pos := 2.5;
> neg := 1.5;
> 
> avg := (pos + neg) / 2.0;
> if (f(avg) >= 0) then
>   pos := avg;\
> else
>   neg := avg;
> fi;
> 
> avg := (pos + neg) / 2.0;
> if (f(avg) >= 0) then
>   pos := avg;
> else
>   neg := avg;
> fi;
# 
# (Figure 9.7)  The example above contained two consecutive identical
# statement blocks.  Here, we make that statement block the body of a
# while statement, which repeatedly executes its body until the
# difference between pos and neg is small.
# 
> f := cow;
> pos := 2.5;
> neg := 1.5;
> 
> while (abs(pos-neg) >= 1e-6) do
>   avg := (pos + neg) / 2.0;
>   if (f(avg) >= 0) then
>     pos := avg;
>   else
>     neg := avg;
>   fi;
> od;
# 
# (9.27)  After the while loop above has run, the value of avg is fairly
# close to being a root of f.
# 
> f(avg);
# 
# (Figure 9.8)  This procedure encapsulates the bisection process.  It
# takes three parameters:  the function to be solved, a positive guess
# such that function(positive)>0, and a negative guess such that
# function(negative)<0.  It returns an approximation to a root of
# function.
# 
> bisection := proc (function, positive, negative)
> 
>   local f, avg, pos, neg;
> 
>   f := function;
>   pos := positive;
>   neg := negative;
>   while (abs(pos-neg) >= 1e-6) do
>     avg := (pos + neg) / 2.0;
>     if (f(avg) >= 0) then
>       pos := avg;
>     else
>       neg := avg;
>     fi;
>   od;
> 
>   RETURN(avg);
> 
> end;
# 
# (9.28)  We use the bisection procedure to approximate a root of cow.
# 
> bisection(cow, 2.5, 1.5);
# 
# (9.29)  We compare the roots to cow found by bisection and fsolve. 
# They differ in the last three digits of their mantissas, which is
# evidence that the convergence criterion used by bisection could be
# improved.
# 
> bisection(cow, 2.5, 1.5);
> fsolve(cow(theta), theta);
# 
# (9.30)  If we try to use bisection to solve this simple linear
# equation, the while loop that it contains will never terminate because
# the positive and negative guesses will never get very close together. 
# If you evaluate this expression, you will need to interrupt the
# execution.
# 
> bisection((x) -> x - 20000., 30000., 10000.);
# 
# (9.31)  If we try to use bisection to solve this linear equation, we
# won't get very close to a root (on a relative basis).  The root is so
# small that the positive and negative guesses get within 10^(-6) very
# quickly.
# 
> bisection((x) -> x - 1e-12, 1.0, 0);
# 
# (9.32)  In this equation, the two guesses are so close together to
# begin with that the while loop within bisection is never executed,
# which means that avg is never given a value, which means that this
# useless symbolic result is reported.
# 
> bisection((x) -> x - 1e-12, 2e-8, 0);
# 
# (9.33)  If we supply faulty positive and negative guesses, bisection
# will converge to one of the guesses, and not to a root.
# 
> bisection((x) -> cos(x) - x, 1.0, 0.5);
# 
# (Figure 9.9)  This is an improved version of bisection.  It makes sure
# that positive and negative have the correct properties, make sure that
# a number is returned as the result under all circumstances, and uses
# an improved convergence criterion.
# 
> bisection := proc (function, positive, negative)
> 
>   local f, avg, pos, neg;
> 
>   f := function;
>   pos := positive;
>   neg := negative;
> 
>   if (f(pos) < 0) then
>     print(`Positive bound to bisection is bad`);
>     RETURN();
>   fi;
> 
>   if (f(neg) > 0) then
>     print(`Negative bound to bisection is bad`);
>     RETURN();
>   fi;
> 
>   avg := pos;
> 
>   while (avg <> (pos+neg)/2.0) do
>     avg := (pos + neg) / 2.0;
>     if (f(avg) >= 0) then
>       pos := avg;
>     else
>       neg := avg;
>     fi;
>   od;
> 
>   RETURN(avg);
> 
> end;
# 
# (9.34)  We compare the improved bisection to fsolve.  Now they differ
# only in the last digit of the mantissa.
# 
> bisection(cow, 2.5, 1.5);
> fsolve(cow(theta), theta);
>