Newton's Amazing Method: Estimate π to How Many Places?

Introduction

OBJECTIVE:  Use Mathematica to implement Newton's Method.

To how many places can you estimate π or e or ? In this module, Isaac Newton and Mathematica team up to help you in this endeavor. Estimates to thousands of places are at your fingertips when you have such an awesome pair to help you. You will also encounter some of the hazards in using Newton's Method.

Technology Guidelines

NOTE:  If you have just finished a module, restart Mathematica or close the Kernel before executing a new module.
TO OPEN CELLS, put your cursor on the right cell bracket and double click.

INITIALIZATION CELLS
When asked if you want to ". . . automatically evaluate all the initialization cells in the     notebook . . . ," respond by pressing the "Yes" button.
TO STOP AN EXECUTION
Select the Kernel pull-down menu and click on Abort Evaluation.
ORDER OF EXECUTION
Execute cells in the order given. Do not skip any Input cells within a given notebook.
SAVING NOTEBOOKS
You can save anytime to any directory you choose, and it is wise to save often.
However, before you do your final save, delete all your output by selecting the
Delete All Output selection under the Kernel pull-down menu.
EXPERIENCING MAJOR PROBLEMS
Save if appropriate, then shut down Mathematica and start it up again.

Part I: Using Newton's Method to Estimate π

Chapter 4, Section 7, Exercise 20

Setting Up the Problem

For those of us who enjoy mathematical challenges, estimating irrational numbers like π, e, and to a large number - no, to a very large number of decimal places - is an exercise of enjoyment. Let's see if we can crash the computer!

To estimate π, we will use Newton's method to estimate the zero of the function f(x) = tan x , for π/2 < x < 3π/2. As you know already, the exact answer is π. First, we define the function f[x_ ]. Clear the symbol x so that there is no conflict in defining the Newton iteration function. Later on you will want to change the function f[x_ ] and clearing x now will avoid complications when you redefine the iteration function for the new f[x_ ].

In[1]:=

">

Now we plot f(x) over the specified domain of interest to determine the approximate location of an x-value that makes the function zero.

In[3]:=

The plot window is cut down by specifying PlotRange → {-2, 2} so that we can see the root clearly.

Applying Newton's Method

We set up Newton's iteration scheme thus.

In[4]:=

We have set the precision at 100 digits for the estimate of the zero. After we go through this once, you will want to come back and change this value, but leave it at 100 for now.

">

Using the graph of the function, we specify a starting value of x=3 for the iterations. We iterate until two consecutive calculated values of x are equal to the precision specified by digits. The symbol n is used to count the number of iterations required to achieve the specified precision. For an explanation of the following commands, click on the "About Mathematica" button below the cell.

In[7]:=

Out[16]//TableForm=

 n x 0 3 1 3.1397077490994629364057777233059473798139974321591021592416756848266555770293221674273344692264211238115048681895887`100. 2 3.1415926491252556944794381440657649099212399776512705648543836226501752887271487430517960635836043249274480038319715`100. 3 3.1415926535897932384626433239544438121837693370155904765056022068985170791187879984488515625960037271773091241873907`100. 4 3.1415926535897932384626433832795028841971693993751058209749445923078164062860698037422630016306266327627447916368344`100. 5 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823`100. 6 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480639842695`100.

">

In the cell above, we first initialize the iteration counter, n, to 0. We want to store each iteration count and the corresponding value of x in a list called xiterates, so we initialize the list with the first pair of values, n (=0) and x (=initial guess). Then we assign the current value of x to the symbol y for later comparison with the next computed value of x. We use the newt[ ] function to calculate the next value of x and increment the iteration count n by 1 with n++ (at least one iteration will always be performed). Now we enter the While[ ] command. Inside the While[ ] command, the values of x and y are tested to see if they are equal. If they are not, the other commands are executed in order: the current value of x is again assigned to y for the next test; a new value of x is computed; n++ increments the iteration counter n by 1; and the new values for n and x are appended to the list of xiterates. This procedure is repeated until the condition is not satisfied, that is, until two consecutive estimates of x are equal to the precision specified in digits.

From the previous results, we estimate π to 100 decimal places.

In[18]:=

Out[18]=

Testing the Results

How do we know that Newton's method is working correctly and that our estimate of π is correct? We can compare it with Mathematica's estimate of π to the same number of digits.

In[19]:=

Out[19]=

If Mathematica's estimate of π is correct to the specified number in digits, then so is ours. Not bad!

In[20]:=

Out[20]=

