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]]
