


From: jypuppy36@hotmail.com
To: r-help-bounces@r-project.org
Subject: R programing help-newton iterations for the square root
Date: Tue, 7 Dec 2010 12:00:01 -0800








NEWTON ITERATIONS FOR THE SQUARE ROOT

Newton iterations to find the root of a real valued function f , i.e. a number x for which f (x) = 0, are of the form

 

Example. To find the square root of a positive number y we can use Newton’s

method to solve the equation f (x) = x^2 - y = 0. Since f '(x) = 2x we

see that

 

Global Convergence. Does Zangwill’s Theorem apply in this example ?

What follows from it ?

We can also use simple calculations here. Show that, if we start with x(0) >

0,



for all k > 0 and



for all k > 0. This implies the sequence decreases monotonically and converges

to from any (positive) starting point.

Speed of Convergence. Does Ostrowski’s Theorem apply in this example ?

What follows from it ?

We can also use simple calculation here to determine the speed of converge.

Show first that



which implies superlinear convergence, and then show



which means convergence is quadratic.

Coding.

(1) Write an R function newton() that applies Newton’s method to a 
general univariate polynomial equation f (x) = 0. Use the standard 
iterative setup with itmax, eps, verbose.

(2) Use the R package polynom to handle polynomials and their derivatives.

If it is not there, install it. You construct a polynomial p from a 
vector a as p <- polynomial(a) (try this out for various a !) and 
then you can use predict (p,x) to evaluate the polynomial at x and you 
can use predict (deriv(p) ,x) to evaluate the derivative at x.

(3) Thus the args of your function are newton ( p , x , eps , i tmax , v
 e r b o s e ) where p is the vector defining the polynomial, and x is 
the starting point for the iterations.

(4) What can you say about the limitations of your newton() function ?

When does it work, when doesn’t it work ? Is there an easy way to make it more reliable ?

(5) What happens for f (x) = x^3, for example ?

(6) Compare and check with the function solve .polynomial() in the 
polynom package, which computes all roots by algebraic methods.





How to use Polynom

>install.packages(“polynom”)

> library(polynom)

> ?polynomial

> polynomial(1:4)

1 + 2*x + 3*x^2 + 4*x^3

> p <- polynomial(c(-2,0,1))         # f( x )

> print(p, decreasing = TRUE)

x^2 - 2

> predict (p,3)                   # f ( 3 )

[1] 7

> predict (deriv(p) ,3)        # f ' ( 3 )

[1] 6

> solve(p) # gives the root

[1] -1.414214 1.414214



Basic function skeleton

MyFunctionName

                  <- function( p, x, itmax = 100, eps = 1e-6, verbose = FALSE)

{

          # INITIAL PROCESS

           repeat{

                        # MAIN PROCESS

                       if( verbose ) {

 # LOG PROCESS

                       }

                      if( ( CONDITION ) || ( itel == itmax ) ){

 return( RETURNING VALUE )

                       }

                      # SAVE OLD VALUE

                      itel <- itel + 1;

                       }

} 		 	   		  
	[[alternative HTML version deleted]]