What is truly amazing is the rate at which Newton's method converges to the root. We can achieve 100 digits of accuracy with so few iterations - Wow! And because it takes so few iterations, it is also very fast. But wait, it gets even better. Go back up and change digits to 1000 or even 10000. See what happens. (Be aware, however, that at some point you will run up against the limitations of your computer, specifically, RAM and cache storage on your hard drive, and it may crash.) Also, when you increase the precision, you have two choices for the initial value of x. If you don't reset it with x=initial value, the iteration will start with the last value calculated (i.e., x to 100 digits). If you reset it, then the first 100 iterations will be repeated on the way to 1000, 10000, 1000000, . . . .

You Try It: Part I

You may also want to do some research to find out what the current best estimate of π is. Research the history of the efforts that mathematicians have put into this problem.

See if you can construct functions so that you can use the above code to estimate .  Begin by defining a function that is zero when  x=. We suggest that you try a polynomial. Replace the terms in red with the appropriate function, and change the bounds over which you plot f(x), if you wish.

In[21]:=

Now set up Newton's iteration scheme. You can set digits to whatever precision you wish.

In[22]:=

Using the graph of your function, specify a starting value of x (in red) for the iterations. We iterate until two consecutive calculated values of x are equal to the precision specified by digits, and we use the symbol n to count the number of iterations required to achieve the specified precision.

In[23]:=

Out[23]=

The output from the code above renders your estimate as the following. Compare it to the , and see if it has given the desired precision. It clearly does not if you did not change the function f[x_] (in red).

In[24]:=

Out[24]=

Out[25]=

Part II: Hazards in Applying Newton's Method

Let's select a polynomial that we can't factor and begin by looking at its graph.

In[26]:=

We set up Newton's iteration scheme and specify 20-digit precision.

In[29]:=

This function has three zeros in the domain over which the plot extends. What happens when we specify a starting value of 2?

We iterate until two consecutive calculated values of x are equal to the precision specified by digits. The symbol n is used to count the number of iterations required to achieve the specified precision. For an explanation of the following commands, click on the "About Mathematica" button below the cell.

In[32]:=

Out[41]//TableForm=

 n x 0 2 1 3.25`20. 2 2.80789955214978214058250856778967691697`20. 3 2.53718280431355258049702600951693577511`20. 4 2.42174162330835861020965218020243064071`20. 5 2.40083733782435009331664440660881704023`20. 6 2.40019741845712916528138454326242828143`20. 7 2.40019683087362685607455054105314563877`20. 8 2.4001968308731317853141320863838139015`20. 9 2.40019683087313178527317722443790160014`20.

What if we had started at x=1.5?

In[43]:=

Out[52]//TableForm=

 n x 0 1.5` 1 0.61471861471861466430510745340143330395`20. 2 0.50283191084360848768813435217893920492`20. 3 0.48291779156501070408890623055014971927`20. 4 0.48228571347382140444169766416310027375`20. 5 0.48228508357646991850221795836066800277`20. 6 0.48228508357584461930303916760655870169`20. 7 0.48228508357584461930303916760520624444`20.

This time the iterations converged to a different value of x. Why did that happen? Can you tell precisely where the break would occur? Does it have anything to do with the derivative of the function?

What if we had started at x=0?

In[54]:=

Out[63]//TableForm=

 n x 0 0 1 ComplexInfinity 2 Indeterminate

This time Newton's method didn't work at all. What is it that has given the expression? It might help to look at the graph again.

In[65]:=

Let's try x = -1.5.

In[66]:=

Out[75]//TableForm=

 n x 0 -1.5` 1 0.87179487179487180625869768846314400434`20. 2 0.59433005145203672412763899507334330629`20. 3 0.49762232946285806650024903798440093934`20. 4 0.48264208414497544610791565970644126932`20. 5 0.4822852842537659608361355121020674189`20. 6 0.48228508357590808617557606435084849605`20. 7 0.48228508357584461930303855775154667099`20. 8 0.48228508357584461930303855774630674893`20.

Why does this converge to the zero to the right instead of to the closer zero to the left?

You Try It: Part II

Here is a cubic polynomial. Study its graph.

In[77]:=

You decide where you should begin your iterations in order to compute each of the zeros of this polynomial. By changing the term in red, you determine which starting values converge to which zero. Check it out with the given starting point of x = 0.

In[83]:=

Out[84]//TableForm=

 n x 0 0 1 5.5`20. 2 5.05185185185185185185185185185540464763`20. 3 4.95954124047054114462870755343416766942`20. 4 4.9556778219555245167016352431626727718`20. 5 4.95567115855749833026146705825486947813`20. 6 4.95567115853769190583509779841672475535`20. 7 4.9556711585376919057825857717336504038`20.

What happened? It seems to have completely skipped over a zero! Why did this happen?