Chapter 3: Worksheet 10 Jack K. Cohen Colorado School of Mines




Newton's Method, Equation Solving with


Suggested Problems Section 3.9: Read section 3.9.




  1. Illustrative Example 2 (the ``cork ball problem'') gives the formula for the volume of a segment of height x and radius r of a sphere whose radius is a as V = ${\frac{{\pi x}}{{6}}}$ (3r2 + x2) (see Figure 1). The same volume was cited on the last worksheet as V = ${\frac{{\pi x^2}}{{3}}}$ (3a - x).
    Figure: Sketch for Problem 1
    \begin{figure}
\epsfysize 110pt
\centerline{\epsffile{ws10p1.eps}}
\end{figure}

    1. Show that these formulas are equivalent. (note: The quantity here denoted by x is also denoted by y and h in other problems and examples in the text.)
    2. Suppose that the density of water is 1 and the density of the ball is ρ. Appeal to Archimede's law to obtain the equation ρVsphere = Vsegment. Using the second form above for the volume of the segment, obtain the cubic equation 4ρa2 = 3ax2 - x3. Show that when, as in the illustrative example, ρ = 1/4 and a = 1, we obtain the equation 3x2 - x3 = 1 in agreement with the text.
    3. Use the Intermediate Value Theorem to prove that the explicit cubic equation, x3 -3x2 + 1 = 0, just derived has at least one root in the interval [0, 1].

  2. One way to solve equations with is to code in your own Newton's method:
    f[x_] := x^3 - 3x^2 + 1
    x = 0.5
    
    x = x - f[x]/f'[x]
    
    Here we have defined the ``cork ball'' function in accordance with syntax and defined a starting value at the midpoint of an interval that we have proved contains a root. The third line computes the ``new'' value of x (xn+1) on the left from the ``old'' value (xn) used on the right (the ``old'' value is 0.5 the first time we do this). We can then ``Enter'' the third line over and over again to generate a sequence of approximations. On this example, you should get successively:
    0.666667
    0.652778
    0.652704
    0.652704
    0.652704
    
    This agrees with the value found in the text. Apply this approach to the following problems:
    1. (3.9.1) x2 -5 = 0;    x0 = 2.5 (i.e. find the positive square root of 5)
    2. (3.9.2) x3 -2 = 0;    x0 = 1.5 (i.e. find the cube root of 2)

  3. One problem with the above method is that we erase each approximation in the process of producing the next one. One simple way to avoid this is to use a ``loop'':
    f[x_] := x^3 - 3x^2 + 1
    x = 0.5;
    Do[
        x = x - f[x]/f'[x];
        Print[x],
    {8}]
    
    which gives the ``cork ball'' output:
    0.666667
    0.652778
    0.652704
    0.652704
    0.652704
    0.652704
    0.652704
    0.652704
    
    The ``8'' that appears in curly brackets at the end of ``Do Loop'' tells how many times to iterate.

    Important and subtle point: By default, tries to compute with exact arithmetic. This can take a long time, so it is usually best to include a decimal point (e.g., 0.5 instead of 1/2) in the starting value to signal to that infinite precision is not required.

    Since Newton's method often converges quickly, let's print the results to greater precision by altering the Print statement:

    f[x_] := x^3 - 3x^2 + 1
    x = 0.5;
    Do[
    	x = x - f[x]/f'[x];
    	Print[ N[x, 16] ],
    {8}]
    
    0.6666666666666666
    0.6527777777777778
    0.6527036468361321
    0.6527036446661393
    0.6527036446661393
    0.6527036446661393
    0.6527036446661393
    0.6527036446661393
    

    Apply this approach to the following problems:

    1. (3.9.3) x5 -100 = 0;    x0 = 5.0 (i.e. find the fifth root of 100)
    2. (3.9.4) x3/2 -10 = 0;    x0 = 2.0 (i.e. compute 102/3)
    3. x1/3 -1;    x0 = 2.0 (i.e. compute 13 the hard way)
    4. x1/3 -1;    x0 = 5.0
    5. Explain why the output in the last two examples differ.

  4. There are several equation solvers directly built in to —be aware that Newton's method is a key element in their construction! However, other considerations are also needed in such professional packages to deal with potential ``divergent'' behavior as occured in the last example. Polynomial equations are markedly simpler than the general case and Solve and NSolve deal with these and with some equations that can be reduced to polynomials by algebraic manipulations (e.g. by squaring). Solve attempts to give exact solutions, while NSolve gives only numerical solutions. As implied, Solve cannot exactly solve all polynomial equations (there is no mathematical method to do this for equations of degree 5 and higher). However, Solve can sometimes handle equations with parameters.

    First a cubic with a parameter—Solve is our only chance:

    equation = x^3 - 2a x^2 - 4 x + 8a == 0;
    solutions = Solve[equation, x]
        {{x -> 2 a}, {x -> 2}, {x -> -2}} (* the solutions written as rules *)
    
    equation /. solutions     (* check solutions by substituting back *)
       {True, True, True}    (* all 3 check exactly *)
    
    Solve works ok here, because the equation was special (an answer that is 17 pages long isn't of much use and that can easily happen with an arbitrary cubic). Here is a quadratic for which we get exact solutions with Solve:
    equation = x^2 - 3x + 4 == 0;
    Solve[equation, x]
               3 + I Sqrt[7]         3 - I Sqrt[7]
        {{x -> -------------}, {x -> -------------}}
                   2                     2
    
    If we only want numerical approximations, it is more efficient to use NSolve:
    equation = x^2 - 3x + 4 == 0;
    NSolve[equation, x]
        {{x -> 1.5 - 1.3228756555323 I}, 
          {x -> 1.5 + 1.3228756555323 I}}
    

    Use NSolve to get all three roots of the ``cork ball'' equation.

    Figure: Diagram for Problem 5.
    \begin{figure}
\epsfysize 100pt
\centerline{\epsffile{ws10p5.eps}}
\end{figure}

  5. A 100-ft tree stands 20 feet from a 10-ft fence. Then the tree is ``broken'' at a height of x feet, so that it just grazes the top of the fence and touches the ground on the other side of the fence with its tip (see Figure 5).
    1. Derive the equation x3 -68x2 + 1100x - 5000 = 0 for x.
    2. Solve this equation using NSolve and decide which roots make physical sense.

  6. NSolve won't work on most non-polynomial equations, for example, it fails on our old friend the ``asteroid'' equation:
    equation = Cos[1000/r] == r/(r + 10);
    NSolve[equation, r]
        Solve::ifun: 
           Warning: Inverse functions are being used by Solve, so
             some solutions may not be found.
         Solve::tdep: 
           The equations appear to involve transcendental
             functions of the variables in an essentially
             non-algebraic way.
                  1000       r
        NSolve[Cos[----] == ------, r]
                    r       10 + r
    
    All those messages are NSolve's way of saying it can't solve the equation. In this case, you must use FindRoot. FindRoot requires additional information—at the minimum a starting guess (as does Newton's method). By previous trial and error, we know there's a root near 50,000, so:
    lhs = Cos[1000/r] - r/(r + 10);
    FindRoot[lhs == 0, {r, 50000}]
        {r -> 50008.3}
    
    solution = N[ FindRoot[lhs, {r, 50000}], 16]
        {r -> 50008.32913584035}
        
    lhs /. solution
                   -11
        -1.10986 10
    
    Notice the techniques for printing more decimal places and for checking the accuracy of the solution. Also observe that it is still an open question of how to get that starting value without trial and error—the next Chapter will give us a technique for doing that.

    1. Plot tan x - x on [0, 4π] to get crude values for the positive roots near the origin.
    2. Use FindRoot to get accurate values for the first 4 positive roots of the equation tan x + x = 0 and check the accuracy of your solutions.