This notebook contains the Mathematica 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 "Simplify" function.

Simplify[A]

(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 "FindRoot" function uses numerical techniques to solve equation.  It often succeeds where "Solve" fails.

angle =
  FindRoot[Sin[theta] - theta*Cos[theta] == Pi/2,
           {theta, 1.0}]

(9.16)  We compare the values of the two sides of the equation from (9.15) when the solution ("angle") is substituted for "theta".  This is the left-hand side ...

ans1 =  ReplaceAll[Sin[theta] - theta*Cos[theta],
                   angle]

(9.17)  ... and this is the right-hand side.

ans2 = N[Pi/2]

(9.18)  The two values look to be the same.  But if we look in more detail at the left-

InputForm[ans1]

(9.19)  and right-hand sides, we see that they are quite different in the low-order digits.

InputForm[ans2]

(9.20)  We can use the "AccuracyGoal" optional parameter to ask "FindRoot" to deliver a more accurate answer.

newangle =
  FindRoot[Sin[theta] - theta*Cos[theta] == Pi/2,
           {theta, 1.0},
           AccuracyGoal->16]

(9.21)  Now we can determine the value of the right-hand side ...

ans3 = ReplaceAll[Sin[theta] - theta*Cos[theta],
                  newangle]

(9.22)  ... and see that it is much improved.

InputForm[ans3]

(9.23)  Because it operates numerically, "FindRoot" cannot be used to find symbolic solutions.

FindRoot[a*x^2 + b*x + c == 0, {x, 1.0}]

(9.24)  We need to find the root of this function to solve Old MacDonald's problem.

cow[theta_] := Sin[theta] - theta*Cos[theta] - Pi/2

(9.25)  This demonstrates that the "cow" function is positive at 2.5 ...

cow[2.5]

(9.26)  ... and negative at 1.5.  Thus, a root of "cow" must lie somewhere in between.

cow[1.5]

(9.27)  This plot verifies that a root of "cow" lies between 1.5 and 2.5.

Plot[cow[theta], {theta, 1.5, 2.5}]

(9.28)  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.29)  This plot verifies that a root of "cow" lies between 1.5 and 2.0.

Plot[cow[theta], {theta, 1.5, 2.0}]

(9.30)  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.31)  It would be nice if we had a bisection function 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 notebook we will work towards developing such a function.

bisection[cow, 2.5, 1.5]

(9.32)  This is an example of an expression sequence.  The expressions are evaluated from left to right, and only the value of the last one is displayed.  Since all three of these expressions are pure expressions, there is no point to arranging them in a sequence.

(1+1; 2+2; 3+3)

(9.33)   The order in which imperative expressions are evaluated can make a difference.  Compare the final value of "y" in this sequence ...

(x=1; x=2; y=x^2; y)

(9.34)  ... and in this sequence.  The only difference is that we have permuted the order of the second and third assignments.

(x=1; y=x^2; x=2; y)

(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 expression 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;
 
 avg)

(Figure 9.6)  Here we have inserted conditional expressions into our expression sequence.  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,
    pos = avg,
    neg = avg];
    
 avg = (pos + neg) / 2.0;
 If[f[avg] >= 0,
    pos = avg,
    neg = avg];
    
 avg)

(Figure 9.7)  The example above contained two consecutive identical expression sequences.  Here, we make that expression sequence the body of a "While" expression, 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] >= 1*^-6,
       (avg = (pos + neg) / 2.0;  
        If[f[avg] >= 0,
           pos = avg,
           neg = avg])];
           
 avg)

(9.38)  After the while loop above has run, the value of "avg" is fairly close to being a root of "cow".

cow[avg]

(Figure 9.8)  We embed the expression from (Figure 9.7) inside a module.  The variables "f", "pos", "neg", and "avg" are now local to the module.

Module[{f, pos, neg, avg},
       (f = cow;
        pos = 2.5;
        neg = 1.5;
        
        While[Abs[pos-neg] >= 1*^-6,
              (avg = (pos + neg) / 2.0;  
               If[f[avg] >= 0,
                  pos = avg,
                  neg = avg])];
                  
        avg)]

(Figure 9.9)    This function 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[function_, positive_, negative_] :=
  Module[{f, pos, neg, avg},
         (f = function;
          pos = positive;
          neg = negative;
          While[Abs[pos-neg] >= 1*^-6,
                (avg = (pos + neg) / 2.0;  
                 If[f[avg] >= 0,
                    pos = avg,
                    neg = avg])];
          avg)]

(9.40)  We use the "bisection" procedure to approximate a root of "cow".

bisection[cow, 2.5, 1.5]

(9.42)  We compare the roots to "cow" found by "FindRoot" ...

InputForm[FindRoot[cow[theta], {theta, 1.0},
          AccuracyGoal->16]]

(9.43)  ... and "bisection".  They differ markedly, which is evidence that the convergence criterion used by "bisection" could be improved.

InputForm[bisection[cow, 2.5, 1.5]]

(9.44)  This function has a single large root.

bigRoot[x_] := x - 2.*^20

(9.45)  If we try to use "bisection" to find the function's root, 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[bigRoot, 3.*^20, 1.*^20]

(9.46)  This function has a single small root.

smallRoot[x_] := x - 1.*^-20

(9.47)  If we try to use "bisection" to find the function's root, 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[smallRoot, 1.0, 0.0]

(9.48)  In this case, 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 a useless symbolic result is reported.

bisection[smallRoot, 2.*^-8, 0]

(9.49)  If we supply faulty positive and negative guesses, "bisection" will converge to one of the guesses, and not to a root.

bisection[bigRoot, 1.*^20, 3.*^20]

(Figure 9.10)  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[function_, positive_, negative_] :=

  Module[{f, pos, neg, avg, delta},
  
         (f = function;
          pos = positive;
          neg = negative;
          
          If[f[pos]<0, Return["Positive guess is bad"]];
          If[f[neg]>0, Return["Negative guess is bad"]];
          
	      delta = Abs[pos-neg] * 2;
	      
          While[Abs[pos-neg] < delta,
                (avg = (pos + neg) / 2.0;
                 delta = Abs[pos-neg];
                 If[f[avg] >= 0,
                    pos = avg,
                    neg = avg])];
                    
          avg)]

(9.51)  We can now compare the improved "bisection" to "FindRoot".

InputForm[bisection[cow, 2.5, 1.5]]